]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Richardson's iteration method as preconditioned solver.
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Apr 1999 11:57:45 +0000 (11:57 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Apr 1999 11:57:45 +0000 (11:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@1177 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver_richardson.h [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h
new file mode 100644 (file)
index 0000000..877f3f2
--- /dev/null
@@ -0,0 +1,141 @@
+/*----------------------------   solver_richardson.h     ---------------------------*/
+/*   $Id$              */
+/*            Ralf Hartmann, University of Heidelberg                               */
+#ifndef __solver_richardson_H
+#define __solver_richardson_H
+/*----------------------------   solver_richardson.h     ---------------------------*/
+
+
+
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+
+
+
+/**
+ * Implementation of the richardson iteration method. The stopping criterion
+ * is the norm of the residual.
+ *
+ * @author Ralf Hartmann
+ */
+template<class Matrix, class Vector>
+class SolverRichardson : public Solver<Matrix, Vector> {
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverRichardson (SolverControl &cn, VectorMemory<Vector> &mem) :
+                   Solver<Matrix,Vector> (cn,mem), omega(1.)
+      {};
+
+                                    /**
+                                     * Solve $Ax=b$ for $x$.
+                                     */
+    virtual ReturnState solve (const Matrix &A,
+                              Vector       &x,
+                              const Vector &b);
+
+                                    /**
+                                     * Set the damping-coefficient.
+                                     * Default is 1., i.e. no damping.
+                                     */
+    void set_omega(double om=1.);
+    
+  protected:
+                                    /**
+                                     * Implementation of the computation of
+                                     * the norm of the residual.
+                                     */
+    virtual double criterion();
+    
+                                    /**
+                                     * Temporary vectors, allocated through
+                                     * the #VectorMemory# object at the start
+                                     * of the actual solution process and
+                                     * deallocated at the end.
+                                     */
+    Vector *Vr, *Vd;
+
+                                    /**
+                                     * Damping-coefficient.
+                                     */
+    double omega;  ;
+    
+                                    /**
+                                     * Within the iteration loop, the
+                                     * square of the residual vector is
+                                     * stored in this variable. The
+                                     * function #criterion# uses this
+                                     * variable to compute the convergence
+                                     * value, which in this class is the
+                                     * norm of the residual vector and thus
+                                     * the square root of the #res2# value.
+                                     */
+    double res2;
+};
+
+
+
+
+/*----------------- Implementation of the Richardson Method ------------------*/
+
+
+template<class Matrix,class Vector>
+Solver<Matrix,Vector>::ReturnState 
+SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
+                               Vector       &x,
+                               const Vector &b) {
+  SolverControl::State conv=SolverControl::iterate;
+
+                                  // Memory allocation
+  Vr  = memory.alloc(); Vector& r  = *Vr; r.reinit(x);
+  Vd  = memory.alloc(); Vector& d  = *Vd; d.reinit(x);
+
+  deallog.push("Richardson");
+
+                                  // Main loop
+  for(int iter=0; conv==SolverControl::iterate; iter++)
+    {
+      A.residual(r,x,b);
+
+      res2 = r*r;
+      conv = control().check (iter, criterion());
+      if (conv != SolverControl::iterate)
+       break;
+
+      A.precondition(d,r);
+      x.add(omega,d);
+    }
+
+                                  // Deallocate Memory
+  memory.free(Vr);
+  memory.free(Vd);
+
+  deallog.pop();
+                                  // Output
+  if (conv == SolverControl::failure)
+    return exceeded;
+  else
+    return success;
+}
+
+
+template<class Matrix,class Vector>
+inline double
+SolverRichardson<Matrix,Vector>::criterion()
+{
+  return sqrt(res2);
+}
+
+template<class Matrix,class Vector>
+inline void
+SolverRichardson<Matrix,Vector>::set_omega(double om)
+{
+  omega=om;
+}
+
+
+/*------------------   solver_richardson.h     ----------------------*/
+/* end of #ifndef __solver_richardson_H */
+#endif
+/*------------------   solver_richardson.h     ----------------------*/

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.