All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
fastSparseOperations.hpp
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 // This file has been modified for use in the OPM project codebase.
11 
12 #ifndef OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
13 #define OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
14 
15 #include <Eigen/Sparse>
16 
17 #include <algorithm>
18 #include <iterator>
19 #include <functional>
20 #include <limits>
21 #include <vector>
22 
23 #include <Eigen/Core>
24 
25 namespace Opm {
26 
27 template < unsigned int depth >
28 struct QuickSort
29 {
30  template <typename T>
31  static inline void sort(T begin, T end)
32  {
33  if (begin != end)
34  {
35  T middle = std::partition (begin, end,
36  std::bind2nd(std::less<typename std::iterator_traits<T>::value_type>(), *begin)
37  );
38  QuickSort< depth-1 >::sort(begin, middle);
39 
40  // std::sort (max(begin + 1, middle), end);
41  T new_middle = begin;
42  QuickSort< depth-1 >::sort(++new_middle, end);
43  }
44  }
45 };
46 
47 template <>
48 struct QuickSort< 0 >
49 {
50  template <typename T>
51  static inline void sort(T begin, T end)
52  {
53  // fall back to standard insertion sort
54  std::sort( begin, end );
55  }
56 };
57 
58 
59 template<typename Lhs, typename Rhs, typename ResultType>
60 void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res)
61 {
62  // initialize result
63  res = ResultType(lhs.rows(), rhs.cols());
64 
65  // if one of the matrices does not contain non zero elements
66  // the result will only contain an empty matrix
67  if( lhs.nonZeros() == 0 || rhs.nonZeros() == 0 )
68  return;
69 
70  typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
71  typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
72 
73  // make sure to call innerSize/outerSize since we fake the storage order.
74  Index rows = lhs.innerSize();
75  Index cols = rhs.outerSize();
76  eigen_assert(lhs.outerSize() == rhs.innerSize());
77 
78  std::vector<bool> mask(rows,false);
79  Eigen::Matrix<Scalar,Eigen::Dynamic,1> values(rows);
80  Eigen::Matrix<Index, Eigen::Dynamic,1> indices(rows);
81 
82  // estimate the number of non zero entries
83  // given a rhs column containing Y non zeros, we assume that the respective Y columns
84  // of the lhs differs in average of one non zeros, thus the number of non zeros for
85  // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
86  // per column of the lhs.
87  // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
88  Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
89 
90  res.setZero();
91  res.reserve(Index(estimated_nnz_prod));
92 
93  //const Scalar epsilon = std::numeric_limits< Scalar >::epsilon();
94  const Scalar epsilon = 0.0;
95 
96  // we compute each column of the result, one after the other
97  for (Index j=0; j<cols; ++j)
98  {
99  Index nnz = 0;
100  for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
101  {
102  const Scalar y = rhsIt.value();
103  for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
104  {
105  const Scalar val = lhsIt.value() * y;
106  if( std::abs( val ) > epsilon )
107  {
108  const Index i = lhsIt.index();
109  if(!mask[i])
110  {
111  mask[i] = true;
112  values[i] = val;
113  indices[nnz] = i;
114  ++nnz;
115  }
116  else
117  values[i] += val;
118  }
119  }
120  }
121 
122  if( nnz > 1 )
123  {
124  // sort indices for sorted insertion to avoid later copying
125  QuickSort< 1 >::sort( indices.data(), indices.data()+nnz );
126  }
127 
128  res.startVec(j);
129  // ordered insertion
130  // still using insertBackByOuterInnerUnordered since we know what we are doing
131  for(Index k=0; k<nnz; ++k)
132  {
133  const Index i = indices[k];
134  res.insertBackByOuterInnerUnordered(j,i) = values[i];
135  mask[i] = false;
136  }
137 
138  }
139  res.finalize();
140 }
141 
142 
143 
144 
145 inline void fastDiagSparseProduct(const std::vector<double>& lhs,
146  const Eigen::SparseMatrix<double>& rhs,
147  Eigen::SparseMatrix<double>& res)
148 {
149  res = rhs;
150 
151  // Multiply rows by diagonal lhs.
152  int n = res.cols();
153  for (int col = 0; col < n; ++col) {
154  typedef Eigen::SparseMatrix<double>::InnerIterator It;
155  for (It it(res, col); it; ++it) {
156  it.valueRef() *= lhs[it.row()];
157  }
158  }
159 }
160 
161 
162 
163 
164 inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
165  const std::vector<double>& rhs,
166  Eigen::SparseMatrix<double>& res)
167 {
168  res = lhs;
169 
170  // Multiply columns by diagonal rhs.
171  int n = res.cols();
172  for (int col = 0; col < n; ++col) {
173  typedef Eigen::SparseMatrix<double>::InnerIterator It;
174  for (It it(res, col); it; ++it) {
175  it.valueRef() *= rhs[col];
176  }
177  }
178 }
179 
180 template<typename Lhs, typename Rhs>
181 inline bool
182 equalSparsityPattern(const Lhs& lhs, const Rhs& rhs)
183 {
184  // if both matrices have equal storage and non zeros match, we can check sparsity pattern
185  bool equal = (Lhs::IsRowMajor == Rhs::IsRowMajor) && (lhs.nonZeros() == rhs.nonZeros());
186 
187  // check complete sparsity pattern
188  if( equal )
189  {
190  typedef std::size_t Index;
191  const Index outerSize = lhs.outerSize();
192  const Index rhsOuterSize = rhs.outerSize();
193  if( outerSize != rhsOuterSize )
194  {
195  return false;
196  }
197 
198  // outer indices
199  const auto rhsOuter = rhs.outerIndexPtr();
200  const auto lhsOuter = lhs.outerIndexPtr();
201  for(Index i=0; i<=outerSize; ++i )
202  {
203  if( lhsOuter[ i ] != rhsOuter[ i ] ) {
204  return false ;
205  }
206  }
207 
208  // inner indices
209  const auto rhsInner = rhs.innerIndexPtr();
210  const auto lhsInner = lhs.innerIndexPtr();
211 
212  const Index nnz = lhs.nonZeros();
213  for( Index i=0; i<nnz; ++i)
214  {
215  if( lhsInner[ i ] != rhsInner[ i ] ) {
216  return false;
217  }
218  }
219  }
220 
221  return equal;
222 }
223 
224 // this function substracts two sparse matrices
225 // if the sparsity pattern is the same a faster add/substract is performed
226 template<typename Lhs, typename Rhs>
227 inline void
228 fastSparseAdd(Lhs& lhs, const Rhs& rhs)
229 {
230  if( equalSparsityPattern( lhs, rhs ) )
231  {
232  typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
233  typedef std::size_t Index;
234 
235  const Index nnz = lhs.nonZeros();
236 
237  // fast add using only the data pointers
238  const Scalar* rhsV = rhs.valuePtr();
239  Scalar* lhsV = lhs.valuePtr();
240 
241  for(Index i=0; i<nnz; ++i )
242  {
243  lhsV[ i ] += rhsV[ i ];
244  }
245  }
246  else
247  {
248  // default Eigen operator+=
249  lhs += rhs;
250  }
251 }
252 
253 // this function substracts two sparse matrices
254 // if the sparsity pattern is the same a faster add/substract is performed
255 template<typename Lhs, typename Rhs>
256 inline void
257 fastSparseSubstract(Lhs& lhs, const Rhs& rhs)
258 {
259  if( equalSparsityPattern( lhs, rhs ) )
260  {
261  typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
262  typedef std::size_t Index;
263 
264  const Index nnz = lhs.nonZeros();
265 
266  // fast add using only the data pointers
267  const Scalar* rhsV = rhs.valuePtr();
268  Scalar* lhsV = lhs.valuePtr();
269 
270  for(Index i=0; i<nnz; ++i )
271  {
272  lhsV[ i ] -= rhsV[ i ];
273  }
274  }
275  else
276  {
277  // default Eigen operator-=
278  lhs -= rhs;
279  }
280 }
281 
282 } // end namespace Opm
283 
284 #endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
Definition: fastSparseOperations.hpp:28
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708