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

Merge pull request #482 from babrodtk/autodiffmatrix_documentation

Browse Source 
Autodiffmatrix beautification
This commit is contained in:
Atgeirr Flø Rasmussen 2015-10-02 12:48:16 +02:00
commit 85d7f6c904
3 changed files with 158 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/AutoDiffBlock.hpp
View File

@ -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));
}

158
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -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<double> 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<double, Eigen::Dynamic>& 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<double>& 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<void>(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<void>(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<class Scalar, int Options, class Index>
void toSparse(Eigen::SparseMatrix<Scalar, Options, Index>& 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 <class V>
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<double>& lhs, const AutoDiffMatrix& rhs, AutoDiffMatrix& res)
{
res = AutoDiffMatrix(lhs) * rhs;
}
/**
* Multiplies an Eigen sparse matrix with an AutoDiffMatrix.
*/
inline AutoDiffMatrix operator*(const Eigen::SparseMatrix<double>& lhs, const AutoDiffMatrix& rhs)
{
AutoDiffMatrix retval;

32
tests/test_autodiffmatrix.cpp
View File

@ -74,9 +74,9 @@ operator ==(const Eigen::SparseMatrix<double>& 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<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3).sparseView());
Eigen::Array<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3).sparseView());
Eigen::Array<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3).sparseView());
Eigen::Array<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3).sparseView());
Eigen::Array<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3).sparseView());
Eigen::Array<double, Eigen::Dynamic, 1> 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<double, Eigen::Dynamic, 1> d1(3);
d1 << 0.2, 1.2, 13.4;