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..c1469cc59 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), @@ -64,20 +74,20 @@ 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); } + /** + * Creates a diagonal matrix from an Eigen diagonal matrix + */ explicit AutoDiffMatrix(const Eigen::DiagonalMatrix& d) : type_(Diagonal), rows_(d.rows()), @@ -90,6 +100,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 +150,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 +207,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 +284,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 +325,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 +366,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 +395,7 @@ namespace Opm - + // Add identity to identity static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Identity); @@ -364,6 +408,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 +421,7 @@ namespace Opm return retval; } + // Add diagonal to diagonal static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -387,6 +433,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 +444,7 @@ namespace Opm return retval; } + // Add sparse to diagonal static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -406,6 +454,7 @@ namespace Opm return retval; } + // Add sparse to sparse static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -418,7 +467,7 @@ namespace Opm - + // Multiply diagonal with diagonal static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -430,6 +479,7 @@ namespace Opm return retval; } + // Multiply diagonal with sparse static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Diagonal); @@ -442,6 +492,7 @@ namespace Opm return retval; } + // Multiply sparse with diagonal static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -454,6 +505,7 @@ namespace Opm return retval; } + // Multiply sparse with sparse static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) { assert(lhs.type_ == Sparse); @@ -467,6 +519,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 +546,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 +590,11 @@ namespace Opm } + + + /** + * Returns element (row, col) in the matrix + */ double coeff(const int row, const int col) const { switch (type_) { @@ -530,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) { /** @@ -545,14 +635,40 @@ namespace Opm return sparse_; } + 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) + + + + /** + * Constructor which sets all members + */ + AutoDiffMatrix(AudoDiffMatrixType type, int rows, int cols, + DiagRep diag=DiagRep(), SparseRep sparse=SparseRep()) + : 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 +690,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; 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;