From: guido Date: Mon, 26 Apr 1999 19:43:05 +0000 (+0000) Subject: precondition now a class of it's own X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b751fc1d698954cce32e0dc7546d84e2a417382b;p=dealii-svn.git precondition now a class of it's own git-svn-id: https://svn.dealii.org/trunk@1207 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h new file mode 100644 index 0000000000..afc432ef27 --- /dev/null +++ b/deal.II/lac/include/lac/precondition.h @@ -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 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 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 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 +void +PreconditionIdentity::operator() (VECTOR& dst, const VECTOR& src) const +{ + dst = src; +} + +template +PreconditionUseMatrix::PreconditionUseMatrix(const MATRIX& M, + function_ptr method) + : + matrix(M), precondition(method) +{} + + +template +void +PreconditionUseMatrix::operator() (VECTOR& dst, + const VECTOR& src) const +{ + (matrix.*precondition)(dst, src); +} + +template +PreconditionRelaxation::PreconditionRelaxation(const MATRIX& M, + function_ptr method, + double omega) + : + matrix(M), precondition(method), omega(omega) +{} + + +template +void +PreconditionRelaxation::operator() (VECTOR& dst, + const VECTOR& src) const +{ + (matrix.*precondition)(dst, src, omega); +} + +#endif diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h index 25fa6d4a73..a57b1567a7 100644 --- a/deal.II/lac/include/lac/solver.h +++ b/deal.II/lac/include/lac/solver.h @@ -13,7 +13,11 @@ template 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 &); - - /** - * 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 + ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition) const; /** * Access to control object. diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index ae7cae0a1f..1b6b5e1da6 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -25,9 +25,11 @@ class SolverBicgstab // /** * Solve primal problem only. */ - virtual ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b); + template + ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition); protected: /** @@ -121,7 +123,8 @@ class SolverBicgstab // /** * The iteration loop itself. */ - ReturnState iterate(); + template + ReturnState iterate(const Preconditioner& precondition); }; @@ -135,7 +138,8 @@ SolverBicgstab::SolverBicgstab(SolverControl &cn, {} -template double +template +double SolverBicgstab::criterion() { res = MA->residual(*Vt, *Vx, *Vb); @@ -143,7 +147,8 @@ SolverBicgstab::criterion() } -template < class Matrix, class Vector > SolverControl::State +template < class Matrix, class Vector > +SolverControl::State SolverBicgstab::start() { res = MA->residual(*Vr, *Vx, *Vb); @@ -154,8 +159,10 @@ SolverBicgstab::start() return state; } -template Solver::ReturnState -SolverBicgstab::iterate() +template +template +Solver::ReturnState +SolverBicgstab::iterate(const Preconditioner& precondition) { SolverControl::State state = SolverControl::iterate; alpha = omega = rho = 1.; @@ -175,7 +182,7 @@ SolverBicgstab::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::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::iterate() } -template Solver::ReturnState +template +template +Solver::ReturnState SolverBicgstab::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::solve(const Matrix &A, { deallog << "Go!" << endl; if (start() == SolverControl::success) break; - state = iterate(); + state = iterate(precondition); } while (state == breakdown); diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 5cd18d9641..30d2ff731c 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -17,7 +17,8 @@ * @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and documentation by Wolfgang Bangerth */ template -class SolverPCG : public Solver { +class SolverPCG : public Solver +{ public: /** @@ -29,9 +30,11 @@ class SolverPCG : public Solver { /** * Solver method. */ - virtual ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b); + template + ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition); protected: /** @@ -76,10 +79,13 @@ double SolverPCG::criterion() template +template Solver::ReturnState SolverPCG::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::solve (const Matrix &A, }; g.scale(-1.); - A.precondition(h,g); + precondition(h,g); d.equ(-1.,h); @@ -137,7 +143,7 @@ SolverPCG::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 index 6b3760898b..0000000000 --- a/deal.II/lac/include/lac/solver_gauss_seidel.h +++ /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 -#include -#include - -/** - * 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 SolverGaussSeidel : public Solver { - public: - /** - * Constructor. Takes a SolverControl and some - VectorMemory. VectorMemory is not used. - */ - SolverGaussSeidel (SolverControl &cn, VectorMemory &mem) : - Solver (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 -SolverGaussSeidel::ReturnState -SolverGaussSeidel::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 -inline double -SolverGaussSeidel::criterion() -{ - return sqrt(res2); -} - -/*---------------------------- solver_gauss_seidel.h ---------------------------*/ -/* end of #ifndef __solver_gauss_seidel_H */ -#endif -/*---------------------------- solver_gauss_seidel.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 5db3b45de1..76f138924e 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -45,9 +45,11 @@ * 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 SolverGMRES : public Solver { +class SolverGMRES : public Solver +{ public: /** * Constructor. @@ -59,9 +61,11 @@ class SolverGMRES : public Solver { /** * Solver method. */ - virtual ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b); + template + ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition); DeclException1 (ExcTooFewTmpVectors, int, @@ -137,10 +141,12 @@ SolverGMRES::givens_rotation (Vector &h, template +template Solver::ReturnState SolverGMRES::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::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::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::solve (const Matrix& A, p = 0.; for (unsigned int i=0; i { /** * Solve $Ax=b$ for $x$. */ - virtual ReturnState solve (const Matrix &A, - Vector &x, - const Vector &b); + template + 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 { template +template Solver::ReturnState SolverRichardson::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::solve (const Matrix &A, if (conv != SolverControl::iterate) break; - A.precondition(d,r); + precondition(d,r); x.add(omega,d); }