From a33cd296be3c0932f4f33dfbf1434546e915404e Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 14 Aug 2001 09:39:53 +0000 Subject: [PATCH] Check the requirements on vectors and matrices in the linear solvers. Minor clean-ups. git-svn-id: https://svn.dealii.org/trunk@4880 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver.h | 96 +++++++++++---------- deal.II/lac/include/lac/solver_bicgstab.h | 3 + deal.II/lac/include/lac/solver_cg.h | 22 +++-- deal.II/lac/include/lac/solver_gmres.h | 23 +++-- deal.II/lac/include/lac/solver_minres.h | 35 ++++---- deal.II/lac/include/lac/solver_qmrs.h | 38 ++++---- deal.II/lac/include/lac/solver_richardson.h | 40 +++++---- 7 files changed, 138 insertions(+), 119 deletions(-) diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h index f638719d1d..1a8b33a7ee 100644 --- a/deal.II/lac/include/lac/solver.h +++ b/deal.II/lac/include/lac/solver.h @@ -19,84 +19,86 @@ class SolverControl; /** - * Base class for iterative solvers. - * -//TODO:[?] * HAS TO BE UPDATED! - * - * This class defines possible + * Base class for iterative solvers. This class defines possible * return states of linear solvers and provides interfaces to a memory * pool and the control object. * * * @sect3{Requirements for template classes} * - * The class is templated to allow for different matrix and vector - * classes, since iterative solvers do not rely on any special structure - * of matrices or the format of storage. However, there are some common - * requirements a matrix or vector type must fulfil to qualify as an - * applicable type for the solvers in this hierarchy. These requirements - * are listed following. The classes do not declare any concrete - * class, they are rather intended to form a `signature' which a concrete - * class has to conform to. + * Since iterative solvers do not rely on any special structure of + * matrices or the format of storage, but only require that matrices + * and vector define certain operations such as matrix-vector + * products, or scalar products between vectors, this class as well as + * the derived classes implementing concrete linear solvers are + * templated on the types of matrices and vectors. However, there are + * some common requirements a matrix or vector type must fulfill to + * qualify as an applicable type for the solvers in this + * hierarchy. These requirements are listed following. The listed + * classes are not any concrete class, they are rather intended to + * form a `signature' which a concrete class has to conform to. Note + * that the matrix and vector classes within this library of course + * conform to this interface. * * @begin{verbatim} * class Matrix * { * public: - * // Application to a Vector - * void vmult (Vector& dst, const Vector& src) const; - * - * // Application of a preconditioner to - * // a Vector, i.e. $dst=\tilde A^(-1) src$, - * // where $\tilde A^(-1)$ is an approximation - * // to the inverse if the matrix stored in - * // this object. - * void precondition (Vector& dst, const Vector& src) const; + * // Application of matrix to vector src. + * // write result into dst + * void vmult (Vector &dst, const Vector &src) const; * * // Application of transpose to a Vector. - * // Only used by special iterative methods. - * void T_vmult (Vector& dst, const Vector& src) const; - * - * // Application of a transposed preconditioner - * // to a Vector. Only used by special - * // iterative methods - * - * void T_precondition (Vector& dst, const Vector& src) const; + * // Only used by certain iterative methods. + * void Tvmult (Vector &dst, const Vector &src) const; * }; * * * class Vector * { * public: + * // resize and/or clear vector. note + * // that the second argument must have + * // a default value equal to false + * void reinit (const unsigned int size, + * bool leave_elements_uninitialized = false); + * * // scalar product - * double operator * (const Vector& v) const; + * double operator * (const Vector &v) const; * * // addition of vectors * // $y = y + x$. - * void add (const Vector& x); + * void add (const Vector &x); + * * // $y = y + ax$. - * void add (double a, const Vector& x); + * void add (const double a, + * const Vector &x); * - * // scaled addition of vectors - * // $y = ay + x$. - * void sadd (double a, - * const Vector& x); * // $y = ay + bx$. - * void sadd (double a, - * double b, const Vector& x); - * // $y = ay + bx + cz$. - * void sadd (double a, - * double b, const Vector& x, - * double c, const Vector& z); + * void sadd (const double a, + * const double b, + * const Vector &x); * * // $y = ax$. - * void equ (double a, const Vector& x); - * // $y = ax + bz$. - * void equ (double a, const Vector& x, - * double b, const Vector& z); + * void equ (const double a, + * const Vector &x); + * + * // scale the elements of the vector + * // by a fixed value + * void scale (const double a); + * + * // return the l2 norm of the vector + * double l2_norm () const; * }; * @end{verbatim} * + * In addition, for some solvers there has to be a global function + * @p{swap(vector &a, vector &b)} that exchanges the values of the two vectors. + * + * The preconditioners used must have the same interface as matrices, + * i.e. in particular they have to provide a member function @p{vmult} + * which denotes the application of the preconditioner. + * * * @sect3{AdditionalData} * diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h index 9bff8e3e00..534720c708 100644 --- a/deal.II/lac/include/lac/solver_bicgstab.h +++ b/deal.II/lac/include/lac/solver_bicgstab.h @@ -22,6 +22,9 @@ /** * Bicgstab algorithm by van der Vorst. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * Like all other solver classes, this class has a local structure called * @p{AdditionalData} which is used to pass additional parameters to the * solver, like damping parameters or the number of temporary vectors. We diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index efacba01da..98c50a0f05 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -25,6 +25,9 @@ /** * Preconditioned cg method. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * Like all other solver classes, this class has a local structure called * @p{AdditionalData} which is used to pass additional parameters to the * solver, like damping parameters or the number of temporary vectors. We @@ -85,14 +88,15 @@ class SolverCG : public Subscriptor, private Solver virtual ~SolverCG (); /** - * Solver method. + * Solve the linear system $Ax=b$ + * for x. */ template void - solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); protected: /** @@ -189,10 +193,10 @@ SolverCG::print_vectors(const unsigned int, template template void -SolverCG::solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition) +SolverCG::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { SolverControl::State conv=SolverControl::iterate; diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index b5a4103754..ef6f3f6de2 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -69,6 +69,9 @@ * of temporary vectors as commented upon above. By default, the number * of these vectors is set to 30. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * @author Wolfgang Bangerth */ template > @@ -107,13 +110,15 @@ class SolverGMRES : public Subscriptor, private Solver const AdditionalData &data=AdditionalData()); /** - * Solver method. + * Solve the linear system $Ax=b$ + * for x. */ template - void solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + void + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); DeclException1 (ExcTooFewTmpVectors, int, @@ -205,10 +210,10 @@ SolverGMRES::givens_rotation (Vector &h, template template void -SolverGMRES::solve (const MATRIX& A, - VECTOR & x, - const VECTOR& b, - const PRECONDITIONER& precondition) +SolverGMRES::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { // this code was written a very // long time ago by people not diff --git a/deal.II/lac/include/lac/solver_minres.h b/deal.II/lac/include/lac/solver_minres.h index 1103da70d6..711c7ba1d4 100644 --- a/deal.II/lac/include/lac/solver_minres.h +++ b/deal.II/lac/include/lac/solver_minres.h @@ -24,6 +24,9 @@ /** * Preconditioned MinRes method. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * Like all other solver classes, this class has a local structure called * @p{AdditionalData} which is used to pass additional parameters to the * solver, like damping parameters or the number of temporary vectors. We @@ -74,14 +77,15 @@ class SolverMinRes : public Subscriptor, private Solver virtual ~SolverMinRes (); /** - * Solver method. + * Solve the linear system $Ax=b$ + * for x. */ template void - solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); /** * Exception @@ -172,10 +176,10 @@ SolverMinRes::print_vectors(const unsigned int, template template typename Solver::ReturnState -SolverMinRes::solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition) +SolverMinRes::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { SolverControl::State conv=SolverControl::iterate; @@ -246,21 +250,20 @@ SolverMinRes::solve (const MATRIX &A, r_l2 = r0; - u[0].reinit(VS,0); + u[0].reinit(VS); delta[0] = 1.; - m[0].reinit(VS,0); - m[1].reinit(VS,0); - m[2].reinit(VS,0); + m[0].reinit(VS); + m[1].reinit(VS); + m[2].reinit(VS); conv = control().check(0,r_l2); while (conv==SolverControl::iterate) - { - + { if (delta[1]!=0) v.scale(1./sqrt(delta[1])); else - v.reinit(VS,0); + v.reinit(VS); A.vmult(u[2],v); u[2].add (-sqrt(delta[1]/delta[0]), u[0]); diff --git a/deal.II/lac/include/lac/solver_qmrs.h b/deal.II/lac/include/lac/solver_qmrs.h index b4c8788594..2bca2aec6c 100644 --- a/deal.II/lac/include/lac/solver_qmrs.h +++ b/deal.II/lac/include/lac/solver_qmrs.h @@ -32,6 +32,9 @@ * preconditioner is used: left preconditioning seems to require the * inverse. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * Like all other solver classes, this class has a local structure called * @p{AdditionalData} which is used to pass additional parameters to the * solver, like damping parameters or the number of temporary vectors. We @@ -106,14 +109,15 @@ class SolverQMRS : public Subscriptor, private Solver const AdditionalData &data=AdditionalData()); /** - * Solver method. + * Solve the linear system $Ax=b$ + * for x. */ template void - solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); /** * Interface for derived class. @@ -218,10 +222,10 @@ SolverQMRS::print_vectors(const unsigned int, template template void -SolverQMRS::solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition) +SolverQMRS::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { deallog.push("QMRS"); @@ -277,8 +281,8 @@ SolverQMRS::solve (const MATRIX &A, template template bool -SolverQMRS::iterate(const MATRIX& A, - const PRECONDITIONER& precondition) +SolverQMRS::iterate(const MATRIX &A, + 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. @@ -319,12 +323,9 @@ SolverQMRS::iterate(const MATRIX& A, precondition.vmult(q,p); tau = v.norm_sqr(); - //deallog << "tau:" << tau << std::endl; rho = q*v; - //deallog << "rho:" << rho << std::endl; - - -while (state == SolverControl::iterate) + + while (state == SolverControl::iterate) { step++; it++; // Step 1 @@ -337,7 +338,6 @@ while (state == SolverControl::iterate) return true; // Step 3 alpha = rho/sigma; - //deallog << "alpha:" << alpha << std::endl; v.add(-alpha,t); // Step 4 @@ -346,10 +346,6 @@ while (state == SolverControl::iterate) psi = 1./(1.+theta); tau *= theta*psi; - //deallog << "psi:" << psi << std::endl; - //deallog << "theta:" << theta << std::endl; - //deallog << "tau:" << tau << std::endl; - d.sadd(psi*theta_old, psi*alpha, p); x.add(d); diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h index 71b87816d6..60682892c0 100644 --- a/deal.II/lac/include/lac/solver_richardson.h +++ b/deal.II/lac/include/lac/solver_richardson.h @@ -22,6 +22,9 @@ * Implementation of the richardson iteration method. The stopping criterion * is the norm of the residual. * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the @ref{Solver} base class. + * * Like all other solver classes, this class has a local structure called * @p{AdditionalData} which is used to pass additional parameters to the * solver, like damping parameters or the number of temporary vectors. We @@ -76,22 +79,25 @@ class SolverRichardson : public Subscriptor, private Solver virtual ~SolverRichardson (); /** - * Solve $Ax=b$ for $x$. + * Solve the linear system $Ax=b$ + * for x. */ template - void solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + void + solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); /** * Solve $A^Tx=b$ for $x$. */ template - void Tsolve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition); + void + Tsolve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition); /** * Set the damping-coefficient. @@ -169,10 +175,10 @@ SolverRichardson::~SolverRichardson() template template void -SolverRichardson::solve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition) +SolverRichardson::solve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { SolverControl::State conv=SolverControl::iterate; @@ -218,10 +224,10 @@ SolverRichardson::solve (const MATRIX &A, template template void -SolverRichardson::Tsolve (const MATRIX &A, - VECTOR &x, - const VECTOR &b, - const PRECONDITIONER& precondition) +SolverRichardson::Tsolve (const MATRIX &A, + VECTOR &x, + const VECTOR &b, + const PRECONDITIONER &precondition) { SolverControl::State conv=SolverControl::iterate; -- 2.39.5