Use x_ from the base class

dune/solvers/iterationsteps/cgstep.cc

@@ -12,7 +12,7 @@ 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_);
@@ -26,11 +26,11 @@ void CGStep<MatrixType, VectorType>::preprocess()
 template <class MatrixType, class VectorType>
 void CGStep<MatrixType, VectorType>::iterate()
 {
-    VectorType q(x_);
+    VectorType q(*x_);
 
     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
+    x_->axpy(alpha, p_);                            // x_1     = x_0 + alpha_0 p_0
     r_.axpy(-alpha, q);                            // r_1     = r_0 - alpha_0 Ap_0
 
     if (preconditioner_)

dune/solvers/iterationsteps/cgstep.hh

@@ -13,6 +13,8 @@ namespace Dune {
         template <class MatrixType, class VectorType>
         class CGStep : public IterationStep<VectorType>
         {
+            using Base = IterationStep<VectorType>;
+
         public:
             CGStep()
                 : preconditioner_(nullptr)
@@ -21,7 +23,7 @@ namespace Dune {
             CGStep(const MatrixType& matrix,
                    VectorType& x,
                    const VectorType& rhs)
-                : p_(rhs.size()), r_(rhs), x_(x), matrix_(stackobject_to_shared_ptr(matrix)),
+                : Base(x), p_(rhs.size()), r_(rhs), matrix_(stackobject_to_shared_ptr(matrix)),
                   preconditioner_(nullptr)
             {}
 
@@ -29,7 +31,7 @@ namespace Dune {
                    VectorType& x,
                    const VectorType& rhs,
                    Preconditioner<MatrixType, VectorType>& preconditioner)
-                : p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix),
+                : Base(x), p_(rhs.size()), r_(rhs), matrix_(matrix),
                   preconditioner_(&preconditioner)
             {}
 
@@ -37,7 +39,7 @@ namespace Dune {
             virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs)
             {
                 matrix_ = stackobject_to_shared_ptr(mat);
-                x_ = x;
+                x_ = &x;
                 r_ = rhs;
             }
 
@@ -45,12 +47,12 @@ namespace Dune {
             void preprocess();
             void iterate();
             // Q: do we really want this interface?
-            VectorType getSol() { return x_; }
+            VectorType getSol() { return *x_; }
 
         private:
             VectorType p_; // search direction
             VectorType r_; // residual
-            VectorType& x_;
+            using Base::x_;
             std::shared_ptr<const MatrixType> matrix_;
             double r_squared_old_;
             Preconditioner<MatrixType, VectorType>* preconditioner_;