Make CGStep inherit from LinearIterationStep

Because, in view of the recent precise definition of
LinearIterationStep, it really is one.

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_);
+    this->mat_->mmv(*x_,r_);
     if (this->hasIgnore())
         for (size_t i=0; i < r_.size(); ++i)
             if (this->ignore()[i].any())
@@ -24,7 +24,7 @@ void CGStep<MatrixType, VectorType>::preprocess()
         if (asCanIgnore != nullptr and this->hasIgnore())
             asCanIgnore->setIgnore(this->ignore());
 
-        preconditioner_->setMatrix(*matrix_);
+        preconditioner_->setMatrix(*this->mat_);
         preconditioner_->apply(p_, r_);
     } else
         p_ = r_;
@@ -41,7 +41,7 @@ void CGStep<MatrixType, VectorType>::iterate()
 
     VectorType q(*x_);
 
-    matrix_->mv(p_, q);                             // q_0     = Ap_0
+    this->mat_->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

@@ -11,9 +11,9 @@ namespace Dune {
 
         //! A conjugate gradient solver
         template <class MatrixType, class VectorType>
-        class CGStep : public IterationStep<VectorType>
+        class CGStep : public LinearIterationStep<MatrixType,VectorType>
         {
-            using Base = IterationStep<VectorType>;
+            using Base = LinearIterationStep<MatrixType,VectorType>;
 
         public:
             CGStep()
@@ -23,7 +23,7 @@ namespace Dune {
             CGStep(const MatrixType& matrix,
                    VectorType& x,
                    const VectorType& rhs)
-                : Base(x), p_(rhs.size()), r_(rhs), matrix_(stackobject_to_shared_ptr(matrix)),
+                : Base(matrix,x), p_(rhs.size()), r_(rhs),
                   preconditioner_(nullptr)
             {}
 
@@ -31,31 +31,30 @@ namespace Dune {
                    VectorType& x,
                    const VectorType& rhs,
                    Preconditioner<MatrixType, VectorType>& preconditioner)
-                : Base(x), p_(rhs.size()), r_(rhs), matrix_(stackobject_to_shared_ptr(matrix)),
+                : Base(matrix,x), p_(rhs.size()), r_(rhs),
                   preconditioner_(&preconditioner)
             {}
 
             //! Set linear operator, solution and right hand side
-            virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs)
+            virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override
             {
-                matrix_ = stackobject_to_shared_ptr(mat);
+                this->setMatrix(mat);
                 setProblem(x);
                 r_ = rhs;
             }
 
-            void check() const;
+            void check() const override;
             void preprocess();
             void iterate();
 
+            using Base::setProblem;
+
         private:
             VectorType p_; // search direction
             VectorType r_; // residual
             using Base::x_;
-            std::shared_ptr<const MatrixType> matrix_;
             double r_squared_old_;
             Preconditioner<MatrixType, VectorType>* preconditioner_;
-
-            using Base::setProblem;
         };
 
 

dune/solvers/iterationsteps/lineariterationstep.hh

@@ -35,6 +35,7 @@ class LinearIterationStep : public IterationStep<VectorType, BitVectorType>,
                                                                  VectorType,
                                                                  BitVectorType>
 {
+    using Base = IterationStep<VectorType, BitVectorType>;
 public:
 
     //! Default constructor
@@ -45,11 +46,17 @@ public:
     /** \brief Destructor */
     virtual ~LinearIterationStep() {}
 
-    //! Constructor being given linear operator, solution and right hand side
+    //! Constructor with given matrix, solution and right hand side
     LinearIterationStep(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
       setProblem(mat, x, rhs);
     }
 
+    //! Constructor with given matrix and solution vector
+    LinearIterationStep(const MatrixType& matrix, VectorType& x)
+    : IterationStep<VectorType,BitVectorType>(x),
+      mat_(stackobject_to_shared_ptr(matrix))
+    {}
+
     //! Set linear operator, solution and right hand side
     virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
         setProblem(x);
@@ -92,15 +99,13 @@ public:
         x = this->getSol();
     }
 
+    using Base::setProblem;
+
     //! The container for the right hand side
     const VectorType* rhs_;
 
     //! The linear operator
     std::shared_ptr<const MatrixType> mat_;
-
-private:
-    using Base = IterationStep<VectorType, BitVectorType>;
-    using Base::setProblem;
 };
 
 }   // namespace Solvers