OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Experiment with fast diagonal-sparse products.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2015-08-26 15:31:48 +02:00 committed by babrodtk
parent 1f32594f79
commit 11a33b3017
2 changed files with 46 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -343,11 +343,13 @@ namespace Opm
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == S); assert(rhs.type_ == S);
AutoDiffMatrix retval; AutoDiffMatrix retval;
Eigen::SparseMatrix<double> diag = spdiag(lhs.d_.diagonal()); // Eigen::SparseMatrix<double> diag = spdiag(lhs.d_.diagonal());
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseProduct(diag, rhs.s_, retval.s_); // fastSparseProduct(diag, rhs.s_, retval.s_);
fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_);
// retval.s_ = lhs.d_ * rhs.s_;
return retval; return retval;
} }
@ -356,11 +358,13 @@ namespace Opm
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval; AutoDiffMatrix retval;
Eigen::SparseMatrix<double> diag = spdiag(rhs.d_.diagonal()); // Eigen::SparseMatrix<double> diag = spdiag(rhs.d_.diagonal());
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseProduct(lhs.s_, diag, retval.s_); // fastSparseProduct(lhs.s_, diag, retval.s_);
fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_);
// retval.s_ = lhs.s_ * rhs.d_;
return retval; return retval;
} }

38
opm/autodiff/fastSparseProduct.hpp
View File

@ -141,6 +141,44 @@ void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res)
inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& lhs,
const Eigen::SparseMatrix<double>& rhs,
Eigen::SparseMatrix<double>& res)
{
res = rhs;
// Multiply rows by diagonal lhs.
int n = res.cols();
for (int col = 0; col < n; ++col) {
typedef Eigen::SparseMatrix<double>::InnerIterator It;
for (It it(res, col); it; ++it) {
it.valueRef() *= rhs.diagonal()(it.row());
}
}
}
inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& rhs,
Eigen::SparseMatrix<double>& res)
{
res = lhs;
// Multiply columns by diagonal rhs.
int n = res.cols();
for (int col = 0; col < n; ++col) {
typedef Eigen::SparseMatrix<double>::InnerIterator It;
for (It it(res, col); it; ++it) {
it.valueRef() *= rhs.diagonal()(col);
}
}
}
} // end namespace Opm } // end namespace Opm
#endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED #endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED