Store the matrix as a std::shared_ptr, instead of a const reference

That makes the cgsteptest compile and pass again.

dune/solvers/iterationsteps/cgstep.cc

@@ -9,10 +9,10 @@ template <class MatrixType, class VectorType>
 void CGStep<MatrixType, VectorType>::preprocess()
 {
     // Compute the residual (r starts out as the rhs)
-    matrix_.mmv(x_,r_);
+    matrix_->mmv(x_,r_);
 
     if (preconditioner_) {
-        preconditioner_->setMatrix(matrix_);
+        preconditioner_->setMatrix(*matrix_);
         preconditioner_->apply(p_, r_);
     } else
         p_ = r_;
@@ -25,7 +25,7 @@ void CGStep<MatrixType, VectorType>::iterate()
 {
     VectorType q(x_);
 
-    matrix_.mv(p_, q);                             // q_0     = Ap_0
+    matrix_->mv(p_, q);                             // q_0     = Ap_0
     const double alpha = r_squared_old_ / (p_ * q); // alpha_0 = r_0*r_0/p_0*Ap_0
     x_.axpy(alpha, p_);                            // x_1     = x_0 + alpha_0 p_0
     r_.axpy(-alpha, q);                            // r_1     = r_0 - alpha_0 Ap_0

dune/solvers/iterationsteps/cgstep.hh

@@ -18,7 +18,7 @@ namespace Dune {
             CGStep(const MatrixType& matrix,
                    VectorType& x,
                    const VectorType& rhs)
-                : p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix),
+                : p_(rhs.size()), r_(rhs), x_(x), matrix_(Dune::stackobject_to_shared_ptr(matrix)),
                   preconditioner_(nullptr)
             {}
 
@@ -33,7 +33,7 @@ namespace Dune {
             //! Set linear operator, solution and right hand side
             virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs)
             {
-                matrix_ = mat;
+                matrix_ = Dune::stackobject_to_shared_ptr(mat);
                 x_ = x;
                 r_ = rhs;
             }
@@ -48,7 +48,7 @@ namespace Dune {
             VectorType p_; // search direction
             VectorType r_; // residual
             VectorType& x_;
-            const MatrixType& matrix_;
+            std::shared_ptr<const MatrixType> matrix_;
             double r_squared_old_;
             Preconditioner<MatrixType, VectorType>* preconditioner_;
         };