ConservativeSparseSparseProduct.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
26 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
32 template<typename Lhs, typename Rhs, typename ResultType>
33 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
34 {
35  typedef typename remove_all<Lhs>::type::Scalar Scalar;
36  typedef typename remove_all<Lhs>::type::Index Index;
37 
38  // make sure to call innerSize/outerSize since we fake the storage order.
39  Index rows = lhs.innerSize();
40  Index cols = rhs.outerSize();
41  eigen_assert(lhs.outerSize() == rhs.innerSize());
42 
43  std::vector<bool> mask(rows,false);
44  Matrix<Scalar,Dynamic,1> values(rows);
45  Matrix<Index,Dynamic,1> indices(rows);
46 
47  // estimate the number of non zero entries
48  // given a rhs column containing Y non zeros, we assume that the respective Y columns
49  // of the lhs differs in average of one non zeros, thus the number of non zeros for
50  // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
51  // per column of the lhs.
52  // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
53  Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
54 
55  res.setZero();
56  res.reserve(Index(estimated_nnz_prod));
57  // we compute each column of the result, one after the other
58  for (Index j=0; j<cols; ++j)
59  {
60 
61  res.startVec(j);
62  Index nnz = 0;
63  for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
64  {
65  Scalar y = rhsIt.value();
66  Index k = rhsIt.index();
67  for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
68  {
69  Index i = lhsIt.index();
70  Scalar x = lhsIt.value();
71  if(!mask[i])
72  {
73  mask[i] = true;
74  values[i] = x * y;
75  indices[nnz] = i;
76  ++nnz;
77  }
78  else
79  values[i] += x * y;
80  }
81  }
82 
83  // unordered insertion
84  for(int k=0; k<nnz; ++k)
85  {
86  int i = indices[k];
87  res.insertBackByOuterInnerUnordered(j,i) = values[i];
88  mask[i] = false;
89  }
90 
91 #if 0
92  // alternative ordered insertion code:
93 
94  int t200 = rows/(log2(200)*1.39);
95  int t = (rows*100)/139;
96 
97  // FIXME reserve nnz non zeros
98  // FIXME implement fast sort algorithms for very small nnz
99  // if the result is sparse enough => use a quick sort
100  // otherwise => loop through the entire vector
101  // In order to avoid to perform an expensive log2 when the
102  // result is clearly very sparse we use a linear bound up to 200.
103  //if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
104  //res.startVec(j);
105  if(true)
106  {
107  if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
108  for(int k=0; k<nnz; ++k)
109  {
110  int i = indices[k];
111  res.insertBackByOuterInner(j,i) = values[i];
112  mask[i] = false;
113  }
114  }
115  else
116  {
117  // dense path
118  for(int i=0; i<rows; ++i)
119  {
120  if(mask[i])
121  {
122  mask[i] = false;
123  res.insertBackByOuterInner(j,i) = values[i];
124  }
125  }
126  }
127 #endif
128 
129  }
130  res.finalize();
131 }
132 
133 
134 } // end namespace internal
135 
136 namespace internal {
137 
138 template<typename Lhs, typename Rhs, typename ResultType,
139  int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
140  int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
141  int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
142 struct conservative_sparse_sparse_product_selector;
143 
144 template<typename Lhs, typename Rhs, typename ResultType>
145 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
146 {
147  typedef typename remove_all<Lhs>::type LhsCleaned;
148  typedef typename LhsCleaned::Scalar Scalar;
149 
150  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
151  {
152  typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
153  typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
154  ColMajorMatrix resCol(lhs.rows(),rhs.cols());
155  internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
156  // sort the non zeros:
157  RowMajorMatrix resRow(resCol);
158  res = resRow;
159  }
160 };
161 
162 template<typename Lhs, typename Rhs, typename ResultType>
163 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
164 {
165  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
166  {
167  typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
168  RowMajorMatrix rhsRow = rhs;
169  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
170  internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
171  res = resRow;
172  }
173 };
174 
175 template<typename Lhs, typename Rhs, typename ResultType>
176 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
177 {
178  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
179  {
180  typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
181  RowMajorMatrix lhsRow = lhs;
182  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
183  internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
184  res = resRow;
185  }
186 };
187 
188 template<typename Lhs, typename Rhs, typename ResultType>
189 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
190 {
191  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
192  {
193  typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
194  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
195  internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
196  res = resRow;
197  }
198 };
199 
200 
201 template<typename Lhs, typename Rhs, typename ResultType>
202 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
203 {
204  typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
205 
206  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
207  {
208  typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
209  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
210  internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
211  res = resCol;
212  }
213 };
214 
215 template<typename Lhs, typename Rhs, typename ResultType>
216 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
217 {
218  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
219  {
220  typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
221  ColMajorMatrix lhsCol = lhs;
222  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
223  internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
224  res = resCol;
225  }
226 };
227 
228 template<typename Lhs, typename Rhs, typename ResultType>
229 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
230 {
231  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
232  {
233  typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
234  ColMajorMatrix rhsCol = rhs;
235  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
236  internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
237  res = resCol;
238  }
239 };
240 
241 template<typename Lhs, typename Rhs, typename ResultType>
242 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
243 {
244  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
245  {
246  typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
247  typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
248  RowMajorMatrix resRow(lhs.rows(),rhs.cols());
249  internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
250  // sort the non zeros:
251  ColMajorMatrix resCol(resRow);
252  res = resCol;
253  }
254 };
255 
256 } // end namespace internal
257 
258 } // end namespace Eigen
259 
260 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H