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