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

Implemented statically allocated memory for SparseMatrix

Browse Source
This commit is contained in:
babrodtk 2015-08-28 14:03:52 +02:00
parent df1d0f795d
commit 0b1f993588
1 changed files with 67 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -24,7 +24,6 @@
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <Eigen/Sparse> #include <Eigen/Sparse>
#include <boost/any.hpp>
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
@ -76,7 +75,7 @@ namespace Opm
: type_(D), : type_(D),
rows_(d.rows()), rows_(d.rows()),
cols_(d.cols()), cols_(d.cols()),
data_(Diag(d.diagonal().array().data(), d.diagonal().array().data() + d.rows())) diag_(Diag(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()))
{ {
} }
@ -85,15 +84,34 @@ namespace Opm
explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s) explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s)
: type_(S), : type_(S),
rows_(s.rows()), rows_(s.rows()),
cols_(s.cols()), cols_(s.cols())
data_(s)
{ {
sparse_[0] = s;
} }
AutoDiffMatrix(const AutoDiffMatrix& other) = default; AutoDiffMatrix(const AutoDiffMatrix& other)
AutoDiffMatrix& operator=(const AutoDiffMatrix& other) = default; {
*this = other;
}
AutoDiffMatrix& operator=(const AutoDiffMatrix& other)
{
type_ = other.type_;
rows_ = other.rows_;
cols_ = other.cols_;
switch(type_) {
case D:
diag_ = other.diag_;
break;
case S:
sparse_[0] = other.sparse_[0];
break;
default:
break;
}
}
@ -118,7 +136,8 @@ namespace Opm
std::swap(type_, other.type_); std::swap(type_, other.type_);
std::swap(rows_, other.rows_); std::swap(rows_, other.rows_);
std::swap(cols_, other.cols_); std::swap(cols_, other.cols_);
data_.swap(other.data_); diag_.swap(other.diag_);
sparse_[0].swap(other.sparse_[0]);
} }
@ -261,14 +280,13 @@ namespace Opm
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = D;
retval.data_ = boost::any(Diag(rows_, rhs)); retval.diag_ = Diag(rows_, rhs);
return retval; return retval;
} }
case D: case D:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
Diag& d = boost::any_cast<Diag&>(retval.data_); for (double& elem : retval.diag_) {
for (double& elem : d) {
elem *= rhs; elem *= rhs;
} }
return retval; return retval;
@ -276,8 +294,7 @@ namespace Opm
case S: case S:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
Sparse& s = boost::any_cast<Sparse&>(retval.data_); retval.sparse_[0] *= rhs;
s *= rhs;
return retval; return retval;
} }
default: default:
@ -299,14 +316,13 @@ namespace Opm
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = D;
boost::any_cast<Diag&>(retval.data_).assign(rows_, 1.0/rhs); retval.diag_.assign(rows_, 1.0/rhs);
return retval; return retval;
} }
case D: case D:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
Diag& d = boost::any_cast<Diag&>(retval.data_); for (double& elem : retval.diag_) {
for (double& elem : d) {
elem /= rhs; elem /= rhs;
} }
return retval; return retval;
@ -314,8 +330,7 @@ namespace Opm
case S: case S:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
Sparse& s = boost::any_cast<Sparse&>(retval.data_); retval.sparse_[0] /= rhs;
s /= rhs;
return retval; return retval;
} }
default: default:
@ -338,13 +353,11 @@ namespace Opm
return rhs; return rhs;
case D: case D:
{ {
const Diag& d = boost::any_cast<const Diag&>(data_); return Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_) * rhs;
return Eigen::Map<const Eigen::VectorXd>(d.data(), rows_) * rhs;
} }
case S: case S:
{ {
const Sparse& s = boost::any_cast<const Sparse&>(data_); return sparse_[0] * rhs;
return s * rhs;
} }
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
@ -364,7 +377,7 @@ namespace Opm
retval.type_ = D; retval.type_ = D;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.data_ = boost::any(Diag(lhs.rows_, 2.0)); retval.diag_.assign(lhs.rows_, 2.0);
return retval; return retval;
} }
@ -374,9 +387,8 @@ namespace Opm
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == I); assert(rhs.type_ == I);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
Diag& d = boost::any_cast<Diag&>(retval.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
d[r] += 1.0; retval.diag_[r] += 1.0;
} }
return retval; return retval;
} }
@ -386,10 +398,8 @@ namespace Opm
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
Diag& d_lhs = boost::any_cast<Diag&>(retval.data_);
const Diag& d_rhs = boost::any_cast<const Diag&>(rhs.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
d_lhs[r] += d_rhs[r]; retval.diag_[r] += rhs.diag_[r];
} }
return retval; return retval;
} }
@ -399,12 +409,11 @@ namespace Opm
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == I); assert(rhs.type_ == I);
AutoDiffMatrix retval; AutoDiffMatrix retval;
Eigen::SparseMatrix<double> ident = spdiag(Eigen::VectorXd::Ones(lhs.rows_));
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
const Sparse& s = boost::any_cast<const Sparse&>(lhs.data_); retval.sparse_[0] = lhs.sparse_[0];
retval.data_ = static_cast<Sparse>(s + ident); retval.sparse_[0] += spdiag(Eigen::VectorXd::Ones(lhs.rows_));
return retval; return retval;
} }
@ -413,12 +422,11 @@ namespace Opm
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval; AutoDiffMatrix retval;
Sparse diag = spdiag(boost::any_cast<const Diag&>(rhs.data_));
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
const Sparse& s = boost::any_cast<const Sparse&>(lhs.data_); retval.sparse_[0] = lhs.sparse_[0];
retval.data_ = static_cast<Sparse>(s + diag); retval.sparse_[0] += spdiag(rhs.diag_);
return retval; return retval;
} }
@ -426,13 +434,8 @@ namespace Opm
{ {
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == S); assert(rhs.type_ == S);
AutoDiffMatrix retval; AutoDiffMatrix retval = lhs;
retval.type_ = S; retval.sparse_[0] += rhs.sparse_[0];
retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_;
const Sparse& s_lhs = boost::any_cast<const Sparse&>(lhs.data_);
const Sparse& s_rhs = boost::any_cast<const Sparse&>(rhs.data_);
retval.data_ = static_cast<Sparse>(s_lhs + s_rhs);
return retval; return retval;
} }
@ -444,16 +447,9 @@ namespace Opm
{ {
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval; AutoDiffMatrix retval = lhs;
retval.type_ = D;
retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_;
retval.data_ = boost::any(Diag(lhs.rows_));
Diag& d = boost::any_cast<Diag&>(retval.data_);
const Diag& d_lhs = boost::any_cast<const Diag&>(lhs.data_);
const Diag& d_rhs = boost::any_cast<const Diag&>(rhs.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
d[r] = d_lhs[r] * d_rhs[r]; retval.diag_[r] *= rhs.diag_[r];
} }
return retval; return retval;
} }
@ -466,11 +462,8 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous? retval.sparse_[0] = Sparse(retval.rows_, retval.cols_);
const Diag& a = boost::any_cast<const Diag&>(lhs.data_); fastDiagSparseProduct(lhs.diag_, rhs.sparse_[0], retval.sparse_[0]);
const Sparse& b = boost::any_cast<const Sparse&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastDiagSparseProduct(a, b, c);
return retval; return retval;
} }
@ -482,11 +475,8 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous? retval.sparse_[0] = Sparse(retval.rows_, retval.cols_);
const Sparse& a = boost::any_cast<const Sparse&>(lhs.data_); fastSparseDiagProduct(lhs.sparse_[0], rhs.diag_, retval.sparse_[0]);
const Diag& b = boost::any_cast<const Diag&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastSparseDiagProduct(a, b, c);
return retval; return retval;
} }
@ -498,11 +488,8 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous? retval.sparse_[0] = Sparse(retval.rows_, retval.cols_);
const Sparse& a = boost::any_cast<const Sparse&>(lhs.data_); fastSparseProduct(lhs.sparse_[0], rhs.sparse_[0], retval.sparse_[0]);
const Sparse& b = boost::any_cast<const Sparse&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastSparseProduct(a, b, c);
return retval; return retval;
} }
@ -518,10 +505,10 @@ namespace Opm
s = spdiag(Eigen::VectorXd::Ones(rows_)); s = spdiag(Eigen::VectorXd::Ones(rows_));
return; return;
case D: case D:
s = spdiag(boost::any_cast<const Diag&>(data_)); s = spdiag(diag_);
return; return;
case S: case S:
s = boost::any_cast<const Sparse&>(data_); s = sparse_[0];
return; return;
} }
} }
@ -547,7 +534,7 @@ namespace Opm
case D: case D:
return rows_; return rows_;
case S: case S:
return boost::any_cast<const Sparse&>(data_).nonZeros(); return sparse_->nonZeros();
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
@ -562,9 +549,9 @@ namespace Opm
case I: case I:
return (row == col) ? 1.0 : 0.0; return (row == col) ? 1.0 : 0.0;
case D: case D:
return (row == col) ? boost::any_cast<const Diag&>(data_)[row] : 0.0; return (row == col) ? diag_[row] : 0.0;
case S: case S:
return boost::any_cast<const Sparse&>(data_).coeff(row, col); return sparse_->coeff(row, col);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
@ -575,22 +562,24 @@ namespace Opm
MatrixType type_; MatrixType type_;
int rows_; int rows_;
int cols_; int cols_;
boost::any data_; Diag diag_;
/*
std::vector<double> d_; /**
Eigen::SparseMatrix<double> s_; * Eigen uses memory allocation within the default constructor, so that
*/ * Sparse a; actually calls malloc. To prevent this, we here use
* Sparse a[1], and only construct the object when needed.
*/
Sparse sparse_[1];
template <class V> template <class V>
static inline static inline
Eigen::SparseMatrix<double> Sparse
spdiag(const V& d) spdiag(const V& d)
{ {
typedef Eigen::SparseMatrix<double> M;
const int n = d.size(); const int n = d.size();
M mat(n, n); Sparse mat(n, n);
mat.reserve(Eigen::ArrayXi::Ones(n, 1)); mat.reserve(Eigen::ArrayXi::Ones(n, 1));
for (M::Index i = 0; i < n; ++i) { for (Sparse::Index i = 0; i < n; ++i) {
if (d[i] != 0.0) { if (d[i] != 0.0) {
mat.insert(i, i) = d[i]; mat.insert(i, i) = d[i];
} }