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 #481 from babrodtk/faster-form-interleaved

Browse Source 
Faster form interleaved
This commit is contained in:
Bård Skaflestad
2015-09-24 09:19:54 +02:00
parent 2c84414641 85cfd9768e
commit a333755a2b
2 changed files with 166 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files

202
opm/autodiff/AutoDiffMatrix.hpp
View File

@@ -38,11 +38,12 @@ namespace Opm
class AutoDiffMatrix class AutoDiffMatrix
{ {
public: public:
typedef std::vector<double> Diag; typedef std::vector<double> DiagRep;
typedef Eigen::SparseMatrix<double> Sparse; typedef Eigen::SparseMatrix<double> SparseRep;
AutoDiffMatrix() AutoDiffMatrix()
: type_(Z), : type_(Zero),
rows_(0), rows_(0),
cols_(0), cols_(0),
diag_(), diag_(),
@@ -53,7 +54,7 @@ namespace Opm
AutoDiffMatrix(const int num_rows, const int num_cols) AutoDiffMatrix(const int num_rows, const int num_cols)
: type_(Z), : type_(Zero),
rows_(num_rows), rows_(num_rows),
cols_(num_cols), cols_(num_cols),
diag_(), diag_(),
@@ -67,7 +68,7 @@ namespace Opm
AutoDiffMatrix(const CreationType t, const int num_rows) AutoDiffMatrix(const CreationType t, const int num_rows)
: type_(t == ZeroMatrix ? Z : I), : type_(t == ZeroMatrix ? Zero : Identity),
rows_(num_rows), rows_(num_rows),
cols_(num_rows), cols_(num_rows),
diag_(), diag_(),
@@ -78,7 +79,7 @@ namespace Opm
explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d) explicit AutoDiffMatrix(const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& d)
: type_(D), : type_(Diagonal),
rows_(d.rows()), rows_(d.rows()),
cols_(d.cols()), cols_(d.cols()),
diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()), diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()),
@@ -90,7 +91,7 @@ namespace Opm
explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s) explicit AutoDiffMatrix(const Eigen::SparseMatrix<double>& s)
: type_(S), : type_(Sparse),
rows_(s.rows()), rows_(s.rows()),
cols_(s.cols()), cols_(s.cols()),
diag_(), diag_(),
@@ -106,7 +107,7 @@ namespace Opm
AutoDiffMatrix(AutoDiffMatrix&& other) AutoDiffMatrix(AutoDiffMatrix&& other)
: type_(Z), : type_(Zero),
rows_(0), rows_(0),
cols_(0), cols_(0),
diag_(), diag_(),
@@ -141,43 +142,43 @@ namespace Opm
assert(rows_ == rhs.rows_); assert(rows_ == rhs.rows_);
assert(cols_ == rhs.cols_); assert(cols_ == rhs.cols_);
switch (type_) { switch (type_) {
case Z: case Zero:
return rhs; return rhs;
case I: case Identity:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Zero:
return *this; return *this;
case I: case Identity:
return addII(*this, rhs); return addII(*this, rhs);
case D: case Diagonal:
return addDI(rhs, *this); return addDI(rhs, *this);
case S: case Sparse:
return addSI(rhs, *this); return addSI(rhs, *this);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
} }
case D: case Diagonal:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Zero:
return *this; return *this;
case I: case Identity:
return addDI(*this, rhs); return addDI(*this, rhs);
case D: case Diagonal:
return addDD(*this, rhs); return addDD(*this, rhs);
case S: case Sparse:
return addSD(rhs, *this); return addSD(rhs, *this);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
} }
case S: case Sparse:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Zero:
return *this; return *this;
case I: case Identity:
return addSI(*this, rhs); return addSI(*this, rhs);
case D: case Diagonal:
return addSD(*this, rhs); return addSD(*this, rhs);
case S: case Sparse:
return addSS(*this, rhs); return addSS(*this, rhs);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
@@ -191,32 +192,32 @@ namespace Opm
{ {
assert(cols_ == rhs.rows_); assert(cols_ == rhs.rows_);
switch (type_) { switch (type_) {
case Z: case Zero:
return AutoDiffMatrix(rows_, rhs.cols_); return AutoDiffMatrix(rows_, rhs.cols_);
case I: case Identity:
return rhs; return rhs;
case D: case Diagonal:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Zero:
return AutoDiffMatrix(rows_, rhs.cols_); return AutoDiffMatrix(rows_, rhs.cols_);
case I: case Identity:
return *this; return *this;
case D: case Diagonal:
return mulDD(*this, rhs); return mulDD(*this, rhs);
case S: case Sparse:
return mulDS(*this, rhs); return mulDS(*this, rhs);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
} }
case S: case Sparse:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Zero:
return AutoDiffMatrix(rows_, rhs.cols_); return AutoDiffMatrix(rows_, rhs.cols_);
case I: case Identity:
return *this; return *this;
case D: case Diagonal:
return mulSD(*this, rhs); return mulSD(*this, rhs);
case S: case Sparse:
return mulSS(*this, rhs); return mulSS(*this, rhs);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
@@ -257,16 +258,16 @@ namespace Opm
AutoDiffMatrix operator*(const double rhs) const AutoDiffMatrix operator*(const double rhs) const
{ {
switch (type_) { switch (type_) {
case Z: case Zero:
return *this; return *this;
case I: case Identity:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = Diagonal;
retval.diag_.assign(rows_, rhs); retval.diag_.assign(rows_, rhs);
return retval; return retval;
} }
case D: case Diagonal:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
for (double& elem : retval.diag_) { for (double& elem : retval.diag_) {
@@ -274,7 +275,7 @@ namespace Opm
} }
return retval; return retval;
} }
case S: case Sparse:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.sparse_ *= rhs; retval.sparse_ *= rhs;
@@ -293,16 +294,16 @@ namespace Opm
AutoDiffMatrix operator/(const double rhs) const AutoDiffMatrix operator/(const double rhs) const
{ {
switch (type_) { switch (type_) {
case Z: case Zero:
return *this; return *this;
case I: case Identity:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = Diagonal;
retval.diag_.assign(rows_, 1.0/rhs); retval.diag_.assign(rows_, 1.0/rhs);
return retval; return retval;
} }
case D: case Diagonal:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
for (double& elem : retval.diag_) { for (double& elem : retval.diag_) {
@@ -310,7 +311,7 @@ namespace Opm
} }
return retval; return retval;
} }
case S: case Sparse:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.sparse_ /= rhs; retval.sparse_ /= rhs;
@@ -330,16 +331,16 @@ namespace Opm
{ {
assert(cols_ == rhs.size()); assert(cols_ == rhs.size());
switch (type_) { switch (type_) {
case Z: case Zero:
return Eigen::VectorXd::Zero(rows_); return Eigen::VectorXd::Zero(rows_);
case I: case Identity:
return rhs; return rhs;
case D: case Diagonal:
{ {
const Eigen::VectorXd d = Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_); const Eigen::VectorXd d = Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_);
return d.cwiseProduct(rhs); return d.cwiseProduct(rhs);
} }
case S: case Sparse:
return sparse_ * rhs; return sparse_ * rhs;
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
@@ -353,10 +354,10 @@ namespace Opm
static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addII(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == I); assert(lhs.type_ == Identity);
assert(rhs.type_ == I); assert(rhs.type_ == Identity);
AutoDiffMatrix retval; AutoDiffMatrix retval;
retval.type_ = D; retval.type_ = Diagonal;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.diag_.assign(lhs.rows_, 2.0); retval.diag_.assign(lhs.rows_, 2.0);
@@ -366,8 +367,8 @@ namespace Opm
static AutoDiffMatrix addDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addDI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
static_cast<void>(rhs); // Silence release-mode warning. static_cast<void>(rhs); // Silence release-mode warning.
assert(lhs.type_ == D); assert(lhs.type_ == Diagonal);
assert(rhs.type_ == I); assert(rhs.type_ == Identity);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.diag_[r] += 1.0; retval.diag_[r] += 1.0;
@@ -377,8 +378,8 @@ namespace Opm
static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == D); assert(lhs.type_ == Diagonal);
assert(rhs.type_ == D); assert(rhs.type_ == Diagonal);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.diag_[r] += rhs.diag_[r]; retval.diag_[r] += rhs.diag_[r];
@@ -389,8 +390,8 @@ namespace Opm
static AutoDiffMatrix addSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addSI(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
static_cast<void>(rhs); // Silence release-mode warning. static_cast<void>(rhs); // Silence release-mode warning.
assert(lhs.type_ == S); assert(lhs.type_ == Sparse);
assert(rhs.type_ == I); assert(rhs.type_ == Identity);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
retval.sparse_ += spdiag(Eigen::VectorXd::Ones(lhs.rows_));; retval.sparse_ += spdiag(Eigen::VectorXd::Ones(lhs.rows_));;
return retval; return retval;
@@ -398,8 +399,8 @@ namespace Opm
static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == S); assert(lhs.type_ == Sparse);
assert(rhs.type_ == D); assert(rhs.type_ == Diagonal);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
retval.sparse_ += spdiag(rhs.diag_); retval.sparse_ += spdiag(rhs.diag_);
return retval; return retval;
@@ -407,8 +408,8 @@ namespace Opm
static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix addSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == S); assert(lhs.type_ == Sparse);
assert(rhs.type_ == S); assert(rhs.type_ == Sparse);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
retval.sparse_ += rhs.sparse_; retval.sparse_ += rhs.sparse_;
return retval; return retval;
@@ -420,8 +421,8 @@ namespace Opm
static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix mulDD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == D); assert(lhs.type_ == Diagonal);
assert(rhs.type_ == D); assert(rhs.type_ == Diagonal);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.diag_[r] *= rhs.diag_[r]; retval.diag_[r] *= rhs.diag_[r];
@@ -431,10 +432,10 @@ namespace Opm
static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix mulDS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == D); assert(lhs.type_ == Diagonal);
assert(rhs.type_ == S); assert(rhs.type_ == Sparse);
AutoDiffMatrix retval; AutoDiffMatrix retval;
retval.type_ = S; retval.type_ = Sparse;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastDiagSparseProduct(lhs.diag_, rhs.sparse_, retval.sparse_); fastDiagSparseProduct(lhs.diag_, rhs.sparse_, retval.sparse_);
@@ -443,10 +444,10 @@ namespace Opm
static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix mulSD(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == S); assert(lhs.type_ == Sparse);
assert(rhs.type_ == D); assert(rhs.type_ == Diagonal);
AutoDiffMatrix retval; AutoDiffMatrix retval;
retval.type_ = S; retval.type_ = Sparse;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseDiagProduct(lhs.sparse_, rhs.diag_, retval.sparse_); fastSparseDiagProduct(lhs.sparse_, rhs.diag_, retval.sparse_);
@@ -455,10 +456,10 @@ namespace Opm
static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs) static AutoDiffMatrix mulSS(const AutoDiffMatrix& lhs, const AutoDiffMatrix& rhs)
{ {
assert(lhs.type_ == S); assert(lhs.type_ == Sparse);
assert(rhs.type_ == S); assert(rhs.type_ == Sparse);
AutoDiffMatrix retval; AutoDiffMatrix retval;
retval.type_ = S; retval.type_ = Sparse;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseProduct(lhs.sparse_, rhs.sparse_, retval.sparse_); fastSparseProduct(lhs.sparse_, rhs.sparse_, retval.sparse_);
@@ -470,16 +471,16 @@ namespace Opm
void toSparse(Eigen::SparseMatrix<Scalar, Options, Index>& s) const void toSparse(Eigen::SparseMatrix<Scalar, Options, Index>& s) const
{ {
switch (type_) { switch (type_) {
case Z: case Zero:
s = Eigen::SparseMatrix<Scalar, Options, Index>(rows_, cols_); s = Eigen::SparseMatrix<Scalar, Options, Index>(rows_, cols_);
return; return;
case I: case Identity:
s = spdiag(Eigen::VectorXd::Ones(rows_)); s = spdiag(Eigen::VectorXd::Ones(rows_));
return; return;
case D: case Diagonal:
s = spdiag(diag_); s = spdiag(diag_);
return; return;
case S: case Sparse:
s = sparse_; s = sparse_;
return; return;
} }
@@ -499,13 +500,13 @@ namespace Opm
int nonZeros() const int nonZeros() const
{ {
switch (type_) { switch (type_) {
case Z: case Zero:
return 0; return 0;
case I: case Identity:
return rows_; return rows_;
case D: case Diagonal:
return rows_; return rows_;
case S: case Sparse:
return sparse_.nonZeros(); return sparse_.nonZeros();
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
@@ -516,36 +517,51 @@ namespace Opm
double coeff(const int row, const int col) const double coeff(const int row, const int col) const
{ {
switch (type_) { switch (type_) {
case Z: case Zero:
return 0.0; return 0.0;
case I: case Identity:
return (row == col) ? 1.0 : 0.0; return (row == col) ? 1.0 : 0.0;
case D: case Diagonal:
return (row == col) ? diag_[row] : 0.0; return (row == col) ? diag_[row] : 0.0;
case S: case Sparse:
return sparse_.coeff(row, col); return sparse_.coeff(row, col);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
} }
const SparseRep& getSparse() const {
if (type_ != Sparse) {
/**
* 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.
*/
SparseRep& mutable_sparse = const_cast<SparseRep&>(sparse_);
toSparse(mutable_sparse);
}
return sparse_;
}
private: private:
enum MatrixType { Z, I, D, S }; enum AudoDiffMatrixType { Zero, Identity, Diagonal, Sparse };
MatrixType type_; AudoDiffMatrixType type_;
int rows_; int rows_;
int cols_; int cols_;
Diag diag_; DiagRep diag_;
Sparse sparse_; SparseRep sparse_;
template <class V> template <class V>
static inline static inline
Sparse SparseRep
spdiag(const V& d) spdiag(const V& d)
{ {
const int n = d.size(); const int n = d.size();
Sparse mat(n, n); SparseRep mat(n, n);
mat.reserve(Eigen::ArrayXi::Ones(n, 1)); mat.reserve(Eigen::ArrayXi::Ones(n, 1));
for (Sparse::Index i = 0; i < n; ++i) { for (SparseRep::Index i = 0; i < n; ++i) {
if (d[i] != 0.0) { if (d[i] != 0.0) {
mat.insert(i, i) = d[i]; mat.insert(i, i) = d[i];
} }

79
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@@ -198,52 +198,87 @@ namespace Opm
namespace detail {
/**
* Simple binary operator that always returns 0.1
* It is used to get the sparsity pattern for our
* interleaved system, and is marginally faster than using
* operator+=.
*/
template<typename Scalar> struct PointOneOp {
EIGEN_EMPTY_STRUCT_CTOR(PointOneOp)
Scalar operator()(const Scalar& a, const Scalar& b) const { return 0.1; }
};
}
void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs, void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs,
Mat& istlA) const Mat& istlA) const
{ {
const int np = eqs.size(); const int np = eqs.size();
// Find sparsity structure as union of basic block sparsity structures, // Find sparsity structure as union of basic block sparsity structures,
// corresponding to the jacobians with respect to pressure. // corresponding to the jacobians with respect to pressure.
// Use addition to get to the union structure. // Use our custom PointOneOp to get to the union structure.
typedef Eigen::SparseMatrix<double> Sp; // Note that we only iterate over the pressure derivatives on purpose.
Sp structure; Eigen::SparseMatrix<double, Eigen::ColMajor> col_major = eqs[0].derivative()[0].getSparse();
eqs[0].derivative()[0].toSparse(structure); detail::PointOneOp<double> point_one;
{ for (int phase = 1; phase < np; ++phase) {
Sp s0; const AutoDiffMatrix::SparseRep& mat = eqs[phase].derivative()[0].getSparse();
for (int phase = 1; phase < np; ++phase) { col_major = col_major.binaryExpr(mat, point_one);
eqs[phase].derivative()[0].toSparse(s0);
structure += s0;
}
} }
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure; // Automatically convert the column major structure to a row-major structure
Eigen::SparseMatrix<double, Eigen::RowMajor> row_major = col_major;
const int size = row_major.rows();
assert(size == row_major.cols());
// Create ISTL matrix with interleaved rows and columns (block structured). // Create ISTL matrix with interleaved rows and columns (block structured).
assert(np == 3); assert(np == 3);
istlA.setSize(s.rows(), s.cols(), s.nonZeros()); istlA.setSize(row_major.rows(), row_major.cols(), row_major.nonZeros());
istlA.setBuildMode(Mat::row_wise); istlA.setBuildMode(Mat::row_wise);
const int* ia = s.outerIndexPtr(); const int* ia = row_major.outerIndexPtr();
const int* ja = s.innerIndexPtr(); const int* ja = row_major.innerIndexPtr();
for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) { for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) {
int ri = row.index(); const int ri = row.index();
for (int i = ia[ri]; i < ia[ri + 1]; ++i) { for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
row.insert(ja[i]); row.insert(ja[i]);
} }
} }
const int size = s.rows(); // Set all blocks to zero.
for (int row = 0; row < size; ++row) { for (int row = 0; row < size; ++row) {
for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) { for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
const int col = ja[col_ix]; const int col = ja[col_ix];
MatrixBlockType block; istlA[row][col] = 0.0;
for (int p1 = 0; p1 < np; ++p1) { }
for (int p2 = 0; p2 < np; ++p2) { }
block[p1][p2] = eqs[p1].derivative()[p2].coeff(row, col);
/**
* Go through all jacobians, and insert in correct spot
*
* The straight forward way to do this would be to run through each
* element in the output matrix, and set all block entries by gathering
* from all "input matrices" (derivatives).
*
* A faster alternative is to instead run through each "input matrix" and
* insert its elements in the correct spot in the output matrix.
*
*/
for (int col = 0; col < size; ++col) {
for (int p1 = 0; p1 < np; ++p1) {
for (int p2 = 0; p2 < np; ++p2) {
// Note that that since these are CSC and not CSR matrices,
// ja contains row numbers instead of column numbers.
const AutoDiffMatrix::SparseRep& s = eqs[p1].derivative()[p2].getSparse();
const int* ia = s.outerIndexPtr();
const int* ja = s.innerIndexPtr();
const double* sa = s.valuePtr();
for (int elem_ix = ia[col]; elem_ix < ia[col + 1]; ++elem_ix) {
const int row = ja[elem_ix];
istlA[row][col][p1][p2] = sa[elem_ix];
} }
} }
istlA[row][col] = block;
} }
} }
} }