diff --git a/dune/solvers/iterationsteps/cgstep.cc b/dune/solvers/iterationsteps/cgstep.cc index aa2264c1c6bae650457a6cb5ca9e02e625f68a2c..9dae98fdfb4d962627ee3fa86dd972774904d55a 100644 --- a/dune/solvers/iterationsteps/cgstep.cc +++ b/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 diff --git a/dune/solvers/iterationsteps/cgstep.hh b/dune/solvers/iterationsteps/cgstep.hh index b794d896f72849fcc666e92f3eb2c5dd842578b1..51eb770da2d343b27e31b08b53579f6eb1e26a24 100644 --- a/dune/solvers/iterationsteps/cgstep.hh +++ b/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; }; diff --git a/dune/solvers/iterationsteps/lineariterationstep.hh b/dune/solvers/iterationsteps/lineariterationstep.hh index 55a4fab7a0e5616e9d63f3d42b71bbc0983e7107..cef10b1d86d6d136f17847aa09badf26faa6b738 100644 --- a/dune/solvers/iterationsteps/lineariterationstep.hh +++ b/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