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