From 3c905845f9ef04157d3b9494c5d359d64723bf23 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Aug 2015 09:33:30 +0200 Subject: [PATCH] Use std::vector instead of DiagonalMatrix. This is because DiagonalMatrix lacks a swap() method. --- opm/autodiff/AutoDiffMatrix.hpp | 871 +++++++++++++++-------------- opm/autodiff/fastSparseProduct.hpp | 10 +- 2 files changed, 448 insertions(+), 433 deletions(-) diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index b86fd16dc..4cf491cd0 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -28,6 +28,8 @@ #include #include +#include + namespace Opm { @@ -35,434 +37,445 @@ namespace Opm class AutoDiffMatrix { public: - AutoDiffMatrix() - : type_(Z), - rows_(0), - cols_(0) - { - } - - - - AutoDiffMatrix(const int rows, const int cols) - : type_(Z), - rows_(rows), - cols_(cols) - { - } - - - - enum CreationType { ZeroMatrix, IdentityMatrix }; - - - AutoDiffMatrix(const CreationType t, const int rows) - : type_(t == ZeroMatrix ? Z : I), - rows_(rows), - cols_(rows) - { - } - - - - explicit AutoDiffMatrix(const Eigen::DiagonalMatrix& d) - : type_(D), - rows_(d.rows()), - cols_(d.cols()), - d_(d) - { - } - - - - explicit AutoDiffMatrix(const Eigen::SparseMatrix& s) - : type_(S), - rows_(s.rows()), - cols_(s.cols()), - s_(s) - { - } - - - - AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const - { - assert(rows_ == rhs.rows_); - assert(cols_ == rhs.cols_); - switch (type_) { - case Z: - return rhs; - case I: - switch (rhs.type_) { - case Z: - return *this; - case I: - return sumII(*this, rhs); - case D: - return rhs + (*this); - case S: - return rhs + (*this); - } - case D: - switch (rhs.type_) { - case Z: - return *this; - case I: - return sumDI(*this, rhs); - case D: - return sumDD(*this, rhs); - case S: - return rhs + (*this); - } - case S: - switch (rhs.type_) { - case Z: - return *this; - case I: - return sumSI(*this, rhs); - case D: - return sumSD(*this, rhs); - case S: - return sumSS(*this, rhs); - } - } - } - - AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const - { - assert(cols_ == rhs.rows_); - switch (type_) { - case Z: - return AutoDiffMatrix(rows_, rhs.cols_); - case I: - switch (rhs.type_) { - case Z: - return rhs; - case I: - return rhs; - case D: - return rhs; - case S: - return rhs; - } - case D: - switch (rhs.type_) { - case Z: - return AutoDiffMatrix(rows_, rhs.cols_); - case I: - return *this; - case D: - return prodDD(*this, rhs); - case S: - 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: - return prodSS(*this, rhs); - } - } - } - - - - - - - - AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs) - { - *this = *this + rhs; - return *this; - } - - - - - - - AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs) - { - *this = *this + rhs * -1.0; - return *this; - } - - - - - - - AutoDiffMatrix operator*(const double rhs) const - { - switch (type_) { - case Z: - return *this; - case I: - { - AutoDiffMatrix retval(*this); - retval.type_ = D; - retval.d_ = rhs * Eigen::VectorXd::Ones(rows_).asDiagonal(); - return retval; - } - case D: - { - AutoDiffMatrix retval(*this); - retval.d_ = rhs * retval.d_; - return retval; - } - case S: - { - AutoDiffMatrix retval(*this); - retval.s_ *= rhs; - return retval; - } - } - } - - - - - - - Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const - { - assert(cols_ == rhs.size()); - switch (type_) { - case Z: - return Eigen::VectorXd::Zero(rows_); - case I: - return rhs; - case D: - return d_ * rhs; - case S: - return s_ * rhs; - } - } - - - - - - - static AutoDiffMatrix sumII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == I); - assert(rhs.type_ == I); - AutoDiffMatrix retval; - 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) - { - assert(lhs.type_ == D); - assert(rhs.type_ == I); - AutoDiffMatrix retval = lhs; - for (int r = 0; r < lhs.rows_; ++r) { - retval.d_.diagonal()(r) += 1.0; - } - return retval; - } - - static AutoDiffMatrix sumDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == D); - assert(rhs.type_ == D); - AutoDiffMatrix retval = lhs; - for (int r = 0; r < lhs.rows_; ++r) { - retval.d_.diagonal()(r) += rhs.d_.diagonal()(r); - } - return retval; - } - - static AutoDiffMatrix sumSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == S); - assert(rhs.type_ == I); - AutoDiffMatrix retval; - Eigen::SparseMatrix ident = spdiag(Eigen::VectorXd::Ones(lhs.rows_)); - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - retval.s_ = lhs.s_ + ident; - return retval; - } - - static AutoDiffMatrix sumSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == S); - assert(rhs.type_ == D); - AutoDiffMatrix retval; - Eigen::SparseMatrix diag = spdiag(rhs.d_.diagonal()); - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - retval.s_ = lhs.s_ + diag; - return retval; - } - - static AutoDiffMatrix sumSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == S); - assert(rhs.type_ == S); - AutoDiffMatrix retval; - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - retval.s_ = lhs.s_ + rhs.s_; - return retval; - } - - - - - - static AutoDiffMatrix prodDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == D); - assert(rhs.type_ == D); - AutoDiffMatrix retval; - retval.type_ = D; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - 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); - assert(rhs.type_ == S); - AutoDiffMatrix retval; - // Eigen::SparseMatrix diag = spdiag(lhs.d_.diagonal()); - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - // fastSparseProduct(diag, rhs.s_, retval.s_); - fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_); - // retval.s_ = lhs.d_ * rhs.s_; - return retval; - } - - static AutoDiffMatrix prodSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == S); - assert(rhs.type_ == D); - AutoDiffMatrix retval; - // Eigen::SparseMatrix diag = spdiag(rhs.d_.diagonal()); - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - // fastSparseProduct(lhs.s_, diag, retval.s_); - fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_); - // retval.s_ = lhs.s_ * rhs.d_; - return retval; - } - - static AutoDiffMatrix prodSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) - { - assert(lhs.type_ == S); - assert(rhs.type_ == S); - AutoDiffMatrix retval; - retval.type_ = S; - retval.rows_ = lhs.rows_; - retval.cols_ = rhs.cols_; - fastSparseProduct(lhs.s_, rhs.s_, retval.s_); - return retval; - } - - - - void toSparse(Eigen::SparseMatrix& s) const - { - switch (type_) { - case Z: - s = Eigen::SparseMatrix(rows_, cols_); - return; - case I: - s = spdiag(Eigen::VectorXd::Ones(rows_)); - return; - case D: - s = spdiag(d_.diagonal()); - return; - case S: - s = s_; - return; - } - } - - - int rows() const - { - return rows_; - } - - int cols() const - { - return cols_; - } - - int nonZeros() const - { - switch (type_) { - case Z: - return 0; - case I: - return rows_; - case D: - return rows_; - case S: - return s_.nonZeros(); - } - } - - - 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_.diagonal()(row) : 0.0; - case S: - return s_.coeff(row, col); - } - } + AutoDiffMatrix() + : type_(Z), + rows_(0), + cols_(0) + { + } + + + + AutoDiffMatrix(const int rows, const int cols) + : type_(Z), + rows_(rows), + cols_(cols) + { + } + + + + enum CreationType { ZeroMatrix, IdentityMatrix }; + + + AutoDiffMatrix(const CreationType t, const int rows) + : type_(t == ZeroMatrix ? Z : I), + rows_(rows), + cols_(rows) + { + } + + + + explicit AutoDiffMatrix(const Eigen::DiagonalMatrix& d) + : type_(D), + rows_(d.rows()), + cols_(d.cols()), + d_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()) + { + } + + + + explicit AutoDiffMatrix(const Eigen::SparseMatrix& s) + : type_(S), + rows_(s.rows()), + cols_(s.cols()), + s_(s) + { + } + + + + void swap(AutoDiffMatrix& other) + { + std::swap(type_, other.type_); + std::swap(rows_, other.rows_); + std::swap(cols_, other.cols_); + d_.swap(other.d_); + s_.swap(other.s_); + } + + + + AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const + { + assert(rows_ == rhs.rows_); + assert(cols_ == rhs.cols_); + switch (type_) { + case Z: + return rhs; + case I: + switch (rhs.type_) { + case Z: + return *this; + case I: + return sumII(*this, rhs); + case D: + return rhs + (*this); + case S: + return rhs + (*this); + } + case D: + switch (rhs.type_) { + case Z: + return *this; + case I: + return sumDI(*this, rhs); + case D: + return sumDD(*this, rhs); + case S: + return rhs + (*this); + } + case S: + switch (rhs.type_) { + case Z: + return *this; + case I: + return sumSI(*this, rhs); + case D: + return sumSD(*this, rhs); + case S: + return sumSS(*this, rhs); + } + } + } + + AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const + { + assert(cols_ == rhs.rows_); + switch (type_) { + case Z: + return AutoDiffMatrix(rows_, rhs.cols_); + case I: + switch (rhs.type_) { + case Z: + return rhs; + case I: + return rhs; + case D: + return rhs; + case S: + return rhs; + } + case D: + switch (rhs.type_) { + case Z: + return AutoDiffMatrix(rows_, rhs.cols_); + case I: + return *this; + case D: + return prodDD(*this, rhs); + case S: + 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: + return prodSS(*this, rhs); + } + } + } + + + + + + + + AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs) + { + *this = *this + rhs; + return *this; + } + + + + + + + AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs) + { + *this = *this + rhs * -1.0; + return *this; + } + + + + + + + AutoDiffMatrix operator*(const double rhs) const + { + switch (type_) { + case Z: + return *this; + case I: + { + AutoDiffMatrix retval(*this); + retval.type_ = D; + retval.d_.assign(rows_, rhs); + return retval; + } + case D: + { + AutoDiffMatrix retval(*this); + for (double& elem : retval.d_) { + elem *= rhs; + } + return retval; + } + case S: + { + AutoDiffMatrix retval(*this); + retval.s_ *= rhs; + return retval; + } + } + } + + + + + + + Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const + { + assert(cols_ == rhs.size()); + switch (type_) { + case Z: + return Eigen::VectorXd::Zero(rows_); + case I: + return rhs; + case D: + return Eigen::Map(d_.data(), rows_) * rhs; + case S: + return s_ * rhs; + } + } + + + + + + + static AutoDiffMatrix sumII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == I); + assert(rhs.type_ == I); + AutoDiffMatrix retval; + retval.type_ = D; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.d_.assign(lhs.rows_, 2.0); + return retval; + } + + static AutoDiffMatrix sumDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + static_cast(rhs); // Silence release-mode warning. + assert(lhs.type_ == D); + assert(rhs.type_ == I); + AutoDiffMatrix retval = lhs; + for (int r = 0; r < lhs.rows_; ++r) { + retval.d_[r] += 1.0; + } + return retval; + } + + static AutoDiffMatrix sumDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == D); + assert(rhs.type_ == D); + AutoDiffMatrix retval = lhs; + for (int r = 0; r < lhs.rows_; ++r) { + retval.d_[r] += rhs.d_[r]; + } + return retval; + } + + static AutoDiffMatrix sumSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == S); + assert(rhs.type_ == I); + AutoDiffMatrix retval; + Eigen::SparseMatrix ident = spdiag(Eigen::VectorXd::Ones(lhs.rows_)); + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.s_ = lhs.s_ + ident; + return retval; + } + + static AutoDiffMatrix sumSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == S); + assert(rhs.type_ == D); + AutoDiffMatrix retval; + Eigen::SparseMatrix diag = spdiag(rhs.d_); + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.s_ = lhs.s_ + diag; + return retval; + } + + static AutoDiffMatrix sumSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == S); + assert(rhs.type_ == S); + AutoDiffMatrix retval; + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.s_ = lhs.s_ + rhs.s_; + return retval; + } + + + + + + static AutoDiffMatrix prodDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == D); + assert(rhs.type_ == D); + AutoDiffMatrix retval; + retval.type_ = D; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.d_.resize(lhs.rows_); + for (int r = 0; r < lhs.rows_; ++r) { + retval.d_[r] = lhs.d_[r] * rhs.d_[r]; + } + return retval; + } + + static AutoDiffMatrix prodDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == D); + assert(rhs.type_ == S); + AutoDiffMatrix retval; + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_); + return retval; + } + + static AutoDiffMatrix prodSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == S); + assert(rhs.type_ == D); + AutoDiffMatrix retval; + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_); + return retval; + } + + static AutoDiffMatrix prodSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) + { + assert(lhs.type_ == S); + assert(rhs.type_ == S); + AutoDiffMatrix retval; + retval.type_ = S; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + fastSparseProduct(lhs.s_, rhs.s_, retval.s_); + return retval; + } + + + + void toSparse(Eigen::SparseMatrix& s) const + { + switch (type_) { + case Z: + s = Eigen::SparseMatrix(rows_, cols_); + return; + case I: + s = spdiag(Eigen::VectorXd::Ones(rows_)); + return; + case D: + s = spdiag(d_); + return; + case S: + s = s_; + return; + } + } + + + int rows() const + { + return rows_; + } + + int cols() const + { + return cols_; + } + + int nonZeros() const + { + switch (type_) { + case Z: + return 0; + case I: + return rows_; + case D: + return rows_; + case S: + return s_.nonZeros(); + } + } + + + 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: - enum MatrixType { Z, I, D, S }; - MatrixType type_; - int rows_; - int cols_; - Eigen::DiagonalMatrix d_; - Eigen::SparseMatrix s_; + enum MatrixType { Z, I, D, S }; + MatrixType type_; + int rows_; + int cols_; + std::vector d_; + Eigen::SparseMatrix s_; - template - static inline - Eigen::SparseMatrix - spdiag(const V& d) - { - typedef Eigen::SparseMatrix M; - const int n = d.size(); - M mat(n, n); - mat.reserve(Eigen::ArrayXi::Ones(n, 1)); - for (M::Index i = 0; i < n; ++i) { - mat.insert(i, i) = d[i]; - } + template + static inline + Eigen::SparseMatrix + spdiag(const V& d) + { + typedef Eigen::SparseMatrix M; + const int n = d.size(); + M mat(n, n); + mat.reserve(Eigen::ArrayXi::Ones(n, 1)); + for (M::Index i = 0; i < n; ++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) { - res = lhs * rhs; + res = lhs * rhs; } inline void fastSparseProduct(const Eigen::SparseMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res) { - res = AutoDiffMatrix(lhs) * rhs; + res = AutoDiffMatrix(lhs) * rhs; } inline AutoDiffMatrix operator*(const Eigen::SparseMatrix& lhs, const AutoDiffMatrix& rhs) { - AutoDiffMatrix retval; - fastSparseProduct(lhs, rhs, retval); - return retval; + AutoDiffMatrix retval; + fastSparseProduct(lhs, rhs, retval); + return retval; } } // namespace Opm diff --git a/opm/autodiff/fastSparseProduct.hpp b/opm/autodiff/fastSparseProduct.hpp index bf94f9804..0a9923e4a 100644 --- a/opm/autodiff/fastSparseProduct.hpp +++ b/opm/autodiff/fastSparseProduct.hpp @@ -142,7 +142,8 @@ void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res) -inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix& lhs, +inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix& lhs, + const std::vector& lhs, const Eigen::SparseMatrix& rhs, Eigen::SparseMatrix& res) { @@ -153,7 +154,7 @@ inline void fastDiagSparseProduct(const Eigen::DiagonalMatrix::InnerIterator 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& lhs, - const Eigen::DiagonalMatrix& rhs, + // const Eigen::DiagonalMatrix& rhs, + const std::vector& rhs, Eigen::SparseMatrix& res) { res = lhs; @@ -172,7 +174,7 @@ inline void fastSparseDiagProduct(const Eigen::SparseMatrix& lhs, for (int col = 0; col < n; ++col) { typedef Eigen::SparseMatrix::InnerIterator It; for (It it(res, col); it; ++it) { - it.valueRef() *= rhs.diagonal()(col); + it.valueRef() *= rhs[col]; // rhs.diagonal()(col); } } }