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();
131 for(Index k=0; k<nnz; ++k)
133 const Index i = indices[k];
134 res.insertBackByOuterInnerUnordered(j,i) = values[i];
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()];
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];
184 #endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
Definition: AdditionalObjectDeleter.hpp:22
void fastSparseDiagProduct(const Eigen::SparseMatrix< double > &lhs, const std::vector< double > &rhs, Eigen::SparseMatrix< double > &res)
Definition: fastSparseProduct.hpp:164
Definition: fastSparseProduct.hpp:28
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Definition: AutoDiffMatrix.hpp:696
static void sort(T begin, T end)
Definition: fastSparseProduct.hpp:51
static void sort(T begin, T end)
Definition: fastSparseProduct.hpp:31
void fastDiagSparseProduct(const std::vector< double > &lhs, const Eigen::SparseMatrix< double > &rhs, Eigen::SparseMatrix< double > &res)
Definition: fastSparseProduct.hpp:145