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 // 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_CONSERVATIVESPARSESPARSEPRODUCT_H
11 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template<typename Lhs, typename Rhs, typename ResultType>
18 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
19 {
20  typedef typename remove_all<Lhs>::type::Scalar Scalar;
21  typedef typename remove_all<Lhs>::type::Index Index;
22 
23  // make sure to call innerSize/outerSize since we fake the storage order.
24  Index rows = lhs.innerSize();
25  Index cols = rhs.outerSize();
26  eigen_assert(lhs.outerSize() == rhs.innerSize());
27 
28  std::vector<bool> mask(rows,false);
29  Matrix<Scalar,Dynamic,1> values(rows);
30  Matrix<Index,Dynamic,1> indices(rows);
31 
32  // estimate the number of non zero entries
33  // given a rhs column containing Y non zeros, we assume that the respective Y columns
34  // of the lhs differs in average of one non zeros, thus the number of non zeros for
35  // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
36  // per column of the lhs.
37  // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
38  Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
39 
40  res.setZero();
41  res.reserve(Index(estimated_nnz_prod));
42  // we compute each column of the result, one after the other
43  for (Index j=0; j<cols; ++j)
44  {
45 
46  res.startVec(j);
47  Index nnz = 0;
48  for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
49  {
50  Scalar y = rhsIt.value();
51  Index k = rhsIt.index();
52  for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
53  {
54  Index i = lhsIt.index();
55  Scalar x = lhsIt.value();
56  if(!mask[i])
57  {
58  mask[i] = true;
59  values[i] = x * y;
60  indices[nnz] = i;
61  ++nnz;
62  }
63  else
64  values[i] += x * y;
65  }
66  }
67 
68  // unordered insertion
69  for(int k=0; k<nnz; ++k)
70  {
71  int i = indices[k];
72  res.insertBackByOuterInnerUnordered(j,i) = values[i];
73  mask[i] = false;
74  }
75 
76 #if 0
77  // alternative ordered insertion code:
78 
79  int t200 = rows/(log2(200)*1.39);
80  int t = (rows*100)/139;
81 
82  // FIXME reserve nnz non zeros
83  // FIXME implement fast sort algorithms for very small nnz
84  // if the result is sparse enough => use a quick sort
85  // otherwise => loop through the entire vector
86  // In order to avoid to perform an expensive log2 when the
87  // result is clearly very sparse we use a linear bound up to 200.
88  //if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
89  //res.startVec(j);
90  if(true)
91  {
92  if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
93  for(int k=0; k<nnz; ++k)
94  {
95  int i = indices[k];
96  res.insertBackByOuterInner(j,i) = values[i];
97  mask[i] = false;
98  }
99  }
100  else
101  {
102  // dense path
103  for(int i=0; i<rows; ++i)
104  {
105  if(mask[i])
106  {
107  mask[i] = false;
108  res.insertBackByOuterInner(j,i) = values[i];
109  }
110  }
111  }
112 #endif
113 
114  }
115  res.finalize();
116 }
117 
118 
119 } // end namespace internal
120 
121 namespace internal {
122 
123 template<typename Lhs, typename Rhs, typename ResultType,
124  int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
125  int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
126  int ResStorageOrder = (traits<ResultType>::Flags&RowMajorBit) ? RowMajor : ColMajor>
128 
129 template<typename Lhs, typename Rhs, typename ResultType>
131 {
133  typedef typename LhsCleaned::Scalar Scalar;
134 
135  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
136  {
139  ColMajorMatrix resCol(lhs.rows(),rhs.cols());
140  internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
141  // sort the non zeros:
142  RowMajorMatrix resRow(resCol);
143  res = resRow;
144  }
145 };
146 
147 template<typename Lhs, typename Rhs, typename ResultType>
149 {
150  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
151  {
153  RowMajorMatrix rhsRow = rhs;
154  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
155  internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
156  res = resRow;
157  }
158 };
159 
160 template<typename Lhs, typename Rhs, typename ResultType>
162 {
163  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
164  {
166  RowMajorMatrix lhsRow = lhs;
167  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
168  internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
169  res = resRow;
170  }
171 };
172 
173 template<typename Lhs, typename Rhs, typename ResultType>
175 {
176  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
177  {
179  RowMajorMatrix resRow(lhs.rows(), rhs.cols());
180  internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
181  res = resRow;
182  }
183 };
184 
185 
186 template<typename Lhs, typename Rhs, typename ResultType>
188 {
190 
191  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
192  {
194  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
195  internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
196  res = resCol;
197  }
198 };
199 
200 template<typename Lhs, typename Rhs, typename ResultType>
202 {
203  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
204  {
206  ColMajorMatrix lhsCol = lhs;
207  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
208  internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
209  res = resCol;
210  }
211 };
212 
213 template<typename Lhs, typename Rhs, typename ResultType>
215 {
216  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
217  {
219  ColMajorMatrix rhsCol = rhs;
220  ColMajorMatrix resCol(lhs.rows(), rhs.cols());
221  internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
222  res = resCol;
223  }
224 };
225 
226 template<typename Lhs, typename Rhs, typename ResultType>
228 {
229  static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
230  {
233  RowMajorMatrix resRow(lhs.rows(),rhs.cols());
234  internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
235  // sort the non zeros:
236  ColMajorMatrix resCol(resRow);
237  res = resCol;
238  }
239 };
240 
241 } // end namespace internal
242 
243 } // end namespace Eigen
244 
245 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
static void conservative_sparse_sparse_product_impl(const Lhs &lhs, const Rhs &rhs, ResultType &res)
A versatible sparse matrix representation.
Definition: SparseMatrix.h:85
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
const unsigned int RowMajorBit
void rhs(const real_t *x, real_t *f)
EIGEN_STRONG_INLINE const Scalar * data() const
#define eigen_assert(x)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:31