Home | History | Annotate | Download | only in SparseCore 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_SPARSEPRODUCT_H
     11 #define EIGEN_SPARSEPRODUCT_H
     12 
     13 namespace Eigen {
     14 
     15 template<typename Lhs, typename Rhs>
     16 struct SparseSparseProductReturnType
     17 {
     18   typedef typename internal::traits<Lhs>::Scalar Scalar;
     19   enum {
     20     LhsRowMajor = internal::traits<Lhs>::Flags & RowMajorBit,
     21     RhsRowMajor = internal::traits<Rhs>::Flags & RowMajorBit,
     22     TransposeRhs = (!LhsRowMajor) && RhsRowMajor,
     23     TransposeLhs = LhsRowMajor && (!RhsRowMajor)
     24   };
     25 
     26   typedef typename internal::conditional<TransposeLhs,
     27     SparseMatrix<Scalar,0>,
     28     typename internal::nested<Lhs,Rhs::RowsAtCompileTime>::type>::type LhsNested;
     29 
     30   typedef typename internal::conditional<TransposeRhs,
     31     SparseMatrix<Scalar,0>,
     32     typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type>::type RhsNested;
     33 
     34   typedef SparseSparseProduct<LhsNested, RhsNested> Type;
     35 };
     36 
     37 namespace internal {
     38 template<typename LhsNested, typename RhsNested>
     39 struct traits<SparseSparseProduct<LhsNested, RhsNested> >
     40 {
     41   typedef MatrixXpr XprKind;
     42   // clean the nested types:
     43   typedef typename remove_all<LhsNested>::type _LhsNested;
     44   typedef typename remove_all<RhsNested>::type _RhsNested;
     45   typedef typename _LhsNested::Scalar Scalar;
     46   typedef typename promote_index_type<typename traits<_LhsNested>::Index,
     47                                          typename traits<_RhsNested>::Index>::type Index;
     48 
     49   enum {
     50     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
     51     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
     52     LhsFlags = _LhsNested::Flags,
     53     RhsFlags = _RhsNested::Flags,
     54 
     55     RowsAtCompileTime    = _LhsNested::RowsAtCompileTime,
     56     ColsAtCompileTime    = _RhsNested::ColsAtCompileTime,
     57     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
     58     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
     59 
     60     InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
     61 
     62     EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
     63 
     64     RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
     65 
     66     Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
     67           | EvalBeforeAssigningBit
     68           | EvalBeforeNestingBit,
     69 
     70     CoeffReadCost = Dynamic
     71   };
     72 
     73   typedef Sparse StorageKind;
     74 };
     75 
     76 } // end namespace internal
     77 
     78 template<typename LhsNested, typename RhsNested>
     79 class SparseSparseProduct : internal::no_assignment_operator,
     80   public SparseMatrixBase<SparseSparseProduct<LhsNested, RhsNested> >
     81 {
     82   public:
     83 
     84     typedef SparseMatrixBase<SparseSparseProduct> Base;
     85     EIGEN_DENSE_PUBLIC_INTERFACE(SparseSparseProduct)
     86 
     87   private:
     88 
     89     typedef typename internal::traits<SparseSparseProduct>::_LhsNested _LhsNested;
     90     typedef typename internal::traits<SparseSparseProduct>::_RhsNested _RhsNested;
     91 
     92   public:
     93 
     94     template<typename Lhs, typename Rhs>
     95     EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs)
     96       : m_lhs(lhs), m_rhs(rhs), m_tolerance(0), m_conservative(true)
     97     {
     98       init();
     99     }
    100 
    101     template<typename Lhs, typename Rhs>
    102     EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs, RealScalar tolerance)
    103       : m_lhs(lhs), m_rhs(rhs), m_tolerance(tolerance), m_conservative(false)
    104     {
    105       init();
    106     }
    107 
    108     SparseSparseProduct pruned(Scalar reference = 0, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) const
    109     {
    110       return SparseSparseProduct(m_lhs,m_rhs,internal::abs(reference)*epsilon);
    111     }
    112 
    113     template<typename Dest>
    114     void evalTo(Dest& result) const
    115     {
    116       if(m_conservative)
    117         internal::conservative_sparse_sparse_product_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result);
    118       else
    119         internal::sparse_sparse_product_with_pruning_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result,m_tolerance);
    120     }
    121 
    122     EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); }
    123     EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); }
    124 
    125     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
    126     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
    127 
    128   protected:
    129     void init()
    130     {
    131       eigen_assert(m_lhs.cols() == m_rhs.rows());
    132 
    133       enum {
    134         ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic
    135                       || _RhsNested::RowsAtCompileTime==Dynamic
    136                       || int(_LhsNested::ColsAtCompileTime)==int(_RhsNested::RowsAtCompileTime),
    137         AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
    138         SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested,_RhsNested)
    139       };
    140       // note to the lost user:
    141       //    * for a dot product use: v1.dot(v2)
    142       //    * for a coeff-wise product use: v1.cwise()*v2
    143       EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
    144         INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
    145       EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
    146         INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
    147       EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
    148     }
    149 
    150     LhsNested m_lhs;
    151     RhsNested m_rhs;
    152     RealScalar m_tolerance;
    153     bool m_conservative;
    154 };
    155 
    156 // sparse = sparse * sparse
    157 template<typename Derived>
    158 template<typename Lhs, typename Rhs>
    159 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseSparseProduct<Lhs,Rhs>& product)
    160 {
    161   product.evalTo(derived());
    162   return derived();
    163 }
    164 
    165 /** \returns an expression of the product of two sparse matrices.
    166   * By default a conservative product preserving the symbolic non zeros is performed.
    167   * The automatic pruning of the small values can be achieved by calling the pruned() function
    168   * in which case a totally different product algorithm is employed:
    169   * \code
    170   * C = (A*B).pruned();             // supress numerical zeros (exact)
    171   * C = (A*B).pruned(ref);
    172   * C = (A*B).pruned(ref,epsilon);
    173   * \endcode
    174   * where \c ref is a meaningful non zero reference value.
    175   * */
    176 template<typename Derived>
    177 template<typename OtherDerived>
    178 inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
    179 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
    180 {
    181   return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
    182 }
    183 
    184 } // end namespace Eigen
    185 
    186 #endif // EIGEN_SPARSEPRODUCT_H
    187