add relaxation methods to preconditioner module

git-svn-id: https://svn.dealii.org/trunk@22826 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/doxygen/headers/preconditioners.h b/deal.II/doc/doxygen/headers/preconditioners.h
index 737f1c0e929df8b51eb7eb5bfc099e3d21802fda..0c7c7ccae0fa8c229070cc5c49c36c192cd900a9 100644 (file)
--- a/deal.II/doc/doxygen/headers/preconditioners.h
+++ b/deal.II/doc/doxygen/headers/preconditioners.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -13,7 +13,9 @@
 
 
 /**
- * @defgroup Preconditioners Preconditioners
+ * @defgroup Preconditioners Preconditioners and Relaxation Operators
+ *
+ * <h3>Preconditioners</h3>
  *
  * Preconditioners are used to accelerate the iterative solution of linear
  * systems. Typical preconditioners are Jacobi, Gauss-Seidel, or SSOR, but the
@@ -21,19 +23,72 @@
  * decompositions (ILU). In addition, sparse direct solvers can be used as
  * preconditioners when available.
  *
- * When talking of preconditioners, we usually expect them to be used
- * in Krylov-space methods. In that case, the act as operators: given
- * a vector $x$, produce the result $y=P^{-1}x$ of the multiplication
- * with the preconditioning operator $P^{-1}$.
+ * Broadly speaking, preconditioners are operators, which are
+ * multiplied with a matrix to improve conditioning. The idea is, that
+ * the preconditioned system <i>P<sup>-1</sup>Ax = P<sup>-1</sup>b</i>
+ * is much easier to solve than the original system <i>Ax = b</i>.What
+ * this means exactly depends on the structure of  the matrix and
+ * cannot be discussed here in generality. For symmetric, positive
+ * definite matrices <i>A</i> and <i>P</i>, it means that the spectral
+ * condition number (the quotient of greatest and smallest eigenvalue)
+ * of <i>P<sup>-1</sup>A</i> is much smaller than the one of <i>A</i>.
  *
- * However, some preconditioners can also be used
- * in the standard linear defect correction iteration,
+ * At hand of the simplest example, Richardson iteration, implemented
+ * in SolverRichardson, the preconditioned iteration looks like
  * @f[
- *  x^{k+1} = x^k - P^{-1} \bigl(A x^k - b\bigr),
+ *  x^{k+1} = x^k - P^{-1} \bigl(A x^k - b\bigr).
  * @f]
- * where <i>P<sup>-1</sup></i> is again the preconditioner. Thus,
- * preconditioning amounts to applying a linear operator to the
- * residual.
+ * Accordingly, preconditioning amounts to applying a linear operator to the
+ * residual, and consequently, the action of the preconditioner
+ * <i>P<sup>-1</sup></i> is implemented as <tt>vmult()</tt>. The
+ * generic interface is like for matrices
+ * @code
+ * class PRECONDITIONER
+ * {
+ *   template <class VECTOR>
+ *   void vmult(VECTOR& dst, const VECTOR& src) const;
+ *
+ *   template <class VECTOR>
+ *   void Tvmult(VECTOR& dst, const VECTOR& src) const;
+ * }
+ * @endcode
+ * It is implemented in all the preconditioner classes in this module.
+ * 
+ * When used
+ * in Krylov space methods, it is up to the method, whether it simply
+ * replaces multiplications with <i>A</i> by those with
+ * <i>P<sup>-1</sup>A</i> (for instance SolverBicgstab), or does more
+ * sophisticated things. SolverCG for instance uses
+ * <i>P<sup>-1</sup></i> to define an inner product, which is the
+ * reason why it requires a symmetric, positive definite operator <i>P</i>.
+ *
+ * <h3>Relaxation methods</h3>
+ *
+ * Many preconditioners rely on an additive splitting <i>A = P - N</i>
+ * into two matrices. In this case, the iteration step of the
+ * Richardson method above can be simplified to
+ * @f[
+ *  x^{k+1} = P^{-1} \bigl(N x^k + b\bigr),
+ * @f]
+ * thus avoiding multiplication with <i>A</i> completely. We call
+ * operators mapping the previous iterate <i>x<sup>k</sup></i> to the
+ * next iterate in this way relaxation operators. Their generic
+ * interface is
+ * @code
+ * class RELAXATION
+ * {
+ *   template <class VECTOR>
+ *   void step(VECTOR& newstep, const VECTOR& oldstep, const VECTOR& rhs) const;
+ *
+ *   template <class VECTOR>
+ *   void Tstep(VECTOR& newstep, const VECTOR& oldstep, const VECTOR& rhs) const;
+ * }
+ * @endcode
+ * The classes with names starting with <tt>Relaxation</tt> in this
+ *   module implement this interface, as well as the preconditioners
+ *   PreconditionJacobi, PreconditionSOR, PreconditionSSORP,
+ *   reconditionBlockJacobi, PreconditionBlockSOR, and
+ *   PreconditionBlockSSOR.
  *
  * <h3>The interface</h3>
  *

diff --git a/deal.II/include/deal.II/lac/precondition_block_base.h b/deal.II/include/deal.II/lac/precondition_block_base.h
index b3c3a6a004d78b515a1ccaf23bd49abeef8deeaa..1777a923226466e1133ae76c6bd333e57c1c2de7 100644 (file)
--- a/deal.II/include/deal.II/lac/precondition_block_base.h
+++ b/deal.II/include/deal.II/lac/precondition_block_base.h
@@ -244,6 +244,15 @@ class PreconditionBlockBase
                                      */
     DeclException0 (ExcDiagonalsNotStored);
 
+                                    /**
+                                     * You are accessing a diagonal
+                                     * block, assuming that it has a certain
+                                     * type. But, the method used for
+                                     * inverting the diagonal blocks
+                                     * does not use this type
+                                     */
+    DeclException0 (ExcInverseNotAvailable);
+    
   protected:
                                     /**
                                      * The method used for inverting blocks.
@@ -565,6 +574,8 @@ inline
 FullMatrix<number>&
 PreconditionBlockBase<number>::inverse(unsigned int i)
 {
+  Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
+  
   if (same_diagonal())
     return var_inverse_full[0];
   
@@ -578,6 +589,8 @@ inline
 Householder<number>&
 PreconditionBlockBase<number>::inverse_householder(unsigned int i)
 {
+  Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
+  
   if (same_diagonal())
     return var_inverse_householder[0];
   
@@ -591,6 +604,8 @@ inline
 LAPACKFullMatrix<number>&
 PreconditionBlockBase<number>::inverse_svd(unsigned int i)
 {
+  Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
+  
   if (same_diagonal())
     return var_inverse_svd[0];
   

diff --git a/deal.II/include/deal.II/lac/relaxation_block.h b/deal.II/include/deal.II/lac/relaxation_block.h
index d56119559329be3532986e943b7d9f0043693100..ca937fdd7d4ac9007a8e032f18d3be168f9f90ad 100644 (file)
--- a/deal.II/include/deal.II/lac/relaxation_block.h
+++ b/deal.II/include/deal.II/lac/relaxation_block.h
@@ -41,6 +41,7 @@ DEAL_II_NAMESPACE_OPEN
  * preconditioner, since its implementation relies on a straight
  * forward implementation of the Gauss-Seidel process.
  *
+ * @ingroup Preconditioners
  * @author Guido Kanschat
  * @date 2010
  */
@@ -240,6 +241,7 @@ class RelaxationBlock :
  * overlapping. On the other hand, this class does not implement the
  * preconditioner interface expected by Solver objects.
  *
+ * @ingroup Preconditioners
  * @author Guido Kanschat
  * @date 2010
  */
@@ -307,13 +309,6 @@ class RelaxationBlockSOR : public virtual Subscriptor,
                                      */
     template <typename number2>
     void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
-
-  protected:
-                                    /**
-                                     * Constructor to be used by
-                                     * RelaxationBlockSSOR.
-                                     */
-//    RelaxationBlockSOR(bool store);
 };
 
 
@@ -329,6 +324,7 @@ class RelaxationBlockSOR : public virtual Subscriptor,
  * not implement the preconditioner interface expected by Solver
  * objects.
  *
+ * @ingroup Preconditioners
  * @author Guido Kanschat
  * @date 2010
  */