diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 330e22e6c..4c62c99c1 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -330,10 +330,11 @@ namespace Opm { assert(lhs.type_ == D); assert(rhs.type_ == D); - AutoDiffMatrix retval = lhs; - for (int r = 0; r < lhs.rows_; ++r) { - retval.d_.diagonal().array() *= rhs.d_.diagonal().array(); - } + AutoDiffMatrix retval; + retval.type_ = D; + retval.rows_ = lhs.rows_; + retval.cols_ = rhs.cols_; + retval.d_ = (lhs.d_.diagonal().array() * rhs.d_.diagonal().array()).matrix().asDiagonal(); return retval; } diff --git a/tests/test_autodiffmatrix.cpp b/tests/test_autodiffmatrix.cpp index 03a7c09a4..9b2bc2d8b 100644 --- a/tests/test_autodiffmatrix.cpp +++ b/tests/test_autodiffmatrix.cpp @@ -173,3 +173,51 @@ BOOST_AUTO_TEST_CASE(AdditionOps) BOOST_CHECK(x == ss); } +BOOST_AUTO_TEST_CASE(MultOps) +{ + // Setup. + Mat z = Mat(AutoDiffMatrix::ZeroMatrix, 3); + Sp zs(3,3); + + Mat i = Mat(AutoDiffMatrix::IdentityMatrix, 3); + Sp is(Eigen::Matrix::Identity(3,3).sparseView()); + + Eigen::Array d1(3); + d1 << 0.2, 1.2, 13.4; + Mat d = Mat(d1.matrix().asDiagonal()); + Sp ds = spdiag(d1); + + Eigen::Matrix s1(3,3); + s1 << + 1.0, 0.0, 2.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 2.0; + Sp ss(s1.sparseView()); + Mat s = Mat(ss); + + // Convert to Eigen::SparseMatrix + Sp x; + z.toSparse(x); + BOOST_CHECK(x == zs); + i.toSparse(x); + BOOST_CHECK(x == is); + d.toSparse(x); + BOOST_CHECK(x == ds); + s.toSparse(x); + BOOST_CHECK(x == ss); + + // Multiply by diagonal matrix. + auto ztd = z * d; + ztd.toSparse(x); + BOOST_CHECK(x == zs*ds); + auto itd = i * d; + itd.toSparse(x); + BOOST_CHECK(x == is*ds); + auto dtd = d * d; + dtd.toSparse(x); + BOOST_CHECK(x == ds*ds); + auto std = s * d; + std.toSparse(x); + BOOST_CHECK(x == ss*ds); +} +