]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
precondition now a class of it's own
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 26 Apr 1999 19:43:05 +0000 (19:43 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 26 Apr 1999 19:43:05 +0000 (19:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@1207 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/precondition.h [new file with mode: 0644]
deal.II/lac/include/lac/solver.h
deal.II/lac/include/lac/solver_bicgstab.h
deal.II/lac/include/lac/solver_cg.h
deal.II/lac/include/lac/solver_gauss_seidel.h [deleted file]
deal.II/lac/include/lac/solver_gmres.h
deal.II/lac/include/lac/solver_richardson.h

diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h
new file mode 100644 (file)
index 0000000..afc432e
--- /dev/null
@@ -0,0 +1,165 @@
+/*-------------------- precondition.h --------------------*/
+//$Id$
+// Guido Kanschat, University of Heidelberg, 1999
+
+#ifndef __lac_precondition_h
+#define __lac_precondition_h
+
+/**
+ * No preconditioning.
+ * @author Guido Kanschat, 1999
+ */
+template<class VECTOR>
+class PreconditionIdentity
+{
+  public:
+                                    /**
+                                     * Execute preconditioning.
+                                     */
+    void operator() (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * Preconditioner using a matrix-builtin function.
+ * This class forms a preconditioner suitable for the LAC solver
+ * classes. Since many preconditioning methods are based on matrix
+ * entries, these have to be implemented as member functions of the
+ * underlying matrix implementation. This class now is intended to
+ * allow easy access to these member functions from LAC solver
+ * classes.
+ *
+ * It seems that all builtin preconditioners have a relaxation
+ * parameter, so please use #PreconditionRelaxation# for these.
+ * @author Guido Kanschat, 1999
+ */ 
+
+template<class MATRIX, class VECTOR>
+class PreconditionUseMatrix
+{
+  public:
+                                    /**
+                                     * Type of the preconditioning
+                                     * function of the matrix.
+                                     */
+    typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&) const;
+    
+                                    /**
+                                     * Constructor.
+                                     * This constructor stores a
+                                     * reference to the matrix object
+                                     * for later use and selects a
+                                     * preconditioning method, which
+                                     * must be a member function of
+                                     * that matrix.
+                                     */
+    PreconditionUseMatrix(const MATRIX & M, function_ptr method);
+    
+                                    /**
+                                     * Execute preconditioning.
+                                     */
+    void operator() (VECTOR&, const VECTOR&) const;
+
+  private:
+                                    /**
+                                     * Pointer to the matrix in use.
+                                     */
+    const MATRIX& matrix;
+                                    /**
+                                     * Pointer to the preconditioning
+                                     * function.
+                                     */
+    const function_ptr precondition;
+};
+
+/**
+ * Preconditioner for builtin relaxation methods.
+ * Uses methods of #SparseMatrix#.
+ * @author Guido Kanschat, 1999
+ */
+template<class MATRIX, class VECTOR>
+class PreconditionRelaxation
+{
+  public:
+                                    /**
+                                     * Type of the preconditioning
+                                     * function of the matrix.
+                                     */
+    typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&,
+                                          typename MATRIX::entry_type) const;
+    
+                                    /**
+                                     * Constructor.
+                                     * This constructor stores a
+                                     * reference to the matrix object
+                                     * for later use and selects a
+                                     * preconditioning method, which
+                                     * must be a member function of
+                                     * that matrix.
+                                     */
+    PreconditionRelaxation(const MATRIX & M, function_ptr method,
+                          double omega = 1.);
+    
+                                    /**
+                                     * Execute preconditioning.
+                                     */
+    void operator() (VECTOR&, const VECTOR&) const;
+
+  private:
+                                    /**
+                                     * Pointer to the matrix in use.
+                                     */
+    const MATRIX& matrix;
+                                    /**
+                                     * Pointer to the preconditioning
+                                     * function.
+                                     */
+    const function_ptr precondition;
+                                    /**
+                                     * Relaxation parameter.
+                                     */
+    double omega;
+};
+
+
+template<class VECTOR>
+void
+PreconditionIdentity<VECTOR>::operator() (VECTOR& dst, const VECTOR& src) const
+{
+  dst = src;
+}
+
+template<class MATRIX, class VECTOR>
+PreconditionUseMatrix<MATRIX, VECTOR>::PreconditionUseMatrix(const MATRIX& M,
+                                                            function_ptr method)
+               :
+               matrix(M), precondition(method)
+{}
+
+
+template<class MATRIX, class VECTOR>
+void
+PreconditionUseMatrix<MATRIX, VECTOR>::operator() (VECTOR& dst,
+                                                  const VECTOR& src) const
+{
+  (matrix.*precondition)(dst, src);
+}
+
+template<class MATRIX, class VECTOR>
+PreconditionRelaxation<MATRIX, VECTOR>::PreconditionRelaxation(const MATRIX& M,
+                                                            function_ptr method,
+                                                              double omega)
+               :
+               matrix(M), precondition(method), omega(omega)
+{}
+
+
+template<class MATRIX, class VECTOR>
+void
+PreconditionRelaxation<MATRIX, VECTOR>::operator() (VECTOR& dst,
+                                                  const VECTOR& src) const
+{
+  (matrix.*precondition)(dst, src, omega);
+}
+
+#endif
index 25fa6d4a73c192ed860cd629898ed58d8d2d3666..a57b1567a738625ac23f384edae6134ba3554093 100644 (file)
@@ -13,7 +13,11 @@ template <class Vector> class VectorMemory;
 
 
 /**
- * Base class for iterative solvers. This class defines possible
+ * Base class for iterative solvers.
+ *
+ * HAS TO BE UPDATED!
+ *
+ * This class defines possible
  * return states of linear solvers and provides interfaces to a memory
  * pool and the control object.
  *
@@ -115,20 +119,18 @@ class Solver
                                      * object.
                                      */
     Solver (SolverControl &, VectorMemory<Vector> &);
-    
-                                    /**
-                                     * Destructor.  This one is virtual,
-                                     * overload it in derived classes if
-                                     * necessary.
-                                     */
-    virtual ~Solver() {}
-    
+
                                     /**
-                                     * Solver procedure.
+                                     * Solver procedure.  This is in
+                                     * fact only a template for the
+                                     * solve-function to be
+                                     * implemented in derived classes
                                      */
-    virtual ReturnState solve (const Matrix &A,
-                              Vector &x,
-                              const Vector &b) = 0;
+    template<class Preconditioner>
+    ReturnState solve (const Matrix &A,
+                      Vector &x,
+                      const Vector &b,
+                      const Preconditioner& precondition) const;
     
                                     /**
                                      * Access to control object.
index ae7cae0a1f1f34811f4fa640f193c495a92e865c..1b6b5e1da6043f2f559e55a6ce1c8f26a3696744 100644 (file)
@@ -25,9 +25,11 @@ class SolverBicgstab //<Matrix,Vector>
                                     /**
                                      * Solve primal problem only.
                                      */
-    virtual ReturnState solve (const Matrix &A,
-                              Vector       &x,
-                              const Vector &b);
+    template<class Preconditioner>
+    ReturnState solve (const Matrix &A,
+                      Vector       &x,
+                      const Vector &b,
+                      const Preconditioner& precondition);
 
   protected:
                                     /**
@@ -121,7 +123,8 @@ class SolverBicgstab //<Matrix,Vector>
                                     /**
                                      * The iteration loop itself.
                                      */
-    ReturnState iterate();
+    template<class Preconditioner>
+    ReturnState iterate(const Preconditioner& precondition);
   
 };
 
@@ -135,7 +138,8 @@ SolverBicgstab<Matrix, Vector>::SolverBicgstab(SolverControl &cn,
 {}
 
 
-template<class Matrix, class Vector> double
+template<class Matrix, class Vector>
+double
 SolverBicgstab<Matrix, Vector>::criterion()
 {
   res = MA->residual(*Vt, *Vx, *Vb);
@@ -143,7 +147,8 @@ SolverBicgstab<Matrix, Vector>::criterion()
 }
 
 
-template < class Matrix, class Vector > SolverControl::State
+template < class Matrix, class Vector >
+SolverControl::State
 SolverBicgstab<Matrix, Vector>::start()
 {
   res = MA->residual(*Vr, *Vx, *Vb);
@@ -154,8 +159,10 @@ SolverBicgstab<Matrix, Vector>::start()
   return state;
 }
 
-template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::iterate()
+template<class Matrix, class Vector>
+template<class Preconditioner>
+Solver<Matrix,Vector>::ReturnState
+SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
 {
   SolverControl::State state = SolverControl::iterate;
   alpha = omega = rho = 1.;
@@ -175,7 +182,7 @@ SolverBicgstab<Matrix, Vector>::iterate()
       beta   = rhobar * alpha / (rho * omega);
       rho    = rhobar;
       p.sadd(beta, 1., r, -beta*omega, v);
-      MA->precondition(y,p);
+      precondition(y,p);
       MA->vmult(v,y);
       rhobar = rbar * v;
 
@@ -183,7 +190,7 @@ SolverBicgstab<Matrix, Vector>::iterate()
     
       alpha = rho/rhobar;
       s.equ(1., r, -alpha, v);
-      MA->precondition(z,s);
+      precondition(z,s);
       MA->vmult(t,z);
       rhobar = t*s;
       omega = rhobar/(t*t);
@@ -198,10 +205,13 @@ SolverBicgstab<Matrix, Vector>::iterate()
 }
 
 
-template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
+template<class Matrix, class Vector>
+template<class Preconditioner>
+Solver<Matrix,Vector>::ReturnState
 SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
                                      Vector       &x,
-                                     const Vector &b)
+                                     const Vector &b,
+                                     const Preconditioner& precondition)
 {
   deallog.push("Bicgstab");
   Vr    = memory.alloc(); Vr->reinit(x);
@@ -225,7 +235,7 @@ SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
     {
       deallog << "Go!" << endl;
       if (start() == SolverControl::success) break;  
-      state = iterate();
+      state = iterate(precondition);
     }
   while (state == breakdown);
 
index 5cd18d96412888a4905e3a7319fd715b2dba09a4..30d2ff731c16fc27905deaac139fe2b39c0363c5 100644 (file)
@@ -17,7 +17,8 @@
  * @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and  documentation by Wolfgang Bangerth
  */
 template<class Matrix, class Vector>
-class SolverPCG : public Solver<Matrix,Vector> {
+class SolverPCG : public Solver<Matrix,Vector>
+{
   public:
 
                                     /**
@@ -29,9 +30,11 @@ class SolverPCG : public Solver<Matrix,Vector> {
                                     /**
                                      * Solver method.
                                      */
-    virtual ReturnState solve (const Matrix &A,
-                              Vector       &x,
-                              const Vector &b);
+    template<class Preconditioner>
+    ReturnState solve (const Matrix &A,
+                      Vector       &x,
+                      const Vector &b,
+                      const Preconditioner& precondition);
 
   protected:
                                     /**
@@ -76,10 +79,13 @@ double SolverPCG<Matrix,Vector>::criterion()
 
 
 template<class Matrix, class Vector>
+template<class Preconditioner>
 Solver<Matrix,Vector>::ReturnState 
 SolverPCG<Matrix,Vector>::solve (const Matrix &A,
                                 Vector       &x,
-                                const Vector &b) {
+                                const Vector &b,
+                                const Preconditioner& precondition)
+{
   SolverControl::State conv=SolverControl::iterate;
   
                                   // Memory allocation
@@ -117,7 +123,7 @@ SolverPCG<Matrix,Vector>::solve (const Matrix &A,
     };
   
   g.scale(-1.);
-  A.precondition(h,g);
+  precondition(h,g);
  
   d.equ(-1.,h);
  
@@ -137,7 +143,7 @@ SolverPCG<Matrix,Vector>::solve (const Matrix &A,
       conv = control().check(it,res);
       if (conv) break;
       
-      A.precondition(h,g);
+      precondition(h,g);
       
       beta = gh;
       gh   = g*h;
diff --git a/deal.II/lac/include/lac/solver_gauss_seidel.h b/deal.II/lac/include/lac/solver_gauss_seidel.h
deleted file mode 100644 (file)
index 6b37608..0000000
+++ /dev/null
@@ -1,132 +0,0 @@
-/*----------------------------   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("lac/solver_gauss_seidel.h");
-  deallog.push("SolverGaussSeidel::solve");
-
-  unsigned int i;
-  unsigned int n = b.size();
-  double res;
-
-  SolverControl::State conv=SolverControl::iterate;
-
-  Vector defect(n);
-
-  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     ---------------------------*/
index 5db3b45de1e61ef4aadf8c076d798ab49804c496..76f138924ea527b4d5dce46121cdb400fa8ee677 100644 (file)
  * computational speed against memory. One of the few other
  * possibilities is to use a good preconditioner.
  *
- * @author Original implementation by the DEAL authors; adapted, cleaned and documented by Wolfgang Bangerth */
+ * @author Wolfgang Bangerth
+ */
 template<class Matrix, class Vector>
-class SolverGMRES : public Solver<Matrix, Vector> {
+class SolverGMRES : public Solver<Matrix, Vector>
+{
   public:
                                     /**
                                      * Constructor.
@@ -59,9 +61,11 @@ class SolverGMRES : public Solver<Matrix, Vector> {
                                     /**
                                      * Solver method.
                                      */
-    virtual ReturnState solve (const Matrix &A,
-                              Vector       &x,
-                              const Vector &b);
+    template<class Preconditioner>
+    ReturnState solve (const Matrix &A,
+                      Vector       &x,
+                      const Vector &b,
+                      const Preconditioner& precondition);
 
     DeclException1 (ExcTooFewTmpVectors,
                    int,
@@ -137,10 +141,12 @@ SolverGMRES<Matrix,Vector>::givens_rotation (Vector &h,
 
 
 template<class Matrix, class Vector>
+template<class Preconditioner>
 Solver<Matrix,Vector>::ReturnState
 SolverGMRES<Matrix,Vector>::solve (const Matrix& A,
                                   Vector      & x,
-                                  const Vector& b)
+                                  const Vector& b,
+                                  const Preconditioner& precondition)
 {
                                   // this code was written by the fathers of
                                   // DEAL. I take absolutely no guarantees
@@ -219,7 +225,7 @@ SolverGMRES<Matrix,Vector>::solve (const Matrix& A,
       if (left_precondition)
        {
          A.residual(p,x,b);
-         A.precondition(v,p);
+         precondition(v,p);
        } else {
          A.residual(v,x,b);
        };
@@ -259,9 +265,9 @@ SolverGMRES<Matrix,Vector>::solve (const Matrix& A,
          if (left_precondition)
            {
              A.vmult(p, *tmp_vectors[inner_iteration]);
-             A.precondition(vv,p);
+             precondition(vv,p);
            } else {
-             A.precondition(p,*tmp_vectors[inner_iteration]);
+             precondition(p,*tmp_vectors[inner_iteration]);
              A.vmult(vv,p);
            };
          
@@ -324,7 +330,7 @@ SolverGMRES<Matrix,Vector>::solve (const Matrix& A,
          p = 0.;
          for (unsigned int i=0; i<dim; ++i)
            p.add(h(i), *tmp_vectors[i]);
-         A.precondition(v,p);
+         precondition(v,p);
          x.add(1.,v);
        };
 
index 877f3f2d357eb83443d5a980e31d95676a458d21..2ccca467a8cf27492f0edd84d48bc8f61da5aed9 100644 (file)
@@ -31,9 +31,11 @@ class SolverRichardson : public Solver<Matrix, Vector> {
                                     /**
                                      * Solve $Ax=b$ for $x$.
                                      */
-    virtual ReturnState solve (const Matrix &A,
-                              Vector       &x,
-                              const Vector &b);
+    template<class Preconditioner>
+    ReturnState solve (const Matrix &A,
+                      Vector       &x,
+                      const Vector &b,
+                      const Preconditioner& precondition);
 
                                     /**
                                      * Set the damping-coefficient.
@@ -81,10 +83,13 @@ class SolverRichardson : public Solver<Matrix, Vector> {
 
 
 template<class Matrix,class Vector>
+template<class Preconditioner>
 Solver<Matrix,Vector>::ReturnState 
 SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
-                               Vector       &x,
-                               const Vector &b) {
+                                       Vector       &x,
+                                       const Vector &b,
+                                       const Preconditioner& precondition)
+{
   SolverControl::State conv=SolverControl::iterate;
 
                                   // Memory allocation
@@ -103,7 +108,7 @@ SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
       if (conv != SolverControl::iterate)
        break;
 
-      A.precondition(d,r);
+      precondition(d,r);
       x.add(omega,d);
     }
 

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.