From 64b17e95105e70b5c0c47b38251c41a2dd287e1f Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 30 Aug 2013 13:04:13 +0000
Subject: [PATCH] Add Dune::Solvers::CG (no preconditioning)
[[Imported from SVN: r11957]]
---
dune/solvers/iterationsteps/Makefile.am | 2 ++
dune/solvers/iterationsteps/cgstep.cc | 25 +++++++++++++++++
dune/solvers/iterationsteps/cgstep.hh | 37 +++++++++++++++++++++++++
3 files changed, 64 insertions(+)
create mode 100644 dune/solvers/iterationsteps/cgstep.cc
create mode 100644 dune/solvers/iterationsteps/cgstep.hh
diff --git a/dune/solvers/iterationsteps/Makefile.am b/dune/solvers/iterationsteps/Makefile.am
index aac9853c..050ac67e 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 00000000..d4e61031
--- /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 00000000..77e74f58
--- /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
--
GitLab