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@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 |
---|