From 4073abeafbadcc653575869f2443fc838681ba41 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Mon, 21 Sep 2015 09:19:56 +0200 Subject: [PATCH 1/3] Added documentation for AutoDiffMatrix --- opm/autodiff/AutoDiffBlock.hpp | 2 +- opm/autodiff/AutoDiffMatrix.hpp | 155 ++++++++++++++++++++++++++++++-- 2 files changed, 147 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index c512a0507..ede98336d 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -158,7 +158,7 @@ namespace Opm } // ... then set the one corrresponding to this variable to identity. assert(blocksizes[index] == num_elem); - jac[index] = M(M::IdentityMatrix, val.size()); + jac[index] = M::createIdentity(val.size()); return AutoDiffBlock(std::move(val), std::move(jac)); } diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 8b1a5b88d..4d6981633 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -35,6 +35,11 @@ namespace Opm { + /** + * AutoDiffMatrix is a wrapper class that optimizes matrix operations. + * Internally, an AutoDiffMatrix can be either Zero, Identity, Diagonal, + * or Sparse, and we utilize this to perform faster matrix operations. + */ class AutoDiffMatrix { public: @@ -42,6 +47,9 @@ namespace Opm typedef Eigen::SparseMatrix SparseRep; + /** + * Creates an empty zero matrix + */ AutoDiffMatrix() : type_(Zero), rows_(0), @@ -52,7 +60,9 @@ namespace Opm } - + /** + * Creates a zero matrix with num_rows x num_cols entries + */ AutoDiffMatrix(const int num_rows, const int num_cols) : type_(Zero), rows_(num_rows), @@ -78,6 +88,9 @@ namespace Opm + /** + * Creates a diagonal matrix from an Eigen diagonal matrix + */ explicit AutoDiffMatrix(const Eigen::DiagonalMatrix& d) : type_(Diagonal), rows_(d.rows()), @@ -90,6 +103,9 @@ namespace Opm + /** + * Creates a sparse matrix from an Eigen sparse matrix + */ explicit AutoDiffMatrix(const Eigen::SparseMatrix& s) : type_(Sparse), rows_(s.rows()), @@ -137,6 +153,12 @@ namespace Opm + /** + * Adds two AutoDiffMatrices. Internally, this function optimizes + * the addition operation based on the structure of the matrix, e.g., + * adding two zero matrices will obviously yield a zero matrix, and + * so on. + */ AutoDiffMatrix operator+(const AutoDiffMatrix& rhs) const { assert(rows_ == rhs.rows_); @@ -188,6 +210,16 @@ namespace Opm } } + + + + + + /** + * Multiplies two AutoDiffMatrices. Internally, this function optimizes + * the multiplication operation based on the structure of the matrix, e.g., + * multiplying M with a zero matrix will obviously yield a zero matrix. + */ AutoDiffMatrix operator*(const AutoDiffMatrix& rhs) const { assert(cols_ == rhs.rows_); @@ -255,6 +287,11 @@ namespace Opm + /** + * Multiplies an AutoDiffMatrix with a scalar. Optimizes internally + * by exploiting that e.g., an identity matrix multiplied by a scalar x + * yields a diagonal matrix with x the diagonal. + */ AutoDiffMatrix operator*(const double rhs) const { switch (type_) { @@ -291,6 +328,11 @@ namespace Opm + /** + * Divides an AutoDiffMatrix by a scalar. Optimizes internally + * by exploiting that e.g., an identity matrix divided by a scalar x + * yields a diagonal matrix with 1/x on the diagonal. + */ AutoDiffMatrix operator/(const double rhs) const { switch (type_) { @@ -327,6 +369,11 @@ namespace Opm + /** + * Multiplies an AutoDiffMatrix with a vector. Optimizes internally + * by exploiting that e.g., an identity matrix multiplied by a vector + * yields the vector itself. + */ Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const { assert(cols_ == rhs.size()); @@ -351,7 +398,7 @@ namespace Opm - + // Add identity to identity static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Identity); @@ -364,6 +411,7 @@ namespace Opm return retval; } + // Add diagonal to identity static AutoDiffMatrix addDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { static_cast(rhs); // Silence release-mode warning. @@ -376,6 +424,7 @@ namespace Opm return retval; } + // Add diagonal to diagonal static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -387,6 +436,7 @@ namespace Opm return retval; } + // Add sparse to identity static AutoDiffMatrix addSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { static_cast(rhs); // Silence release-mode warning. @@ -397,6 +447,7 @@ namespace Opm return retval; } + // Add sparse to diagonal static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -406,6 +457,7 @@ namespace Opm return retval; } + // Add sparse to sparse static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -418,7 +470,7 @@ namespace Opm - + // Multiply diagonal with diagonal static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -430,6 +482,7 @@ namespace Opm return retval; } + // Multiply diagonal with sparse static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -442,6 +495,7 @@ namespace Opm return retval; } + // Multiply sparse with diagonal static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -454,6 +508,7 @@ namespace Opm return retval; } + // Multiply sparse with sparse static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -467,6 +522,13 @@ namespace Opm } + + + /** + * Converts the AutoDiffMatrix to an Eigen SparseMatrix.This might be + * an expensive operation to perform for e.g., an identity matrix or a + * diagonal matrix. + */ template void toSparse(Eigen::SparseMatrix& s) const { @@ -487,16 +549,33 @@ namespace Opm } + + /** + * Returns number of rows in the matrix + */ int rows() const { return rows_; } + + + /** + * Returns number of columns in the matrix + */ int cols() const { return cols_; } + + + /** + * Returns number of non-zero elements in the matrix. Optimizes internally + * by exploiting that e.g., an n*n identity matrix has n non-zeros. + * Note that an n*n diagonal matrix is defined to have n non-zeros, even though + * several diagonal elements might be 0.0. + */ int nonZeros() const { switch (type_) { @@ -514,6 +593,11 @@ namespace Opm } + + + /** + * Returns element (row, col) in the matrix + */ double coeff(const int row, const int col) const { switch (type_) { @@ -547,12 +631,57 @@ namespace Opm private: enum AudoDiffMatrixType { Zero, Identity, Diagonal, Sparse }; - AudoDiffMatrixType type_; - int rows_; - int cols_; - DiagRep diag_; - SparseRep sparse_; + AudoDiffMatrixType type_; //< Type of matrix + int rows_; //< Number of rows + int cols_; //< Number of columns + DiagRep diag_; //< Diagonal representation (only if type==Diagonal) + SparseRep sparse_; //< Sparse representation (only if type==Sparse) + + /** + * Returns the sparse representation of this matrix. Note that this might + * be an expensive operation to perform if the internal structure is not + * sparse. + */ + const Sparse& getSparse() const { + if (type_ != S) { + /** + * If we are not a sparse matrix, our internal variable sparse_ + * is undefined, and hence changing it so that it happens to be + * a sparse representation of our true data does not change our + * true data, and hence justifies that we do not really violate + * the const qualifier. + */ + Sparse& mutable_sparse = const_cast(sparse_); + toSparse(mutable_sparse); + } + return sparse_; + } + + + + /** + * Constructor which sets all members + */ + AutoDiffMatrix(MatrixType type, int rows, int cols, + Diag diag=Diag(), Sparse sparse=Sparse()) + : type_(type), + rows_(rows), + cols_(cols), + diag_(diag), + sparse_(sparse) + { + } + + + + + + /** + * Creates a sparse diagonal matrix from d. + * Typical use is to convert a standard vector to an + * Eigen sparse matrix. + */ template static inline SparseRep @@ -574,19 +703,27 @@ namespace Opm - + /** + * Utility function to lessen code changes required elsewhere. + */ inline void fastSparseProduct(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res) { res = lhs * rhs; } + /** + * Utility function to lessen code changes required elsewhere. + */ inline void fastSparseProduct(const Eigen::SparseMatrix& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res) { res = AutoDiffMatrix(lhs) * rhs; } + /** + * Multiplies an Eigen sparse matrix with an AutoDiffMatrix. + */ inline AutoDiffMatrix operator*(const Eigen::SparseMatrix& lhs, const AutoDiffMatrix& rhs) { AutoDiffMatrix retval; From 7f9175e0463c8c68d0180cb262f514c41156e055 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Mon, 21 Sep 2015 09:20:28 +0200 Subject: [PATCH 2/3] Created named constructor for identity matrices --- opm/autodiff/AutoDiffMatrix.hpp | 13 +++++-------- tests/test_autodiffmatrix.cpp | 32 ++++++++++++++++---------------- 2 files changed, 21 insertions(+), 24 deletions(-) diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 4d6981633..1ec4d9b22 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -74,16 +74,13 @@ namespace Opm - enum CreationType { ZeroMatrix, IdentityMatrix }; - - AutoDiffMatrix(const CreationType t, const int num_rows) - : type_(t == ZeroMatrix ? Zero : Identity), - rows_(num_rows), - cols_(num_rows), - diag_(), - sparse_() + /** + * Creates an identity matrix with num_rows_cols x num_rows_cols entries + */ + static AutoDiffMatrix createIdentity(const int num_rows_cols) { + return AutoDiffMatrix(Identity, num_rows_cols, num_rows_cols); } diff --git a/tests/test_autodiffmatrix.cpp b/tests/test_autodiffmatrix.cpp index b2c4b9231..25431f911 100644 --- a/tests/test_autodiffmatrix.cpp +++ b/tests/test_autodiffmatrix.cpp @@ -74,9 +74,9 @@ operator ==(const Eigen::SparseMatrix& A, BOOST_AUTO_TEST_CASE(Initialization) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Eigen::Array d1(3); d1 << 0.2, 1.2, 13.4; @@ -93,9 +93,9 @@ BOOST_AUTO_TEST_CASE(Initialization) BOOST_AUTO_TEST_CASE(EigenConversion) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Eigen::Array d1(3); d1 << 0.2, 1.2, 13.4; @@ -128,10 +128,10 @@ BOOST_AUTO_TEST_CASE(EigenConversion) BOOST_AUTO_TEST_CASE(AdditionOps) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); Sp zs(3,3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Sp is(Eigen::Matrix::Identity(3,3).sparseView()); Eigen::Array d1(3); @@ -176,10 +176,10 @@ BOOST_AUTO_TEST_CASE(AdditionOps) BOOST_AUTO_TEST_CASE(MultOps) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); Sp zs(3,3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Sp is(Eigen::Matrix::Identity(3,3).sparseView()); Eigen::Array d1(3); @@ -267,10 +267,10 @@ BOOST_AUTO_TEST_CASE(MultOps) BOOST_AUTO_TEST_CASE(MultOpsDouble) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); Sp zs(3,3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Sp is(Eigen::Matrix::Identity(3,3).sparseView()); Eigen::Array d1(3); @@ -307,10 +307,10 @@ BOOST_AUTO_TEST_CASE(MultOpsDouble) BOOST_AUTO_TEST_CASE(DivOpsDouble) { // Setup. - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); Sp zs(3,3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Sp is(Eigen::Matrix::Identity(3,3).sparseView()); Eigen::Array d1(3); @@ -345,10 +345,10 @@ BOOST_AUTO_TEST_CASE(DivOpsDouble) BOOST_AUTO_TEST_CASE(MultVectorXd) { - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); Sp zs(3,3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Sp is(Eigen::Matrix::Identity(3,3).sparseView()); Eigen::Array d1(3); @@ -398,9 +398,9 @@ BOOST_AUTO_TEST_CASE(Coeff) BOOST_AUTO_TEST_CASE(nonZeros) { - Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Mat z = Mat(3, 3); - Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Mat i = Mat::createIdentity(3); Eigen::Array d1(3); d1 << 0.2, 1.2, 13.4; From 654e4a81bf84483d6bfec093e8c0d030a782b808 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Fri, 2 Oct 2015 10:33:55 +0200 Subject: [PATCH 3/3] Fixed errors from rebasing --- opm/autodiff/AutoDiffMatrix.hpp | 34 ++++++++++++--------------------- 1 file changed, 12 insertions(+), 22 deletions(-) diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 1ec4d9b22..c1469cc59 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -611,6 +611,15 @@ namespace Opm } } + + + + + /** + * Returns the sparse representation of this matrix. Note that this might + * be an expensive operation to perform if the internal structure is not + * sparse. + */ const SparseRep& getSparse() const { if (type_ != Sparse) { /** @@ -626,6 +635,7 @@ namespace Opm return sparse_; } + private: enum AudoDiffMatrixType { Zero, Identity, Diagonal, Sparse }; @@ -635,33 +645,13 @@ namespace Opm DiagRep diag_; //< Diagonal representation (only if type==Diagonal) SparseRep sparse_; //< Sparse representation (only if type==Sparse) - /** - * Returns the sparse representation of this matrix. Note that this might - * be an expensive operation to perform if the internal structure is not - * sparse. - */ - const Sparse& getSparse() const { - if (type_ != S) { - /** - * If we are not a sparse matrix, our internal variable sparse_ - * is undefined, and hence changing it so that it happens to be - * a sparse representation of our true data does not change our - * true data, and hence justifies that we do not really violate - * the const qualifier. - */ - Sparse& mutable_sparse = const_cast(sparse_); - toSparse(mutable_sparse); - } - return sparse_; - } - /** * Constructor which sets all members */ - AutoDiffMatrix(MatrixType type, int rows, int cols, - Diag diag=Diag(), Sparse sparse=Sparse()) + AutoDiffMatrix(AudoDiffMatrixType type, int rows, int cols, + DiagRep diag=DiagRep(), SparseRep sparse=SparseRep()) : type_(type), rows_(rows), cols_(cols),