From c3dda5c6eb18e113909f7164e2856885cc3996e8 Mon Sep 17 00:00:00 2001 From: Oliver Sander <oliver.sander@tu-dresden.de> Date: Tue, 17 Jan 2017 16:34:50 +0100 Subject: [PATCH] Make CGStep inherit from LinearIterationStep Because, in view of the recent precise definition of LinearIterationStep, it really is one. --- dune/solvers/iterationsteps/cgstep.cc | 6 +++--- dune/solvers/iterationsteps/cgstep.hh | 19 +++++++++---------- .../iterationsteps/lineariterationstep.hh | 15 ++++++++++----- 3 files changed, 22 insertions(+), 18 deletions(-) diff --git a/dune/solvers/iterationsteps/cgstep.cc b/dune/solvers/iterationsteps/cgstep.cc index aa2264c1..9dae98fd 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 b794d896..51eb770d 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 55a4fab7..cef10b1d 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 -- GitLab