diff --git a/dune/solvers/iterationsteps/Makefile.am b/dune/solvers/iterationsteps/Makefile.am index aac9853c9a5c3026577394e6dccad78064518865..050ac67ee659965850172392b7e4138f2e0c2fb1 100644 --- a/dune/solvers/iterationsteps/Makefile.am +++ b/dune/solvers/iterationsteps/Makefile.am @@ -4,6 +4,8 @@ iterationstepsdir = $(includedir)/dune/solvers/iterationsteps iterationsteps_HEADERS = amgstep.hh \ blockgsstep.hh \ blockgsstep.cc \ + cgstep.hh \ + cgstep.cc \ iterationstep.hh \ lineariterationstep.hh \ linegsstep.hh \ diff --git a/dune/solvers/iterationsteps/cgstep.cc b/dune/solvers/iterationsteps/cgstep.cc new file mode 100644 index 0000000000000000000000000000000000000000..d4e61031e4819efe79cc3a399bc0298f927b7d5c --- /dev/null +++ b/dune/solvers/iterationsteps/cgstep.cc @@ -0,0 +1,25 @@ +template <class MatrixType, class VectorType> +void CGStep<MatrixType, VectorType>::preprocess() +{ + // Compute the residual (r starts out as the rhs) + matrix_.mmv(x_,r_); + p_ = r_; + r_squared_old = p_*r_; +} + +template <class MatrixType, class VectorType> +void CGStep<MatrixType, VectorType>::iterate() +{ + 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 + r_.axpy(-alpha, q); // r_1 = r_0 - alpha_0 Ap_0 + + const double r_squared = r_ * r_; + const double beta = r_squared / r_squared_old; // beta_0 = r_1*r_1/ (r_0*r_0) + p_ *= beta; // p_1 = r_1 + beta_0 p_0 + p_ += r_; + r_squared_old = r_squared; +} diff --git a/dune/solvers/iterationsteps/cgstep.hh b/dune/solvers/iterationsteps/cgstep.hh new file mode 100644 index 0000000000000000000000000000000000000000..77e74f58465de6874f3fc5d1ad2a5a5756a0f1b9 --- /dev/null +++ b/dune/solvers/iterationsteps/cgstep.hh @@ -0,0 +1,37 @@ +#ifndef DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH +#define DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH + +#include <dune/solvers/iterationsteps/iterationstep.hh> + +namespace Dune { + namespace Solvers { + //! A conjugate gradient solver + template <class MatrixType, class VectorType> + class CGStep : public IterationStep<VectorType> + { + public: + CGStep(const MatrixType& matrix, + VectorType& x, + const VectorType& rhs) + : p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix) + {} + + void check() const; + void preprocess(); + void iterate(); + // Q: do we really want this interface? + VectorType getSol() { return x_; } + + private: + VectorType p_; // search direction + VectorType r_; // residual + VectorType& x_; + const MatrixType& matrix_; + double r_squared_old; + }; + +#include "cgstep.cc" + } +} + +#endif