Store a copy of the right-hand side

dune-istl's AMG wants a non-const right-hand side to store the
residual, however to conform to the interface we can only get a const
reference passed in.  So we use a non-const copy instead.

dune/solvers/iterationsteps/amgstep.hh

@@ -47,7 +47,7 @@ public:
 
         fop_ = std::unique_ptr<Operator>(new Operator(*stiffnessMatrix));
         amg_ = std::unique_ptr<AMG>(new AMG(*fop_, coarseningCriterion, smootherArgs, 1, 1, 1, false));
-        amg_->pre(*this->x_,*this->rhs_);
+        amg_->pre(*this->x_, residual_);
     }
 
     /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
@@ -66,6 +66,7 @@ public:
         LinearIterationStep<MatrixType, VectorType>::setProblem(stiffnessMatrix, x, rhs);
 
         fop_ = std::unique_ptr<Operator>(new Operator(stiffnessMatrix));
+        residual_ = rhs;
     }
 
     /** \brief Sets up an algebraic hierarchy
@@ -91,19 +92,21 @@ private:
 
     /** \brief The dune-istl AMG step */
     std::unique_ptr<AMG> amg_;
+
+    VectorType residual_;
 };
 
 template <class MatrixType, class VectorType>
 void AMGStep<MatrixType,VectorType>::preprocess()
 {
     amg_ = std::unique_ptr<AMG>(new AMG(*fop_, coarseningCriterion_, smootherArgs_, 1, 1, 1, false));
-    amg_->pre(*this->x_,*this->rhs_);
+    amg_->pre(*this->x_, residual_);
 }
 
 template <class MatrixType, class VectorType>
 void AMGStep<MatrixType,VectorType>::iterate()
 {
-    amg_->apply(*this->x_,*this->rhs_);
+    amg_->apply(*this->x_, residual_);
 }
 
 template <class MatrixType, class VectorType>