From: kanschat Date: Sat, 20 Nov 2010 00:00:06 +0000 (+0000) Subject: add relaxation methods to preconditioner module X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ef000df1b4e32d6fe1baf951fec554120807772;p=dealii-svn.git 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 737f1c0e92..0c7c7ccae0 100644 --- 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 + * + *

Preconditioners

* * 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 P-1Ax = P-1b + * is much easier to solve than the original system Ax = b.What + * this means exactly depends on the structure of the matrix and + * cannot be discussed here in generality. For symmetric, positive + * definite matrices A and P, it means that the spectral + * condition number (the quotient of greatest and smallest eigenvalue) + * of P-1A is much smaller than the one of A. * - * 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 P-1 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 + * P-1 is implemented as vmult(). The + * generic interface is like for matrices + * @code + * class PRECONDITIONER + * { + * template + * void vmult(VECTOR& dst, const VECTOR& src) const; + * + * template + * 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 A by those with + * P-1A (for instance SolverBicgstab), or does more + * sophisticated things. SolverCG for instance uses + * P-1 to define an inner product, which is the + * reason why it requires a symmetric, positive definite operator P. + * + *

Relaxation methods

+ * + * Many preconditioners rely on an additive splitting A = P - N + * 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 A completely. We call + * operators mapping the previous iterate xk to the + * next iterate in this way relaxation operators. Their generic + * interface is + * @code + * class RELAXATION + * { + * template + * void step(VECTOR& newstep, const VECTOR& oldstep, const VECTOR& rhs) const; + * + * template + * void Tstep(VECTOR& newstep, const VECTOR& oldstep, const VECTOR& rhs) const; + * } + * @endcode + * The classes with names starting with Relaxation in this + * module implement this interface, as well as the preconditioners + * PreconditionJacobi, PreconditionSOR, PreconditionSSORP, + * reconditionBlockJacobi, PreconditionBlockSOR, and + * PreconditionBlockSSOR. * *

The interface

* 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 b3c3a6a004..1777a92322 100644 --- 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& PreconditionBlockBase::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& PreconditionBlockBase::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& PreconditionBlockBase::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 d561195593..ca937fdd7d 100644 --- 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 void Tstep (Vector& dst, const Vector& 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 */