Add Dune::Solvers::CG (no preconditioning)

[[Imported from SVN: r11957]]

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 \

dune/solvers/iterationsteps/cgstep.cc

+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;
+}

dune/solvers/iterationsteps/cgstep.hh

+#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