]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new preconditioner interface
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 Jun 2000 20:03:23 +0000 (20:03 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 Jun 2000 20:03:23 +0000 (20:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@3031 0785d39b-7218-0410-832d-ea1e28bc413d

23 files changed:
deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/vectors.cc
deal.II/lac/include/lac/eigen.h
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/precondition.h
deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/solver_bicgstab.h
deal.II/lac/include/lac/solver_cg.h
deal.II/lac/include/lac/solver_gmres.h
deal.II/lac/include/lac/solver_minres.h
deal.II/lac/include/lac/solver_qmrs.h
deal.II/lac/include/lac/solver_richardson.h
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_vanka.h
deal.II/lac/include/lac/sparse_vanka.templates.h
deal.II/lac/source/precondition_block.cc
deal.II/lac/source/sparse_vanka.cc
tests/big-tests/multigrid/multigrid.cc
tests/deal.II/block_matrices.cc
tests/deal.II/mg.cc

index b53705d6b5ebeb72636d81c3181420c5b24cfd8a..0037d7e679a84f7d3665ca7677ec17ec64144a52 100644 (file)
@@ -443,8 +443,8 @@ class PreconditionMG
                                      * This is the operator used by
                                      * LAC iterative solvers.
                                      */
-    void operator() (VECTOR       &dst,
-                    const VECTOR &src) const;
+    void vmult (VECTOR       &dst,
+               const VECTOR &src) const;
     
   private:
                                     /**
@@ -583,8 +583,8 @@ PreconditionMG<MG, VECTOR>::PreconditionMG(MG                       &mg,
 
 template<class MG, class VECTOR>
 void
-PreconditionMG<MG,VECTOR>::operator() (VECTOR       &dst,
-                                      const VECTOR &src) const
+PreconditionMG<MG,VECTOR>::vmult (VECTOR       &dst,
+                                 const VECTOR &src) const
 {
   multigrid->copy_to_mg(src);
   multigrid->vcycle(*pre, *post, *coarse);
index 523441d82c7b7d2cf2506a4eea98625f9b4f2991..b3281538f3814949314f1bef257bb7cf7a2b5eb5 100644 (file)
@@ -154,10 +154,8 @@ void ProblemBase<dim>::solve ()
   PrimitiveVectorMemory<> memory;
   SolverCG<>              cg(control,memory);
 
-  PreconditionRelaxation<>
-    prec(system_matrix,
-        &SparseMatrix<double>::template precondition_SSOR<double>,
-        1.2);
+  PreconditionSSOR<> prec;
+  prec.initialize (system_matrix, 1.2);
 
                                   // solve
   cg.solve (system_matrix, solution, right_hand_side, prec);
index f5728bf73735616dc4351192d0a260db06cae4ee..2d2062a377080cb1642fcc6eb5cdfbc82819b1b4 100644 (file)
@@ -444,10 +444,8 @@ void VectorTools::project (const DoFHandler<dim>    &dof,
   PrimitiveVectorMemory<> memory;
   SolverCG<>              cg(control,memory);
 
-  PreconditionRelaxation<>
-    prec(mass_matrix,
-        &SparseMatrix<double>::template precondition_SSOR<double>,
-        1.2);
+  PreconditionSSOR<> prec;
+  prec.initialize(mass_matrix, 1.2);
                                   // solve
   cg.solve (mass_matrix, vec, tmp, prec);
   
@@ -763,10 +761,8 @@ VectorTools::project_boundary_values (const DoFHandler<dim>    &dof,
   PrimitiveVectorMemory<> memory;
   SolverCG<>              cg(control,memory);
 
-  PreconditionRelaxation<>
-    prec(mass_matrix,
-        &SparseMatrix<double>::template precondition_SSOR<double>,
-        1.2);
+  PreconditionSSOR<> prec;
+  prec.initialize(mass_matrix, 1.2);
                                   // solve
   cg.solve (mass_matrix, boundary_projection, rhs, prec);
 
index 6293e808632d84a4c89fcaca326a798fe1b4f7dd..5e1ef124ae053db2b243e305cba6d84385dc59a9 100644 (file)
@@ -357,6 +357,7 @@ EigenPower<VECTOR>::solve (double       &value,
     }
   
   memory.free(Vy);
+  memory.free(Vr);
 
   deallog.pop();
                                   // Output
@@ -465,7 +466,8 @@ EigenInverse<VECTOR>::solve (double       &value,
     }
 
   memory.free(Vy);
-
+  memory.free(Vr);
+  
   deallog.pop();
 
                                   // Output
index e44678438a89c224a1b14e25ecbfd9078734fae1..0a4a9bf4b6c834cfe69b0daf70d521bfd4ba7096 100644 (file)
@@ -560,13 +560,6 @@ class FullMatrix : public Subscriptor
                                     /**
                                      * Exception
                                      */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The given index " << arg1
-                   << " should be less than " << arg2 << ".");
-                                    /**
-                                     * Exception
-                                     */
     DeclException2 (ExcDimensionMismatch,
                    int, int,
                    << "The two dimensions " << arg1 << " and " << arg2
@@ -699,8 +692,8 @@ template <typename number>
 inline number
 FullMatrix<number>::operator() (const unsigned int i, const unsigned int j) const
 {  
-  Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
-  Assert (j<dim_range, ExcInvalidIndex (j, dim_range));
+  Assert (i<dim_image, ExcIndexRange (i, 0, dim_image));
+  Assert (j<dim_range, ExcIndexRange (j, 0, dim_range));
   return el(i,j);
 };
 
@@ -710,8 +703,8 @@ template <typename number>
 inline number &
 FullMatrix<number>::operator() (const unsigned int i, const unsigned int j)
 {
-  Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
-  Assert (j<dim_range, ExcInvalidIndex (j, dim_range));
+  Assert (i<dim_image, ExcIndexRange (i, 0, dim_image));
+  Assert (j<dim_range, ExcIndexRange (j, 0, dim_range));
   return el(i,j);
 };
 
index 02dd2809c5c201daa8eb4433ff6365b1d93fd2c1..9a2f13c4939299bedd4902fa35f9f37437d714ab 100644 (file)
@@ -24,7 +24,6 @@
 #include <iomanip>
 #include <algorithm>
 
-
 template <typename number>
 FullMatrix<number>::FullMatrix (const unsigned int n) :
                val (0),
index 97e65c802006a961971b15c28b24f9e734cde370..410465e54aae9b7d6c2cce5e249c62abb246d237 100644 (file)
@@ -47,60 +47,6 @@ class PreconditionIdentity
                                      */
     template<class VECTOR>
     void Tvmult (VECTOR&, const VECTOR&) const;
-                                    /**
-                                     * Apply preconditioner.
-                                     */
-    template<class VECTOR>
-    void operator () (VECTOR&, const VECTOR&) const;
-};
-
-
-/**
- * Jacobi preconditioner using matrix built-in function.  The MATRIX
- * class used is required to have a function
- * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double}
- *
- * @author Guido Kanschat, 2000
- */
-template <class MATRIX = SparseMatrix<double> >
-class PreconditionJacobi
-{
-  public:
-                                    /**
-                                     * Initialize matrix and
-                                     * relaxation parameter. The
-                                     * matrix is just stored in the
-                                     * preconditioner object. The
-                                     * relaxation parameter should be
-                                     * larger than zero and smaller
-                                     * than 2 for numerical
-                                     * reasons. It defaults to 1.
-                                     */
-    void initialize (const MATRIX& A, const double omega = 1.);
-                                    /**
-                                     * Apply preconditioner.
-                                     */
-    template<class VECTOR>
-    void vmult (VECTOR&, const VECTOR&) const;
-                                    /**
-                                     * Apply transpose
-                                     * preconditioner. Since this is
-                                     * a symmetric preconditioner,
-                                     * this function is the same as
-                                     * @ref{vmult}.
-                                     */
-    template<class VECTOR>
-    void Tvmult (VECTOR&, const VECTOR&) const;
-
-  private:
-                                    /**
-                                     * Pointer to the matrix object.
-                                     */
-    SmartPointer<const MATRIX> A;
-                                    /**
-                                     * Relaxation parameter.
-                                     */
-    double omega;
 };
 
 
@@ -174,8 +120,8 @@ class PreconditionUseMatrix
                                      * of this object with the two
                                      * arguments given here.
                                      */
-    void operator() (VECTOR       &dst,
-                    const VECTOR &src) const;
+    void vmult (VECTOR       &dst,
+               const VECTOR &src) const;
 
   private:
                                     /**
@@ -192,94 +138,120 @@ class PreconditionUseMatrix
 
 
 /**
- * Preconditioner for builtin relaxation methods.
- * Application of this preconditioner includes
- * use of the #precondition_...# methods of #SparseMatrix#.
+ * Base class for other preconditioners.
+ * Here, only some common features Jacobi, SOR and SSOR preconditioners
+ * are implemented. For preconditioning, refer to derived classes.
  *
- * \subsection{Use}
- * You will usually not want to create a named object of this type,
- * although possible. The most common use is like this:
- * \begin{verbatim}
- *    SolverGMRES<SparseMatrix<double>,
- *                Vector<double> >      gmres(control,memory,500);
- *
- *    gmres.solve (matrix, solution, right_hand_side,
- *                PreconditionRelaxation<SparseMatrix<double>,Vector<double> >
- *                (matrix,&SparseMatrix<double>::template precondition_Jacobi,
- *                   0.5));
- * \end{verbatim}
- * This creates an unnamed object to be passed as the fourth parameter to
- * the solver function of the #SolverGMRES# class. It assumes that the
- * #SparseMatrix# class has a function #precondition_Jacobi# taking two
- * vectors (source and destination) and a relaxation value as parameters. (Unlike
- * for the #PreconditionUseMatrix# class, this time it should work, with
- * relaxation parameter $0.5$.)
- *
- * Note that due to the default template parameters, the above example
- * could be written shorter as follows:
- * \begin{verbatim}
- *    ...
- *    gmres.solve (matrix, solution, right_hand_side,
- *                PreconditionRelaxation<>
- *                  (matrix,
- *                    &SparseMatrix<double>::template precondition_Jacobi,
- *                    0.5));
- * \end{verbatim}
- *
- * @author Guido Kanschat, Wolfgang Bangerth, 1999
+ * @author Guido Kanschat, 2000
  */
-template<class MATRIX = SparseMatrix<double>, class VECTOR = Vector<double> >
+template<class MATRIX = SparseMatrix<double> >
 class PreconditionRelaxation
 {
   public:
                                     /**
-                                     * Type of the preconditioning
-                                     * function of the matrix.
+                                     * Initialize matrix and
+                                     * relaxation parameter. The
+                                     * matrix is just stored in the
+                                     * preconditioner object. The
+                                     * relaxation parameter should be
+                                     * larger than zero and smaller
+                                     * than 2 for numerical
+                                     * reasons. It defaults to 1.
                                      */
-    typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&,
-                                          typename MATRIX::value_type) const;
-    
+    void initialize (const MATRIX& A, const double omega = 1.);
+  protected:
                                     /**
-                                     * 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.
+                                     * Pointer to the matrix object.
                                      */
-    PreconditionRelaxation(const MATRIX      &M,
-                          const function_ptr method,
-                          const double       omega = 1.);
-    
+    SmartPointer<const MATRIX> A;
                                     /**
-                                     * Execute preconditioning. Calls the
-                                     * function passed to the constructor
-                                     * of this object with the two
-                                     * arguments given here, and the
-                                     * relaxation parameter passed to the
-                                     * constructor.
+                                     * Relaxation parameter.
                                      */
-    void operator() (VECTOR&, const VECTOR&) const;
+    double omega;
+};
 
-  private:
+/**
+ * Jacobi preconditioner using matrix built-in function.  The MATRIX
+ * class used is required to have a function
+ * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double}
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionJacobi : public PreconditionRelaxation<MATRIX>
+{
+  public:
                                     /**
-                                     * Pointer to the matrix in use.
+                                     * Apply preconditioner.
                                      */
-    const MATRIX& matrix;
-    
+    template<class VECTOR>
+    void vmult (VECTOR&, const VECTOR&) const;
                                     /**
-                                     * Pointer to the preconditioning
-                                     * function.
+                                     * Apply transpose
+                                     * preconditioner. Since this is
+                                     * a symmetric preconditioner,
+                                     * this function is the same as
+                                     * @ref{vmult}.
                                      */
-    const function_ptr precondition;
-    
+    template<class VECTOR>
+    void Tvmult (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * SOR preconditioner using matrix built-in function.  The MATRIX
+ * class used is required to have functions
+ * @p{precondition_SOR(VECTOR&, const VECTOR&, double)} and
+ * @p{precondition_TSOR(VECTOR&, const VECTOR&, double)}.
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionSOR : public PreconditionRelaxation<MATRIX>
+{
+  public:
                                     /**
-                                     * Relaxation parameter.
+                                     * Apply preconditioner.
                                      */
-    double omega;
+    template<class VECTOR>
+    void vmult (VECTOR&, const VECTOR&) const;
+                                    /**
+                                     * Apply transpose
+                                     * preconditioner.
+                                     */
+    template<class VECTOR>
+    void Tvmult (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * SSOR preconditioner using matrix built-in function.  The MATRIX
+ * class used is required to have a function
+ * @p{precondition_SSOR(VECTOR&, const VECTOR&, double}
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
+{
+  public:
+                                    /**
+                                     * Apply preconditioner.
+                                     */
+    template<class VECTOR>
+    void vmult (VECTOR&, const VECTOR&) const;
+                                    /**
+                                     * Apply transpose
+                                     * preconditioner. Since this is
+                                     * a symmetric preconditioner,
+                                     * this function is the same as
+                                     * @ref{vmult}.
+                                     */
+    template<class VECTOR>
+    void Tvmult (VECTOR&, const VECTOR&) const;
 };
 
+
 /**
  * Preconditioner using an iterative solver.  This preconditioner uses
  * a fully initialized LAC iterative solver for the approximate
@@ -307,7 +279,7 @@ class PreconditionLACSolver
                                      * Execute preconditioning.
                                      */
     template<class VECTOR>
-    void operator() (VECTOR&, const VECTOR&) const;
+    void vmult (VECTOR&, const VECTOR&) const;
 
   private:
                                     /**
@@ -351,10 +323,17 @@ class PreconditionedMatrix
                          VectorMemory<VECTOR>&  mem);
 
                                     /**
-                                     * Preconditioned matrix-vector-product.
+                                     * Preconditioned
+                                     * matrix-vector-product.
                                      */
     void vmult (VECTOR& dst, const VECTOR& src) const;
 
+                                    /**
+                                     * Transposed preconditioned
+                                     * matrix-vector-product.
+                                     */
+    void Tvmult (VECTOR& dst, const VECTOR& src) const;
+
                                     /**
                                      * Residual $b-PAx$.
                                      */
@@ -380,13 +359,6 @@ class PreconditionedMatrix
 
 /* ---------------------------------- Inline functions ------------------- */
 
-template<class VECTOR>
-inline void
-PreconditionIdentity::operator () (VECTOR& dst, const VECTOR& src) const
-{
-  dst = src;
-}
-
 template<class VECTOR>
 inline void
 PreconditionIdentity::vmult (VECTOR& dst, const VECTOR& src) const
@@ -405,56 +377,95 @@ PreconditionIdentity::Tvmult (VECTOR& dst, const VECTOR& src) const
 
 template <class MATRIX>
 inline void
-PreconditionJacobi<MATRIX>::initialize (const MATRIX& rA, double o)
+PreconditionRelaxation<MATRIX>::initialize (const MATRIX& rA, double o)
 {
   A = &rA;
   omega = o;
 }
 
+//----------------------------------------------------------------------//
+
 template <class MATRIX>
 template<class VECTOR>
 inline void
 PreconditionJacobi<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
 {
-//TODO: Assert object not initialized
+  Assert (A!=0, ExcNotInitialized());
+  A->precondition_Jacobi (dst, src, omega);
+}
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionJacobi<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (A!=0, ExcNotInitialized());
   A->precondition_Jacobi (dst, src, omega);
 }
 
 
 //----------------------------------------------------------------------//
 
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (A!=0, ExcNotInitialized());
+  A->precondition_SOR (dst, src, omega);
+}
 
-template<class MATRIX, class VECTOR>
-PreconditionUseMatrix<MATRIX,VECTOR>::PreconditionUseMatrix(const MATRIX& M,
-                                                            function_ptr method)
-               :
-               matrix(M), precondition(method)
-{}
 
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (A!=0, ExcNotInitialized());
+  A->precondition_TSOR (dst, src, omega);
+}
 
-template<class MATRIX, class VECTOR>
-void
-PreconditionUseMatrix<MATRIX,VECTOR>::operator() (VECTOR& dst,
-                                                  const VECTOR& src) const
+
+//----------------------------------------------------------------------//
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
 {
-  (matrix.*precondition)(dst, src);
+  Assert (A!=0, ExcNotInitialized());
+  A->precondition_SSOR (dst, src, omega);
+}
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (A!=0, ExcNotInitialized());
+  A->precondition_SSOR (dst, src, omega);
 }
 
+
+//----------------------------------------------------------------------//
+
+
 template<class MATRIX, class VECTOR>
-PreconditionRelaxation<MATRIX,VECTOR>::PreconditionRelaxation(const MATRIX& M,
-                                                            function_ptr method,
-                                                              double omega)
+PreconditionUseMatrix<MATRIX,VECTOR>::PreconditionUseMatrix(const MATRIX& M,
+                                                            function_ptr method)
                :
-               matrix(M), precondition(method), omega(omega)
+               matrix(M), precondition(method)
 {}
 
 
 template<class MATRIX, class VECTOR>
 void
-PreconditionRelaxation<MATRIX,VECTOR>::operator() (VECTOR& dst,
-                                                  const VECTOR& src) const
+PreconditionUseMatrix<MATRIX,VECTOR>::vmult (VECTOR& dst,
+                                            const VECTOR& src) const
 {
-  (matrix.*precondition)(dst, src, omega);
+  (matrix.*precondition)(dst, src);
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -471,8 +482,8 @@ PreconditionLACSolver<SOLVER,MATRIX,PRECONDITION>
 template<class SOLVER, class MATRIX, class PRECONDITION>
 template<class VECTOR>
 void
-PreconditionLACSolver<SOLVER,MATRIX,PRECONDITION>::operator() (VECTOR& dst,
-                                                                const VECTOR& src) const
+PreconditionLACSolver<SOLVER,MATRIX,PRECONDITION>::vmult (VECTOR& dst,
+                                                         const VECTOR& src) const
 {
   solver.solve(matrix, dst, src, precondition);
 }
@@ -499,10 +510,25 @@ PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
   VECTOR* h = mem.alloc();
   h->reinit(src);
   A.vmult(*h, src);
-  P(dst, *h);
+  P.vmult(dst, *h);
   mem.free(h);
 }
 
+
+template<class MATRIX, class PRECOND, class VECTOR>
+inline void
+PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
+::Tvmult (VECTOR& dst,
+        const VECTOR& src) const
+{
+  VECTOR* h = mem.alloc();
+  h->reinit(src);
+  A.Tvmult(*h, src);
+  P.Tvmult(dst, *h);
+  mem.free(h);
+}
+
+
 template<class MATRIX, class PRECOND, class VECTOR>
 inline double
 PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
@@ -513,7 +539,7 @@ PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
   VECTOR* h = mem.alloc();
   h->reinit(src);
   A.vmult(*h, src);
-  P(dst, *h);
+  P.vmult(dst, *h);
   mem.free(h);
   dst.sadd(-1.,1.,rhs);
   return dst.l2_norm ();
index bc45f9e90782084883fc976b5df29e0e549b09eb..af9d1be944d4f7fb3b5e39505480d1253f57c5ef 100644 (file)
@@ -41,7 +41,7 @@ template <typename number> class Vector;
  * For all matrices that are empty above and below the diagonal
  * blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner
  * is a direct solver. For all matrices that are empty only above the diagonal blocks
- * (e.g. the matrices one gets by the DG method with downstream numbering) the
+ * (e.g. the matrices one gets by the DG method with downstream numbering)
  * #BlockSOR# is a direct solver.
  * 
  * This first implementation of the #PreconditionBlock# assumes the
@@ -84,6 +84,18 @@ class PreconditionBlock: public Subscriptor
                                      */
     virtual ~PreconditionBlock();
 
+                                    /**
+                                     * Initialize matrix and block
+                                     * size.  We store the matrix and
+                                     * the block size in the
+                                     * preconditioner object. In a
+                                     * second step, the inverses of
+                                     * the diagonal blocks may be
+                                     * computed.
+                                     */
+    void initialize (const SparseMatrix<number>& A,
+                    const unsigned int block_size);
+
                                     /**
                                      * Deletes the inverse diagonal block
                                      * matrices if existent, sets the
@@ -95,27 +107,20 @@ class PreconditionBlock: public Subscriptor
     virtual void clear();
 
                                     /**
-                                     * Takes the matrix that should be used
-                                     * for the preconditioning. A reference
-                                     * to it is stored within this class,
-                                     * but ownership of the matrix remains
-                                     * with the caller of this function.
+                                     * Use only the inverse of the
+                                     * first diagonal block to save
+                                     * memory and computation time.
+                                     *
+                                     * Possible applications:
+                                     * computing on a cartesian grid,
+                                     * all diagonal blocks are the
+                                     * same or all diagonal blocks
+                                     * are at least similar and
+                                     * inversion of one of them still
+                                     * yields a preconditioner.
                                      */
-    void use_matrix(const SparseMatrix<number> &M);
+    void set_same_diagonal ();
     
-                                    /**
-                                     * Set the right block size before calling
-                                     * #invert_diagblocks# or calling the
-                                     * #operator ()# function
-                                     * of #PreconditionBlockSOR#
-                                     * or #PreconditionBlockJacobi#.
-                                     * If #block_size==1# BlockSOR or 
-                                     * BlockJacobi are equal to the
-                                     * standard SOR or Jacobi 
-                                     * preconditioner.
-                                     */
-    void set_block_size (const unsigned int bsize);
-
                                     /**
                                      * Stores the inverse of
                                      * the diagonal blocks
@@ -147,6 +152,12 @@ class PreconditionBlock: public Subscriptor
                                      */
     unsigned int block_size () const;
 
+                                    /**
+                                     * Determine, whether inverses
+                                     * have been computed.
+                                     */
+    bool inverses_ready () const;
+    
                                     /**
                                      * Exception
                                      */
@@ -170,29 +181,37 @@ class PreconditionBlock: public Subscriptor
                                      */
     DeclException0 (ExcInverseMatricesAlreadyExist);
 
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcBlockSizeNotSet);
-
                                     /**
                                      * Exception
                                      */
     DeclException0 (ExcMatrixNotSquare);
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNoMatrixGivenToUse);
    
   protected:
+                                    /**
+                                     * Access to the inverse diagonal blocks.
+                                     * 
+                                     */
+    const FullMatrix<inverse_type>& inverse (unsigned int i) const;
+    
                                     /**
                                      * Size of the blocks. Each diagonal
                                      * block is assumed to be of the
                                      * same size.
                                      */
     unsigned int blocksize;
-    
+
+                                    /**
+                                     * Pointer to the matrix. Make sure that
+                                     * the matrix exists as long as this class
+                                     * needs it, i.e. until calling #invert_diagblocks#,
+                                     * or (if the inverse matrices should not be
+                                     * stored) until the last call of the 
+                                     * preconditoining #vmult# function of the
+                                     * derived classes.
+                                     */
+    SmartPointer<const SparseMatrix<number> > A;
+
+  private:
                                     /**
                                      * Storage of the inverse matrices of
                                      * the diagonal blocks matrices as
@@ -201,18 +220,13 @@ class PreconditionBlock: public Subscriptor
                                      * using #inverse_type=float# saves memory
                                      * in comparison with #inverse_type=double#.
                                      */
-    vector<FullMatrix<inverse_type> > inverse;
+    vector<FullMatrix<inverse_type> > _inverse;
 
                                     /**
-                                     * Pointer to the matrix. Make sure that
-                                     * the matrix exists as long as this class
-                                     * needs it, i.e. until calling #invert_diagblocks#,
-                                     * or (if the inverse matrices should not be
-                                     * stored) until the last call of the 
-                                     * preconditoining #operator()# function of the
-                                     * derived classes.
+                                     * Flag for diagonal compression.
+                                     * @see set_same_diagonal()
                                      */
-    SmartPointer<const SparseMatrix<number> > A;
+    bool same_diagonal;
 };
 
 
@@ -226,11 +240,6 @@ class PreconditionBlockJacobi : public PreconditionBlock<number,inverse_type>
   public:
                                     /**
                                      * Execute block Jacobi preconditioning.
-                                     * Make sure that the right
-                                     * block size
-                                     * of the matrix is set by
-                                     * #set_block_size#
-                                     * before calling this function.
                                      *
                                      * This function will automatically use the
                                      * inverse matrices if they exist, if not
@@ -239,7 +248,13 @@ class PreconditionBlockJacobi : public PreconditionBlock<number,inverse_type>
                                      * matrices in each preconditioning step.
                                      */
     template <typename number2>
-    void operator() (Vector<number2>&, const Vector<number2>&) const;
+    void vmult (Vector<number2>&, const Vector<number2>&) const;
+
+                                    /**
+                                     * Same as #vmult#, since Jacobi is symmetric.
+                                     */
+    template <typename number2>
+    void Tvmult (Vector<number2>&, const Vector<number2>&) const;
 };
 
 
@@ -273,11 +288,6 @@ class PreconditionBlockSOR : public PreconditionBlock<number,inverse_type>
     
                                     /**
                                      * Execute block SOR preconditioning.
-                                     * Make sure that the right
-                                     * block size
-                                     * of the matrix is set by
-                                     * #set_block_size#
-                                     * before calling this function.
                                      *
                                      * This function will automatically use the
                                      * inverse matrices if they exist, if not
@@ -290,7 +300,7 @@ class PreconditionBlockSOR : public PreconditionBlock<number,inverse_type>
                                      * BlockSOR is a direct solver.
                                      */
     template <typename number2>
-    void operator() (Vector<number2>&, const Vector<number2>&) const;
+    void vmult (Vector<number2>&, const Vector<number2>&) const;
 
   private:
                                     /**
index ccf17906352117a377bc3f0043e6aa73492b20a0..04ac0056629d9bac29da768c6e3d8d2ee8c54d1b 100644 (file)
@@ -14,6 +14,7 @@
 #define __deal2__precondition_block_templates_h
 
 
+#include <base/exceptions.h>
 #include <lac/precondition_block.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 template <typename number, typename inverse_type>
 PreconditionBlock<number,inverse_type>::PreconditionBlock ():
                blocksize(0),
-               A(0)  {};
+               A(0),
+               same_diagonal(false)
+{};
 
 
 template <typename number, typename inverse_type>
 PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
 {
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
+  if (_inverse.size()!=0)
+    _inverse.erase(_inverse.begin(), _inverse.end());
 }
 
 
 template <typename number, typename inverse_type>
 void PreconditionBlock<number,inverse_type>::clear ()
 {
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
-  blocksize=0;
+  if (_inverse.size()!=0)
+    _inverse.erase(_inverse.begin(), _inverse.end());
+  blocksize     = 0;
+  same_diagonal = false;  
 }
 
 
 template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::use_matrix(
-  const SparseMatrix<number> &M)
+void PreconditionBlock<number,inverse_type>
+::initialize (const SparseMatrix<number> &M, unsigned int bsize)
 {
+  Assert (M.m() == M.n(), ExcMatrixNotSquare());
   A = &M;
+  Assert (bsize>0, ExcIndexRange(bsize, 1, M.m()));
+  Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
+  blocksize=bsize;
 }
 
 
 template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::set_block_size(unsigned int bsize) {
-  blocksize=bsize;
+const FullMatrix<inverse_type>&
+PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
+{
+  if (same_diagonal)
+    return _inverse[0];
+  
+  Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size()));
+  return _inverse[i];
 }
 
 
 template <typename number, typename inverse_type>
-unsigned int PreconditionBlock<number,inverse_type>::block_size() const {
+unsigned int PreconditionBlock<number,inverse_type>::block_size() const
+{
   return blocksize;
 }
 
 
+template <typename number, typename inverse_type>
+void
+PreconditionBlock<number,inverse_type>::set_same_diagonal()
+{
+  Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+  same_diagonal = true;
+}
+
+
+template <typename number, typename inverse_type>
+bool
+PreconditionBlock<number,inverse_type>::inverses_ready() const
+{
+  return (_inverse.size() != 0);
+}
+
+
 template <typename number, typename inverse_type>
 void PreconditionBlock<number,inverse_type>::invert_diagblocks()
 {
-  Assert (A!=0, ExcNoMatrixGivenToUse());
-  const SparseMatrix<number> &M=*A;
-  Assert (M.m() == M.n(), ExcMatrixNotSquare());
-  Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist());
+  Assert (A!=0, ExcNotInitialized());
+  Assert (blocksize!=0, ExcNotInitialized());
 
-  Assert (blocksize!=0, ExcBlockSizeNotSet());
-  Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
+  const SparseMatrix<number> &M=*A;
+  Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist());
 
   const unsigned int n_cells = M.m()/blocksize;
 
-                                  // cell_row, cell_column are the
-                                  // numbering of the blocks (cells).
-                                  // row_cell, column_cell are the local
-                                  // numbering of the unknowns in the
-                                  // blocks.
-                                  // row, column are the global numbering
-                                  // of the unkowns.
-
-                                  // set the #inverse# array to the right
-                                  // size. we could do it like this:
-                                  // inverse = vector<>(n_cells,FullMatrix<>())
-                                  // but this would involve copying many
-                                  // FullMatrix objects.
-                                  //
-                                  // the following is a neat trick which
-                                  // avoids copying
-  if (true)
-    {
-      vector<FullMatrix<inverse_type> > tmp(n_cells,
-                                           FullMatrix<inverse_type>(blocksize));
-      inverse.swap (tmp);
-    };
-  
   FullMatrix<inverse_type> M_cell(blocksize);
-  
-  for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
+
+  if (same_diagonal)
     {
-      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-       for (unsigned int column_cell=0, column=cell*blocksize;
-            column_cell<blocksize; ++column_cell, ++column)
-         M_cell(row_cell,column_cell)=M.el(row,column);
-//      try
-//     {
+                                      // Invert only the first block
+                                      // This is a copy of the code in the
+                                      // 'else' part, stripped of the outer loop
+      _inverse.resize(1);
+
+      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+       for (unsigned int column_cell=0; column_cell<blocksize; ++column_cell)
+         M_cell(row_cell,column_cell)=M.el(row_cell,column_cell);
+
       if (blocksize <=4)
        {  
-         inverse[cell].invert(M_cell);
+         _inverse[0].invert(M_cell);
        }
       else
        {
          M_cell.gauss_jordan();
-         inverse[cell]=M_cell;
+         _inverse[0]=M_cell;
        }
+    }
+  else
+    {
+                                      // cell_row, cell_column are the
+                                      // numbering of the blocks (cells).
+                                      // row_cell, column_cell are the local
+                                      // numbering of the unknowns in the
+                                      // blocks.
+                                      // row, column are the global numbering
+                                      // of the unkowns.
+
+                                      // set the #_inverse# array to the right
+                                      // size. we could do it like this:
+                                      // _inverse = vector<>(n_cells,FullMatrix<>())
+                                      // but this would involve copying many
+                                      // FullMatrix objects.
+                                      //
+                                      // the following is a neat trick which
+                                      // avoids copying
+      vector<FullMatrix<inverse_type> > tmp(n_cells,
+                                           FullMatrix<inverse_type>(blocksize));
+      _inverse.swap (tmp);
+
+      M_cell.clear ();
+      
+      for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
+       {
+         for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+           for (unsigned int column_cell=0, column=cell*blocksize;
+                column_cell<blocksize; ++column_cell, ++column)
+             M_cell(row_cell,column_cell)=M.el(row,column);
+//      try
+//     {
+         if (blocksize <=4)
+           {  
+             _inverse[cell].invert(M_cell);
+           }
+         else
+           {
+             M_cell.gauss_jordan();
+             _inverse[cell]=M_cell;
+           }
 //      }
 //      catch (ExcNotImplemented &)
+       }
     }
 }
 
@@ -130,8 +183,8 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
 template <typename number, typename inverse_type>
 template <typename number2>
 void PreconditionBlockJacobi<number,inverse_type>
-::operator() (Vector<number2>       &dst,
-             const Vector<number2> &src) const
+::vmult (Vector<number2>       &dst,
+        const Vector<number2> &src) const
 {
                                   // introduce the following typedef
                                   // since in the use of exceptions,
@@ -150,20 +203,11 @@ void PreconditionBlockJacobi<number,inverse_type>
                                   // exceptions that do not take
                                   // args...
   typedef PreconditionBlock<number,inverse_type> BaseClass;
-  Assert(A!=0, typename BaseClass::ExcNoMatrixGivenToUse());
+  Assert(A!=0, ExcNotInitialized());
   
   const SparseMatrix<number> &M=*A;
   const unsigned int n_cells=M.m()/blocksize;
 
-  Assert (M.m() == M.n(),
-         typename BaseClass::ExcMatrixNotSquare());
-  Assert (blocksize!=0,
-         typename BaseClass::ExcBlockSizeNotSet());
-  Assert (M.m()%blocksize==0,
-         typename BaseClass::ExcWrongBlockSize(blocksize, M.m()));
-  Assert (inverse.size()==0 || inverse.size()==n_cells,
-         typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells));
-
   Vector<number2> b_cell(blocksize), x_cell(blocksize);
 
                                       // cell_row, cell_column are the
@@ -175,7 +219,7 @@ void PreconditionBlockJacobi<number,inverse_type>
                                       // of the unkowns.
   unsigned int row, row_cell, begin_diag_block=0;
 
-  if (inverse.size()==0)
+  if (!inverses_ready())
     {
       FullMatrix<number> M_cell(blocksize);
       for (unsigned int cell=0; cell<n_cells; ++cell)
@@ -203,7 +247,7 @@ void PreconditionBlockJacobi<number,inverse_type>
          {
            b_cell(row_cell)=src(row);
          }
-       inverse[cell].vmult(x_cell, b_cell);
+       inverse(cell).vmult(x_cell, b_cell);
                                         // distribute x_cell to dst
        for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
          dst(row)=x_cell(row_cell);
@@ -211,15 +255,29 @@ void PreconditionBlockJacobi<number,inverse_type>
        begin_diag_block+=blocksize;
       }
 }
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::Tvmult (Vector<number2>       &dst,
+        const Vector<number2> &src) const
+{
+  vmult(dst, src);
+}
+
+
 /*--------------------- PreconditionBlockSOR -----------------------*/
 
 template<typename number, typename inverse_type>
 PreconditionBlockSOR<number,inverse_type>::PreconditionBlockSOR(const number omega):
-               omega(omega)  {}
+               omega(omega)
+{}
 
 
 template<typename number, typename inverse_type>
-PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR(){}
+PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR()
+{}
 
 
 
@@ -231,12 +289,12 @@ PreconditionBlockSOR<number,inverse_type>::set_omega(number om)
 }
 
 
-
+//TODO: implement Tvmult
 
 template <typename number, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2>       &dst,
-                                                           const Vector<number2> &src) const
+void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &dst,
+                                                      const Vector<number2> &src) const
 {
                                   // introduce the following typedef
                                   // since in the use of exceptions,
@@ -256,20 +314,11 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2>
                                   // args...
   typedef PreconditionBlock<number,inverse_type> BaseClass;
 
-  Assert (A!=0, typename BaseClass::ExcNoMatrixGivenToUse());
+  Assert (A!=0, ExcNotInitialized());
   
   const SparseMatrix<number> &M=*A;
   const unsigned int n_cells=M.m()/blocksize;
 
-  Assert (M.m() == M.n(),
-         typename BaseClass::ExcMatrixNotSquare());
-  Assert (blocksize!=0,
-         typename BaseClass::ExcBlockSizeNotSet());
-  Assert (M.m()%blocksize==0,
-         typename BaseClass::ExcWrongBlockSize(blocksize, M.m()));
-  Assert (inverse.size()==0 || inverse.size()==n_cells,
-         typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells));
-
   const SparsityPattern &spars    = M.get_sparsity_pattern();
   const unsigned int    *rowstart = spars.get_rowstart_indices();
   const unsigned int    *columns  = spars.get_column_numbers();
@@ -286,7 +335,7 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2>
   unsigned int row, column, row_cell, begin_diag_block=0;
   number2 b_cell_row;
 
-  if (inverse.size()==0)
+  if (!inverses_ready())
     {
       FullMatrix<number> M_cell(blocksize);
       for (unsigned int cell=0; cell<n_cells; ++cell)
@@ -324,7 +373,7 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2>
                }
            b_cell(row_cell)=b_cell_row;
          }
-       inverse[cell].vmult(x_cell, b_cell);
+       inverse(cell).vmult(x_cell, b_cell);
                                         // distribute x_cell to dst
        for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
          dst(row)=omega*x_cell(row_cell);
index 182b07c136d80c010b53b52f14a0009e1e50b141..dd23b001f6ab416c50ae92039e204e77f96a3826 100644 (file)
@@ -63,12 +63,12 @@ class SolverBicgstab : private Solver<VECTOR>
                                     /**
                                      * Solve primal problem only.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
-          const Preconditioner& precondition);
+          const PRECONDITIONER& precondition);
 
   protected:
                                     /**
@@ -257,7 +257,7 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
       beta   = rhobar * alpha / (rho * omega);
       rho    = rhobar;
       p.sadd(beta, 1., r, -beta*omega, v);
-      precondition(y,p);
+      precondition.vmult(y,p);
       A.vmult(v,y);
       rhobar = rbar * v;
 
@@ -269,7 +269,7 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
        return typename Solver<VECTOR>::ReturnState(breakdown);
     
       s.equ(1., r, -alpha, v);
-      precondition(z,s);
+      precondition.vmult(z,s);
       A.vmult(t,z);
       rhobar = t*s;
       omega = rhobar/(t*t);
index 1b7f41e5e67e7b86059a486e49b0265c65b61796..0b2d344ce1d410e73e62f58665840b970cd7f682 100644 (file)
@@ -70,12 +70,12 @@ class SolverCG : private Solver<VECTOR>
                                     /**
                                      * Solver method.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
-          const Preconditioner& precondition);
+          const PRECONDITIONER& precondition);
 
   protected:
                                     /**
@@ -164,12 +164,12 @@ SolverCG<VECTOR>::print_vectors(const unsigned int,
 
 
 template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
 SolverCG<VECTOR>::solve (const MATRIX &A,
                         VECTOR       &x,
                         const VECTOR &b,
-                        const Preconditioner& precondition)
+                        const PRECONDITIONER& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -210,7 +210,7 @@ SolverCG<VECTOR>::solve (const MATRIX &A,
     };
   
   g.scale(-1.);
-  precondition(h,g);
+  precondition.vmult(h,g);
  
   d.equ(-1.,h);
  
@@ -234,7 +234,7 @@ SolverCG<VECTOR>::solve (const MATRIX &A,
       if (conv)
        break;
       
-      precondition(h,g);
+      precondition.vmult(h,g);
       
       beta = gh;
       gh   = g*h;
index 2dbf439fe807be7ebcb5f5f73b3563dc933bb42a..d264fdb8ad9f1f4956b9659771ee258dd73f488a 100644 (file)
@@ -271,7 +271,7 @@ unsigned int dim = 0;
       if (left_precondition)
        {
          A.residual(p,x,b);
-         precondition(v,p);
+         precondition.vmult(v,p);
        } else {
          A.residual(v,x,b);
        };
@@ -311,9 +311,9 @@ unsigned int dim = 0;
          if (left_precondition)
            {
              A.vmult(p, *tmp_vectors[inner_iteration]);
-             precondition(vv,p);
+             precondition.vmult(vv,p);
            } else {
-             precondition(p,*tmp_vectors[inner_iteration]);
+             precondition.vmult(p,*tmp_vectors[inner_iteration]);
              A.vmult(vv,p);
            };
          
@@ -376,7 +376,7 @@ unsigned int dim = 0;
          p = 0.;
          for (unsigned int i=0; i<dim; ++i)
            p.add(h(i), *tmp_vectors[i]);
-         precondition(v,p);
+         precondition.vmult(v,p);
          x.add(1.,v);
        };
 
index bff486943c8dd96cbcf9b8ae8c1c5c52f9bee7c5..100c05740c9db6e55e0e49c99ba53250e3929b7d 100644 (file)
@@ -74,12 +74,12 @@ class SolverMinRes : private Solver<VECTOR>
                                     /**
                                      * Solver method.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState
     solve (const MATRIX &A,
           VECTOR       &x,
           const VECTOR &b,
-          const Preconditioner& precondition);
+          const PRECONDITIONER& precondition);
 
                                     /**
                                      * Exception
@@ -168,12 +168,12 @@ SolverMinRes<VECTOR>::print_vectors(const unsigned int,
 
 
 template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
 SolverMinRes<VECTOR>::solve (const MATRIX &A,
                             VECTOR       &x,
                             const VECTOR &b,
-                            const Preconditioner& precondition)
+                            const PRECONDITIONER& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -234,7 +234,7 @@ SolverMinRes<VECTOR>::solve (const MATRIX &A,
                                   // positiv definite and symmetric
 
                                   // M v = u[1]
-  precondition (v,u[1]);
+  precondition.vmult (v,u[1]);
   
   delta[1] = v * u[1];
                                   // Preconditioner positive
@@ -270,7 +270,7 @@ SolverMinRes<VECTOR>::solve (const MATRIX &A,
                                       // precondition: solve M v = u[2]
                                       // Preconditioner has to be positiv
                                       // definite and symmetric.
-      precondition(v,u[2]);
+      precondition.vmult(v,u[2]);
  
       delta[2] = v * u[2];
 
index 20ce79596fc78123fd398f9f138c941ba463a44b..31946803a8218a6176432dfe1d61cb4c79e59fd0 100644 (file)
@@ -102,12 +102,12 @@ class SolverQMRS : private Solver<VECTOR>
                                     /**
                                      * Solver method.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState
     solve (const MATRIX &A,
                       VECTOR       &x,
                       const VECTOR &b,
-                      const Preconditioner& precondition);
+                      const PRECONDITIONER& precondition);
 
                                     /**
                                      * Interface for derived class.
@@ -168,9 +168,9 @@ class SolverQMRS : private Solver<VECTOR>
                                     /**
                                      * The iteration loop itself.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState 
-    iterate(const MATRIX& A, const Preconditioner& precondition);
+    iterate(const MATRIX& A, const PRECONDITIONER& precondition);
                                     /**
                                      * The current iteration step.
                                      */
@@ -210,12 +210,12 @@ SolverQMRS<VECTOR>::print_vectors(const unsigned int,
 
 
 template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
 SolverQMRS<VECTOR>::solve (const MATRIX &A,
                           VECTOR       &x,
                           const VECTOR &b,
-                          const Preconditioner& precondition)
+                          const PRECONDITIONER& precondition)
 {
   deallog.push("QMRS");
   
@@ -266,10 +266,10 @@ SolverQMRS<VECTOR>::solve (const MATRIX &A,
 
 
 template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
 SolverQMRS<VECTOR>::iterate(const MATRIX& A,
-                           const Preconditioner& precondition)
+                           const PRECONDITIONER& precondition)
 {
 /* Remark: the matrix A in the article is the preconditioned matrix.
  * Therefore, we have to precondition x before we compute the first residual.
@@ -296,7 +296,7 @@ SolverQMRS<VECTOR>::iterate(const MATRIX& A,
   d.reinit(x);
   
                                   // Apply right preconditioning to x
-  precondition(q,x);  
+  precondition.vmult(q,x);  
                                   // Preconditioned residual
   res = A.residual(v, q, b);
 
@@ -305,7 +305,7 @@ SolverQMRS<VECTOR>::iterate(const MATRIX& A,
 
   p = v;
   
-  precondition(q,p);
+  precondition.vmult(q,p);
  
   tau = v.norm_sqr();
   //deallog << "tau:" << tau << endl;
@@ -359,12 +359,12 @@ while (state == SolverControl::iterate)
        return breakdown;
                                       // Step 7
       rho_old = rho;
-      precondition(q,v);
+      precondition.vmult(q,v);
       rho = q*v;
 
       beta = rho/rho_old;
       p.sadd(beta,v);
-      precondition(q,p);
+      precondition.vmult(q,p);
     }
   return exceeded;
 }
index eb5a6332e955a86fe87804d4823aef9b68185499..fc6ad64a472cb77ebe34dc449fbbf3dafc19cd58 100644 (file)
@@ -73,11 +73,11 @@ class SolverRichardson : private Solver<VECTOR>
                                     /**
                                      * Solve $Ax=b$ for $x$.
                                      */
-    template<class MATRIX, class Preconditioner>
+    template<class MATRIX, class PRECONDITIONER>
     typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
                                                       VECTOR       &x,
                                                       const VECTOR &b,
-                                                      const Preconditioner& precondition);
+                                                      const PRECONDITIONER& precondition);
 
                                     /**
                                      * Set the damping-coefficient.
@@ -146,12 +146,12 @@ SolverRichardson<VECTOR>::SolverRichardson(SolverControl &cn,
 
 
 template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
 SolverRichardson<VECTOR>::solve (const MATRIX &A,
                                 VECTOR       &x,
                                 const VECTOR &b,
-                                const Preconditioner& precondition)
+                                const PRECONDITIONER& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -170,7 +170,7 @@ SolverRichardson<VECTOR>::solve (const MATRIX &A,
       if (conv != SolverControl::iterate)
        break;
 
-      precondition(d,r);
+      precondition.vmult(d,r);
       x.add(additional_data.omega,d);
       print_vectors(iter,x,r,d);
     }
index e092baa40a953eaa4eb23f5c2a100d420b377b62..2dae4fe8237258844bb6147483887aa67dd46051 100644 (file)
@@ -192,8 +192,8 @@ class SparseILU : protected SparseMatrix<number>
                                      * Same as #apply_decomposition#, format for LAC.
                                      */
     template <typename somenumber>
-    void operator() (Vector<somenumber>       &dst,
-                    const Vector<somenumber> &src) const;
+    void vmult (Vector<somenumber>       &dst,
+               const Vector<somenumber> &src) const;
 
                                     /**
                                      * Exception
@@ -218,8 +218,8 @@ class SparseILU : protected SparseMatrix<number>
 template <typename number>
 template <typename somenumber>
 void
-SparseILU<number>::operator() (Vector<somenumber>       &dst,
-                              const Vector<somenumber> &src) const
+SparseILU<number>::vmult (Vector<somenumber>       &dst,
+                         const Vector<somenumber> &src) const
 {
   apply_decomposition(dst, src);
 }
index cb318d9b9468029c53dd581d9220ebf01b4f5f89..bffe5898ac123960ae0ea9253618c2c4345c90ed 100644 (file)
@@ -203,8 +203,8 @@ class SparseVanka
                                      * update vector in #dst#.
                                      */
     template<typename number2>
-    void operator() (Vector<number2>       &dst,
-                    const Vector<number2> &src) const;
+    void vmult (Vector<number2>       &dst,
+               const Vector<number2> &src) const;
 
                                     /**
                                      * Exception
@@ -256,7 +256,7 @@ class SparseVanka
                                      * function with a
                                      * #vector<bool>(n_dofs,true)#.
                                      *
-                                     * The #operator()# of this class
+                                     * The #vmult# of this class
                                      * of course calls this function
                                      * with a null pointer
                                      */
@@ -516,7 +516,7 @@ class SparseBlockVanka : public SparseVanka<number>
                                      * Apply the preconditioner.
                                      */
     template<typename number2>
-    void operator() (Vector<number2>       &dst,
+    void vmult (Vector<number2>       &dst,
                     const Vector<number2> &src) const;
     
   private:
index bb77852398769cca80d6884957158328e9299e1c..b3d0df796d8587c3dc42c81b0ccd11d415aed935 100644 (file)
@@ -186,13 +186,13 @@ SparseVanka<number>::compute_inverse (const unsigned int    row,
   inverses[row]->gauss_jordan();
 };
 
-
+//TODO: implement Tvmult
 
 template<typename number>
 template<typename number2>
 void
-SparseVanka<number>::operator ()(Vector<number2>       &dst,
-                                const Vector<number2> &src) const
+SparseVanka<number>::vmult (Vector<number2>       &dst,
+                           const Vector<number2> &src) const
 {
                                   // first set output vector to zero
   dst.clear ();
@@ -543,11 +543,12 @@ SparseBlockVanka<number>::compute_dof_masks (const SparseMatrix<number> &M,
     };
 };
 
+//TODO: implement Tvmult
 
 template <typename number>
 template <typename number2>
-void SparseBlockVanka<number>::operator() (Vector<number2>       &dst,
-                                          const Vector<number2> &src) const
+void SparseBlockVanka<number>::vmult (Vector<number2>       &dst,
+                                     const Vector<number2> &src) const
 {
   dst.clear ();
 
index d94b57ac28e1902988db9fa7aee5423a9d5ab627..426c93ce768235700abc3b836b4b3b7e8625bd65 100644 (file)
@@ -35,25 +35,36 @@ template class PreconditionBlock<double, double>;
 // explicit instantiations for "float" PreconditionBlock
 template class PreconditionBlockJacobi<float, float>;
 
-template void PreconditionBlockJacobi<float, float>::operator() (
-  Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::operator() (
-  Vector<double> &, const Vector<double> &) const;
-
+template void PreconditionBlockJacobi<float, float>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
 
 template class PreconditionBlockJacobi<double, float>;
 
-template void PreconditionBlockJacobi<double, float>::operator() (
-  Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::operator() (
-  Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
 
 template class PreconditionBlockJacobi<double, double>;
 
-template void PreconditionBlockJacobi<double, double>::operator() (
-  Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::operator() (
-  Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
 
 /*--------------------- PreconditionBlockGaussSeidel -----------------------*/
 
@@ -61,9 +72,9 @@ template void PreconditionBlockJacobi<double, double>::operator() (
 // explicit instantiations for "float" PreconditionBlock
 template class PreconditionBlockSOR<float, float>;
 
-template void PreconditionBlockSOR<float, float>::operator() (
+template void PreconditionBlockSOR<float, float>::vmult (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::operator() (
+template void PreconditionBlockSOR<float, float>::vmult (
   Vector<double> &, const Vector<double> &) const;
 
 
@@ -76,16 +87,16 @@ template void PreconditionBlockSOR<float, float>::operator() (
 template class PreconditionBlockSOR<double, float>;
 
 
-template void PreconditionBlockSOR<double, float>::operator() (
+template void PreconditionBlockSOR<double, float>::vmult (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::operator() (
+template void PreconditionBlockSOR<double, float>::vmult (
   Vector<double> &, const Vector<double> &) const;
 
 template class PreconditionBlockSOR<double, double>;
 
-template void PreconditionBlockSOR<double, double>::operator() (
+template void PreconditionBlockSOR<double, double>::vmult (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::operator() (
+template void PreconditionBlockSOR<double, double>::vmult (
   Vector<double> &, const Vector<double> &) const;
 
 
index b23b455b4dece7ecb68f309c4905858a7ca9ca76..9828eea8ddb3717e28014411fef01dee1d46df08 100644 (file)
 template class SparseVanka<float>;
 template class SparseVanka<double>;
 
-template void SparseVanka<double>::operator () (Vector<float>       &dst,
-                                               const Vector<float> &src) const;
-template void SparseVanka<double>::operator () (Vector<double>       &dst,
-                                               const Vector<double> &src) const;
+template void SparseVanka<double>::vmult (Vector<float>       &dst,
+                                         const Vector<float> &src) const;
+template void SparseVanka<double>::vmult (Vector<double>       &dst,
+                                         const Vector<double> &src) const;
 
 
 template class SparseBlockVanka<float>;
 template class SparseBlockVanka<double>;
 
-template void SparseBlockVanka<double>::operator () (Vector<float>       &dst,
-                                                    const Vector<float> &src) const;
-template void SparseBlockVanka<double>::operator () (Vector<double>       &dst,
-                                                    const Vector<double> &src) const;
+template void SparseBlockVanka<double>::vmult (Vector<float>       &dst,
+                                              const Vector<float> &src) const;
+template void SparseBlockVanka<double>::vmult (Vector<double>       &dst,
+                                              const Vector<double> &src) const;
index 891b8cbe5829a1bf85da40f70376fabcc510bbd0..d28e210a3e5ee3fdaa02e32ab759237bd5d92e1f 100644 (file)
@@ -154,8 +154,8 @@ MGSmootherLAC::smooth (const unsigned int level,
   SolverControl control(2,1.e-300,false,false);
   PrimitiveVectorMemory<> mem;
   SolverRichardson<> rich(control, mem);
-  PreconditionRelaxation<>
-    prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+  PreconditionSSOR<> prec;
+  prec.initialize((*matrices)[level], 1.);
 
   rich.solve((*matrices)[level], u, rhs, prec);
 }
@@ -377,10 +377,8 @@ void LaplaceProblem<dim>::solve ()
       PrimitiveVectorMemory<> vector_memory;
       SolverCG<>              cg (solver_control, vector_memory);
 
-      PreconditionRelaxation<>
-       preconditioner(global_system_matrix,
-                      &SparseMatrix<double>::template precondition_SSOR<double>,
-                      1.2);
+      PreconditionSSOR<> preconditioner;
+      preconditioner.initialize(global_system_matrix, 1.2);
       
       solution.clear ();
       cg.solve (global_system_matrix, solution, system_rhs,
index 7eb3c1c726bfee027f4b3d5ab0155776b6eac8d1..d3e75ea57204c77ae9932a2828115bc26b93c75b 100644 (file)
@@ -266,8 +266,8 @@ void LaplaceProblem<Vector,Matrix,Sparsity>::solve ()
   PrimitiveVectorMemory<Vector> vector_memory;
   SolverCG<Vector>        cg (solver_control, vector_memory);
 
-  PreconditionRelaxation<Matrix,Vector> preconditioner
-    (system_matrix, &Matrix::precondition_Jacobi, 0.8);
+  PreconditionJacobi<Matrix> preconditioner;
+  preconditioner.initialize (system_matrix, 0.8);
   
   cg.solve (system_matrix, solution, system_rhs,
            preconditioner);
index 84cfb6f38e7faf18faec73807509f61c183e9bfc..52241020bf3f6e9710d35ecd7aac98c330142280 100644 (file)
@@ -180,8 +180,8 @@ MGSmootherLAC::smooth (const unsigned int level,
   SolverControl control(2,1.e-300,false,false);
   PrimitiveVectorMemory<> mem;
   SolverRichardson<> rich(control, mem);
-  PreconditionRelaxation<>
-    prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+  PreconditionSSOR<> prec;
+  prec.initialize((*matrices)[level], 1.);
 
   rich.solve((*matrices)[level], u, rhs, prec);
 }

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.