]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implementation of the Gauss-Seidel iterative solver
authornauber <nauber@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Apr 1999 08:39:39 +0000 (08:39 +0000)
committernauber <nauber@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Apr 1999 08:39:39 +0000 (08:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@1184 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/lac/include/lac/solver_gauss_seidel.h b/deal.II/lac/include/lac/solver_gauss_seidel.h
new file mode 100644 (file)
index 0000000..19c8a41
--- /dev/null
@@ -0,0 +1,136 @@
+/*----------------------------   solver_gauss_seidel.h     ---------------------------*/
+/*      $Id$             */
+#ifndef __solver_gauss_seidel_H
+#define __solver_gauss_seidel_H
+/*----------------------------   solver_gauss_seidel.h     ---------------------------*/
+
+#include <base/trace.h>
+#include <lac/solver_control.h>
+#include <lac/solver.h>
+
+/**
+ * Implementation of the Gauss-Seidel method. The stopping criterion
+ * is the norm of the defect, i.e. if x is the iteration vector, b the
+ * rhs and A the matrix, the $res = \| A x - b \|$. 
+ */
+template<class Matrix, class Vector>
+class SolverGaussSeidel : public Solver<Matrix, Vector> {
+ public:
+  /**
+   * Constructor.  Takes a SolverControl and some
+   VectorMemory. VectorMemory is not used. 
+   */
+  SolverGaussSeidel (SolverControl &cn, VectorMemory<Vector> &mem) :
+    Solver<Matrix,Vector> (cn,mem)
+    {};
+  
+  /**
+   * Solver method. Just specify the vectors A, x and b and let it
+   * run. x should contain a reasonable starting value.  
+   */
+    virtual ReturnState solve (const Matrix &A,
+                              Vector       &x,
+                              const Vector &b);
+
+  protected:
+                                    /**
+                                     * Implementation of the computation of
+                                     * the norm of the residual.
+                                     */
+    virtual double criterion();
+    
+                                    /**
+                                     * 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 Gauss-Seidel Method ------------------------*/
+
+template<class Matrix,class Vector> 
+SolverGaussSeidel<Matrix,Vector>::ReturnState 
+SolverGaussSeidel<Matrix,Vector>::solve (const Matrix &A,
+                                        Vector       &x,
+                                        const Vector &b) {
+
+  deallog.push(__FILE__);
+  deallog.push("SolverGaussSeidel::solve");
+
+  TRACEMSG("Initialize Variables");
+
+  unsigned int i;
+  unsigned int n = b.size();
+  double res;
+
+  SolverControl::State conv=SolverControl::iterate;
+
+  Vector defect(n);
+
+  TRACEMSG("Entering main loop");
+
+  deallog.push("Main loop");
+
+  // Main loop
+  for(int iter=0; conv==SolverControl::iterate; iter++)
+    {
+      // Calculate defect, i.e. defect = A x - b
+      A.vmult(defect,x);      
+      defect.add(-1,b);    
+      
+      // Calculate residual
+      res2 = defect * defect;
+     
+      // Apply Gauss-Seidel preconditioner
+      for (i=0;i<n;i++)
+       {
+         defect(i) = defect(i) / A(i,i);
+       }
+
+      // Correct defect
+      x.add(-1,defect);
+
+      // Check residual
+      res = criterion();
+      conv = (control().check(iter, res));
+      if (conv != SolverControl::iterate)
+       break;
+    }
+  
+  // Output
+  
+  deallog.pop();
+
+  if (conv == SolverControl::failure)
+    {
+      TRACEMSG("*** exceeded maximum number of iterations");
+      deallog.pop();
+      deallog.pop();
+      return exceeded;
+    }
+  else
+    {
+      TRACEMSG("success");
+      deallog.pop();
+      deallog.pop();
+      return success;
+    }
+} 
+
+template<class Matrix,class Vector>
+inline double
+SolverGaussSeidel<Matrix,Vector>::criterion()
+{
+  return sqrt(res2);
+}
+
+/*----------------------------   solver_gauss_seidel.h     ---------------------------*/
+/* end of #ifndef __solver_gauss_seidel_H */
+#endif
+/*----------------------------   solver_gauss_seidel.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.