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

Use std::vector instead of DiagonalMatrix.

Browse Source 
This is because DiagonalMatrix lacks a swap() method.
This commit is contained in:
Atgeirr Flø Rasmussen
2015-08-27 09:33:30 +02:00
committed by babrodtk
parent 87f677af02
commit 3c905845f9
2 changed files with 448 additions and 433 deletions
Download Patch File Download Diff File Expand all files Collapse all files

871
opm/autodiff/AutoDiffMatrix.hpp
View File

@@ -28,6 +28,8 @@
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <opm/autodiff/fastSparseProduct.hpp> #include <opm/autodiff/fastSparseProduct.hpp>
#include <vector>
namespace Opm namespace Opm
{ {
@@ -35,434 +37,445 @@ namespace Opm
class AutoDiffMatrix class AutoDiffMatrix
{ {
public: public:
AutoDiffMatrix() AutoDiffMatrix()
: type_(Z), : type_(Z),
rows_(0), rows_(0),
cols_(0) cols_(0)
{ {
} }
AutoDiffMatrix(const int rows, const int cols) AutoDiffMatrix(const int rows, const int cols)
: type_(Z), : type_(Z),
rows_(rows), rows_(rows),
cols_(cols) cols_(cols)
{ {
} }
enum CreationType { ZeroMatrix, IdentityMatrix }; enum CreationType { ZeroMatrix, IdentityMatrix };
AutoDiffMatrix(const CreationType t, const int rows) AutoDiffMatrix(const CreationType t, const int rows)
: type_(t == ZeroMatrix ? Z : I), : type_(t == ZeroMatrix ? Z : I),
rows_(rows), rows_(rows),
cols_(rows) cols_(rows)
{ {
} }
explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d) explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d)
: type_(D), : type_(D),
rows_(d.rows()), rows_(d.rows()),
cols_(d.cols()), cols_(d.cols()),
d_(d) d_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows())
{ {
} }
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()),
s_(s) s_(s)
{ {
} }
AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const void swap(AutoDiffMatrix& other)
{ {
assert(rows_ == rhs.rows_); std::swap(type_, other.type_);
assert(cols_ == rhs.cols_); std::swap(rows_, other.rows_);
switch (type_) { std::swap(cols_, other.cols_);
case Z: d_.swap(other.d_);
return rhs; s_.swap(other.s_);
case I: }
switch (rhs.type_) {
case Z:
return *this;
case I: AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const
return sumII(*this, rhs); {
case D: assert(rows_ == rhs.rows_);
return rhs + (*this); assert(cols_ == rhs.cols_);
case S: switch (type_) {
return rhs + (*this); case Z:
} return rhs;
case D: case I:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Z:
return *this; return *this;
case I: case I:
return sumDI(*this, rhs); return sumII(*this, rhs);
case D: case D:
return sumDD(*this, rhs); return rhs + (*this);
case S: case S:
return rhs + (*this); return rhs + (*this);
} }
case S: case D:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Z:
return *this; return *this;
case I: case I:
return sumSI(*this, rhs); return sumDI(*this, rhs);
case D: case D:
return sumSD(*this, rhs); return sumDD(*this, rhs);
case S: case S:
return sumSS(*this, rhs); return rhs + (*this);
} }
} case S:
} switch (rhs.type_) {
case Z:
AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const return *this;
{ case I:
assert(cols_ == rhs.rows_); return sumSI(*this, rhs);
switch (type_) { case D:
case Z: return sumSD(*this, rhs);
return AutoDiffMatrix(rows_, rhs.cols_); case S:
case I: return sumSS(*this, rhs);
switch (rhs.type_) { }
case Z: }
return rhs; }
case I:
return rhs; AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const
case D: {
return rhs; assert(cols_ == rhs.rows_);
case S: switch (type_) {
return rhs; case Z:
} return AutoDiffMatrix(rows_, rhs.cols_);
case D: case I:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Z:
return AutoDiffMatrix(rows_, rhs.cols_); return rhs;
case I: case I:
return *this; return rhs;
case D: case D:
return prodDD(*this, rhs); return rhs;
case S: case S:
return prodDS(*this, rhs); return rhs;
} }
case S: case D:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Z:
return AutoDiffMatrix(rows_, rhs.cols_); return AutoDiffMatrix(rows_, rhs.cols_);
case I: case I:
return *this; return *this;
case D: case D:
return prodSD(*this, rhs); return prodDD(*this, rhs);
case S: case S:
return prodSS(*this, rhs); return prodDS(*this, rhs);
} }
} case S:
} switch (rhs.type_) {
case Z:
return AutoDiffMatrix(rows_, rhs.cols_);
case I:
return *this;
case D:
return prodSD(*this, rhs);
case S:
AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs) return prodSS(*this, rhs);
{ }
*this = *this + rhs; }
return *this; }
}
AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs) AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs)
{ {
*this = *this + rhs * -1.0; *this = *this + rhs;
return *this; return *this;
} }
AutoDiffMatrix operator*(const double rhs) const AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs)
{ {
switch (type_) { *this = *this + rhs * -1.0;
case Z: return *this;
return *this; }
case I:
{
AutoDiffMatrix retval(*this);
retval.type_ = D;
retval.d_ = rhs * Eigen::VectorXd::Ones(rows_).asDiagonal();
return retval;
} AutoDiffMatrix operator*(const double rhs) const
case D: {
{ switch (type_) {
AutoDiffMatrix retval(*this); case Z:
retval.d_ = rhs * retval.d_; return *this;
return retval; case I:
} {
case S: AutoDiffMatrix retval(*this);
{ retval.type_ = D;
AutoDiffMatrix retval(*this); retval.d_.assign(rows_, rhs);
retval.s_ *= rhs; return retval;
return retval; }
} case D:
} {
} AutoDiffMatrix retval(*this);
for (double& elem : retval.d_) {
elem *= rhs;
}
return retval;
}
case S:
Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const {
{ AutoDiffMatrix retval(*this);
assert(cols_ == rhs.size()); retval.s_ *= rhs;
switch (type_) { return retval;
case Z: }
return Eigen::VectorXd::Zero(rows_); }
case I: }
return rhs;
case D:
return d_ * rhs;
case S:
return s_ * rhs;
}
} Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const
{
assert(cols_ == rhs.size());
switch (type_) {
case Z:
return Eigen::VectorXd::Zero(rows_);
case I:
static AutoDiffMatrix sumII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) return rhs;
{ case D:
assert(lhs.type_ == I); return Eigen::Map<const Eigen::VectorXd>(d_.data(), rows_) * rhs;
assert(rhs.type_ == I); case S:
AutoDiffMatrix retval; return s_ * rhs;
retval.type_ = D; }
retval.rows_ = lhs.rows_; }
retval.cols_ = rhs.cols_;
retval.d_ = Eigen::VectorXd::Constant(lhs.rows_, 2.0).asDiagonal();
return retval;
}
static AutoDiffMatrix sumDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ static AutoDiffMatrix sumII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(lhs.type_ == D); {
assert(rhs.type_ == I); assert(lhs.type_ == I);
AutoDiffMatrix retval = lhs; assert(rhs.type_ == I);
for (int r = 0; r < lhs.rows_; ++r) { AutoDiffMatrix retval;
retval.d_.diagonal()(r) += 1.0; retval.type_ = D;
} retval.rows_ = lhs.rows_;
return retval; retval.cols_ = rhs.cols_;
} retval.d_.assign(lhs.rows_, 2.0);
return retval;
static AutoDiffMatrix sumDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) }
{
assert(lhs.type_ == D); static AutoDiffMatrix sumDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(rhs.type_ == D); {
AutoDiffMatrix retval = lhs; static_cast<void>(rhs); // Silence release-mode warning.
for (int r = 0; r < lhs.rows_; ++r) { assert(lhs.type_ == D);
retval.d_.diagonal()(r) += rhs.d_.diagonal()(r); assert(rhs.type_ == I);
} AutoDiffMatrix retval = lhs;
return retval; for (int r = 0; r < lhs.rows_; ++r) {
} retval.d_[r] += 1.0;
}
static AutoDiffMatrix sumSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) return retval;
{ }
assert(lhs.type_ == S);
assert(rhs.type_ == I); static AutoDiffMatrix sumDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
AutoDiffMatrix retval; {
Eigen::SparseMatrix<double> ident = spdiag(Eigen::VectorXd::Ones(lhs.rows_)); assert(lhs.type_ == D);
retval.type_ = S; assert(rhs.type_ == D);
retval.rows_ = lhs.rows_; AutoDiffMatrix retval = lhs;
retval.cols_ = rhs.cols_; for (int r = 0; r < lhs.rows_; ++r) {
retval.s_ = lhs.s_ + ident; retval.d_[r] += rhs.d_[r];
return retval; }
} return retval;
}
static AutoDiffMatrix sumSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ static AutoDiffMatrix sumSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(lhs.type_ == S); {
assert(rhs.type_ == D); assert(lhs.type_ == S);
AutoDiffMatrix retval; assert(rhs.type_ == I);
Eigen::SparseMatrix<double> diag = spdiag(rhs.d_.diagonal()); AutoDiffMatrix retval;
retval.type_ = S; Eigen::SparseMatrix<double> ident = spdiag(Eigen::VectorXd::Ones(lhs.rows_));
retval.rows_ = lhs.rows_; retval.type_ = S;
retval.cols_ = rhs.cols_; retval.rows_ = lhs.rows_;
retval.s_ = lhs.s_ + diag; retval.cols_ = rhs.cols_;
return retval; retval.s_ = lhs.s_ + ident;
} return retval;
}
static AutoDiffMatrix sumSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ static AutoDiffMatrix sumSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(lhs.type_ == S); {
assert(rhs.type_ == S); assert(lhs.type_ == S);
AutoDiffMatrix retval; assert(rhs.type_ == D);
retval.type_ = S; AutoDiffMatrix retval;
retval.rows_ = lhs.rows_; Eigen::SparseMatrix<double> diag = spdiag(rhs.d_);
retval.cols_ = rhs.cols_; retval.type_ = S;
retval.s_ = lhs.s_ + rhs.s_; retval.rows_ = lhs.rows_;
return retval; retval.cols_ = rhs.cols_;
} retval.s_ = lhs.s_ + diag;
return retval;
}
static AutoDiffMatrix sumSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{
static AutoDiffMatrix prodDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) assert(lhs.type_ == S);
{ assert(rhs.type_ == S);
assert(lhs.type_ == D); AutoDiffMatrix retval;
assert(rhs.type_ == D); retval.type_ = S;
AutoDiffMatrix retval; retval.rows_ = lhs.rows_;
retval.type_ = D; retval.cols_ = rhs.cols_;
retval.rows_ = lhs.rows_; retval.s_ = lhs.s_ + rhs.s_;
retval.cols_ = rhs.cols_; return retval;
retval.d_ = (lhs.d_.diagonal().array() * rhs.d_.diagonal().array()).matrix().asDiagonal(); }
return retval;
}
static AutoDiffMatrix prodDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{
assert(lhs.type_ == D); static AutoDiffMatrix prodDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(rhs.type_ == S); {
AutoDiffMatrix retval; assert(lhs.type_ == D);
// Eigen::SparseMatrix<double> diag = spdiag(lhs.d_.diagonal()); assert(rhs.type_ == D);
retval.type_ = S; AutoDiffMatrix retval;
retval.rows_ = lhs.rows_; retval.type_ = D;
retval.cols_ = rhs.cols_; retval.rows_ = lhs.rows_;
// fastSparseProduct(diag, rhs.s_, retval.s_); retval.cols_ = rhs.cols_;
fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_); retval.d_.resize(lhs.rows_);
// retval.s_ = lhs.d_ * rhs.s_; for (int r = 0; r < lhs.rows_; ++r) {
return retval; retval.d_[r] = lhs.d_[r] * rhs.d_[r];
} }
return retval;
static AutoDiffMatrix prodSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) }
{
assert(lhs.type_ == S); static AutoDiffMatrix prodDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
assert(rhs.type_ == D); {
AutoDiffMatrix retval; assert(lhs.type_ == D);
// Eigen::SparseMatrix<double> diag = spdiag(rhs.d_.diagonal()); assert(rhs.type_ == S);
retval.type_ = S; AutoDiffMatrix retval;
retval.rows_ = lhs.rows_; retval.type_ = S;
retval.cols_ = rhs.cols_; retval.rows_ = lhs.rows_;
// fastSparseProduct(lhs.s_, diag, retval.s_); retval.cols_ = rhs.cols_;
fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_); fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_);
// retval.s_ = lhs.s_ * rhs.d_; return retval;
return retval; }
}
static AutoDiffMatrix prodSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
static AutoDiffMatrix prodSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) {
{ assert(lhs.type_ == S);
assert(lhs.type_ == S); assert(rhs.type_ == D);
assert(rhs.type_ == S); AutoDiffMatrix retval;
AutoDiffMatrix retval; retval.type_ = S;
retval.type_ = S; retval.rows_ = lhs.rows_;
retval.rows_ = lhs.rows_; retval.cols_ = rhs.cols_;
retval.cols_ = rhs.cols_; fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_);
fastSparseProduct(lhs.s_, rhs.s_, retval.s_); return retval;
return retval; }
}
static AutoDiffMatrix prodSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{
assert(lhs.type_ == S);
void toSparse(Eigen::SparseMatrix<double>& s) const assert(rhs.type_ == S);
{ AutoDiffMatrix retval;
switch (type_) { retval.type_ = S;
case Z: retval.rows_ = lhs.rows_;
s = Eigen::SparseMatrix<double>(rows_, cols_); retval.cols_ = rhs.cols_;
return; fastSparseProduct(lhs.s_, rhs.s_, retval.s_);
case I: return retval;
s = spdiag(Eigen::VectorXd::Ones(rows_)); }
return;
case D:
s = spdiag(d_.diagonal());
return; void toSparse(Eigen::SparseMatrix<double>& s) const
case S: {
s = s_; switch (type_) {
return; case Z:
} s = Eigen::SparseMatrix<double>(rows_, cols_);
} return;
case I:
s = spdiag(Eigen::VectorXd::Ones(rows_));
int rows() const return;
{ case D:
return rows_; s = spdiag(d_);
} return;
case S:
int cols() const s = s_;
{ return;
return cols_; }
} }
int nonZeros() const
{ int rows() const
switch (type_) { {
case Z: return rows_;
return 0; }
case I:
return rows_; int cols() const
case D: {
return rows_; return cols_;
case S: }
return s_.nonZeros();
} int nonZeros() const
} {
switch (type_) {
case Z:
double coeff(const int row, const int col) const return 0;
{ case I:
switch (type_) { return rows_;
case Z: case D:
return 0.0; return rows_;
case I: case S:
return (row == col) ? 1.0 : 0.0; return s_.nonZeros();
case D: }
return (row == col) ? d_.diagonal()(row) : 0.0; }
case S:
return s_.coeff(row, col);
} double coeff(const int row, const int col) const
} {
switch (type_) {
case Z:
return 0.0;
case I:
return (row == col) ? 1.0 : 0.0;
case D:
return (row == col) ? d_[row] : 0.0;
case S:
return s_.coeff(row, col);
}
}
private: private:
enum MatrixType { Z, I, D, S }; enum MatrixType { Z, I, D, S };
MatrixType type_; MatrixType type_;
int rows_; int rows_;
int cols_; int cols_;
Eigen::DiagonalMatrix<double, Eigen::Dynamic> d_; std::vector<double> d_;
Eigen::SparseMatrix<double> s_; Eigen::SparseMatrix<double> s_;
template <class V> template <class V>
static inline static inline
Eigen::SparseMatrix<double> Eigen::SparseMatrix<double>
spdiag(const V& d) spdiag(const V& d)
{ {
typedef Eigen::SparseMatrix<double> M; typedef Eigen::SparseMatrix<double> M;
const int n = d.size(); const int n = d.size();
M mat(n, n); M 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 (M::Index i = 0; i < n; ++i) {
mat.insert(i, i) = d[i]; mat.insert(i, i) = d[i];
} }
return mat; return mat;
} }
}; };
@@ -471,21 +484,21 @@ namespace Opm
inline void fastSparseProduct(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res) inline void fastSparseProduct(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
{ {
res = lhs * rhs; res = lhs * rhs;
} }
inline void fastSparseProduct(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res) inline void fastSparseProduct(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
{ {
res = AutoDiffMatrix(lhs) * rhs; res = AutoDiffMatrix(lhs) * rhs;
} }
inline AutoDiffMatrix operator*(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs) inline AutoDiffMatrix operator*(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs)
{ {
AutoDiffMatrix retval; AutoDiffMatrix retval;
fastSparseProduct(lhs, rhs, retval); fastSparseProduct(lhs, rhs, retval);
return retval; return retval;
} }
} // namespace Opm } // namespace Opm

10
opm/autodiff/fastSparseProduct.hpp
View File

@@ -142,7 +142,8 @@ void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res)
inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& lhs, inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& lhs,
const std::vector<double>& lhs,
const Eigen::SparseMatrix<double>& rhs, const Eigen::SparseMatrix<double>& rhs,
Eigen::SparseMatrix<double>& res) Eigen::SparseMatrix<double>& res)
{ {
@@ -153,7 +154,7 @@ inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix<double, Eigen::Dyn
for (int col = 0; col < n; ++col) { for (int col = 0; col < n; ++col) {
typedef Eigen::SparseMatrix<double>::InnerIterator It; typedef Eigen::SparseMatrix<double>::InnerIterator It;
for (It it(res, col); it; ++it) { for (It it(res, col); it; ++it) {
it.valueRef() *= lhs.diagonal()(it.row()); it.valueRef() *= lhs[it.row()]; // lhs.diagonal()(it.row());
} }
} }
} }
@@ -162,7 +163,8 @@ inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix<double, Eigen::Dyn
inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs, inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& rhs, // const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& rhs,
const std::vector<double>& rhs,
Eigen::SparseMatrix<double>& res) Eigen::SparseMatrix<double>& res)
{ {
res = lhs; res = lhs;
@@ -172,7 +174,7 @@ inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
for (int col = 0; col < n; ++col) { for (int col = 0; col < n; ++col) {
typedef Eigen::SparseMatrix<double>::InnerIterator It; typedef Eigen::SparseMatrix<double>::InnerIterator It;
for (It it(res, col); it; ++it) { for (It it(res, col); it; ++it) {
it.valueRef() *= rhs.diagonal()(col); it.valueRef() *= rhs[col]; // rhs.diagonal()(col);
} }
} }
} }