use correct types of SparseMatrices.

This commit is contained in:
Robert K 2014-12-02 10:12:27 +01:00
parent a266e98bac
commit 2ac6a211b2
2 changed files with 20 additions and 16 deletions

View File

@ -189,11 +189,11 @@ namespace {
template <typename Scalar, class IntVec>
Eigen::SparseMatrix<Scalar>
typename AutoDiffBlock<Scalar>::M
constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
{
const int subset_size = indices.size();
Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
typename AutoDiffBlock<Scalar>::M mat(full_size, subset_size);
mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
for (int i = 0; i < subset_size; ++i) {
mat.insert(indices[i], i) = 1;
@ -204,18 +204,6 @@ namespace {
} // anon namespace
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices)
{
Eigen::SparseMatrix<Scalar> sub
= constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose();
return sub * x;
}
/// Returns x(indices).
template <typename Scalar, class IntVec>
@ -223,9 +211,25 @@ Eigen::Array<Scalar, Eigen::Dynamic, 1>
subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices)
{
return (constructSupersetSparseMatrix<Scalar>(x.size(), indices).transpose() * x.matrix()).array();
typedef typename Eigen::Array<Scalar, Eigen::Dynamic, 1>::Index Index;
const Index size = indices.size();
Eigen::Array<Scalar, Eigen::Dynamic, 1> ret( size );
for( Index i=0; i<size; ++i )
ret[ i ] = x[ indices[ i ] ];
return std::move(ret);
}
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices)
{
const typename AutoDiffBlock<Scalar>::M sub
= constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose();
return sub * x;
}
/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.

View File

@ -260,7 +260,7 @@ namespace Opm
#endif
M id(Jn[n].rows(), Jn[n].cols());
id.setIdentity();
const M Di = solver.solve(id);
const Eigen::SparseMatrix<typename M::Scalar, Eigen::ColMajor> Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());