12 #ifndef OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED 13 #define OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED 15 #include <Eigen/Sparse> 27 template <
unsigned int depth >
31 static inline void sort(T begin, T end)
35 T middle = std::partition (begin, end,
36 std::bind2nd(std::less<
typename std::iterator_traits<T>::value_type>(), *begin)
51 static inline void sort(T begin, T end)
54 std::sort( begin, end );
59 template<
typename Lhs,
typename Rhs,
typename ResultType>
63 res = ResultType(lhs.rows(), rhs.cols());
67 if( lhs.nonZeros() == 0 || rhs.nonZeros() == 0 )
70 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
71 typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
74 Index rows = lhs.innerSize();
75 Index cols = rhs.outerSize();
76 eigen_assert(lhs.outerSize() == rhs.innerSize());
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);
88 Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
91 res.reserve(Index(estimated_nnz_prod));
94 const Scalar epsilon = 0.0;
97 for (Index j=0; j<cols; ++j)
100 for (
typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
102 const Scalar y = rhsIt.value();
103 for (
typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
105 const Scalar val = lhsIt.value() * y;
106 if( std::abs( val ) > epsilon )
108 const Index i = lhsIt.index();
125 QuickSort< 1 >::sort( indices.data(), indices.data()+nnz );
131 for(Index k=0; k<nnz; ++k)
133 const Index i = indices[k];
134 res.insertBackByOuterInnerUnordered(j,i) = values[i];
145 inline void fastDiagSparseProduct(
const std::vector<double>& lhs,
146 const Eigen::SparseMatrix<double>& rhs,
147 Eigen::SparseMatrix<double>& res)
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()];
164 inline void fastSparseDiagProduct(
const Eigen::SparseMatrix<double>& lhs,
165 const std::vector<double>& rhs,
166 Eigen::SparseMatrix<double>& res)
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];
180 template<
typename Lhs,
typename Rhs>
182 equalSparsityPattern(
const Lhs& lhs,
const Rhs& rhs)
185 bool equal = (Lhs::IsRowMajor == Rhs::IsRowMajor) && (lhs.nonZeros() == rhs.nonZeros());
190 typedef std::size_t Index;
191 const Index outerSize = lhs.outerSize();
192 const Index rhsOuterSize = rhs.outerSize();
193 if( outerSize != rhsOuterSize )
199 const auto rhsOuter = rhs.outerIndexPtr();
200 const auto lhsOuter = lhs.outerIndexPtr();
201 for(Index i=0; i<=outerSize; ++i )
203 if( lhsOuter[ i ] != rhsOuter[ i ] ) {
209 const auto rhsInner = rhs.innerIndexPtr();
210 const auto lhsInner = lhs.innerIndexPtr();
212 const Index nnz = lhs.nonZeros();
213 for( Index i=0; i<nnz; ++i)
215 if( lhsInner[ i ] != rhsInner[ i ] ) {
226 template<
typename Lhs,
typename Rhs>
228 fastSparseAdd(Lhs& lhs,
const Rhs& rhs)
230 if( equalSparsityPattern( lhs, rhs ) )
232 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
233 typedef std::size_t Index;
235 const Index nnz = lhs.nonZeros();
238 const Scalar* rhsV = rhs.valuePtr();
239 Scalar* lhsV = lhs.valuePtr();
241 for(Index i=0; i<nnz; ++i )
243 lhsV[ i ] += rhsV[ i ];
255 template<
typename Lhs,
typename Rhs>
257 fastSparseSubstract(Lhs& lhs,
const Rhs& rhs)
259 if( equalSparsityPattern( lhs, rhs ) )
261 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
262 typedef std::size_t Index;
264 const Index nnz = lhs.nonZeros();
267 const Scalar* rhsV = rhs.valuePtr();
268 Scalar* lhsV = lhs.valuePtr();
270 for(Index i=0; i<nnz; ++i )
272 lhsV[ i ] -= rhsV[ i ];
284 #endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED Definition: fastSparseOperations.hpp:28
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708