From ec6bf4675217bb6217559da85d8848a9aacd006f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 3 Aug 1998 11:29:35 +0000 Subject: [PATCH] Move the iterative solver classes and helpers from the indepenant MIA library to the lac library. The reason is missing support of MIA and a structure which hardly seems designed but rather hacked together... git-svn-id: https://svn.dealii.org/trunk@469 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver.h | 235 +++++++++++++++++++++++ deal.II/lac/include/lac/solver_cg.h | 160 +++++++++++++++ deal.II/lac/include/lac/solver_control.h | 189 ++++++++++++++++++ deal.II/lac/include/lac/vector_memory.h | 80 ++++++++ deal.II/lac/source/solver_control.cc | 86 +++++++++ 5 files changed, 750 insertions(+) create mode 100644 deal.II/lac/include/lac/solver.h create mode 100644 deal.II/lac/include/lac/solver_cg.h create mode 100644 deal.II/lac/include/lac/solver_control.h create mode 100644 deal.II/lac/include/lac/vector_memory.h create mode 100644 deal.II/lac/source/solver_control.cc diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h new file mode 100644 index 0000000000..687ad83cd3 --- /dev/null +++ b/deal.II/lac/include/lac/solver.h @@ -0,0 +1,235 @@ +/*---------------------------- solver.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_H +#define __solver_H +/*---------------------------- solver.h ---------------------------*/ + + + +// forward declaration +class SolverControl; +template class VectorMemory; + + + +/** + * 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. + * + * 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. + * \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 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; + * }; + * + * + * class Vector + * { + * public: + * // scalar product + * double operator * (const Vector& v) const; + * + * // addition of vectors + * // $y = y + x$. + * void add (const Vector& x); + * // $y = y + ax$. + * void add (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); + * + * // $y = ax$. + * void equ (double a, const Vector& x); + * // $y = ax + bz$. + * void equ (double a, const Vector& x, + * double b, const Vector& z); + * }; + * \end{verbatim} + */ +template +class Solver +{ + public: + /** + * Declare possible return values of a + * solver object. + */ + enum ReturnState { + success=0, exceeded, breakdown + }; + + + + /** + * Constructor. Assign a control + * object which stores the required + * precision and an object to provide + * memory. + * + * Of both objects, a reference is + * stored, so it is the user's + * responsibility to guarantee that the + * lifetime of the two arguments is at + * least as long as that of the solver + * object. + */ + Solver (SolverControl &, VectorMemory &); + + /** + * Destructor. This one is virtual, + * overload it in derived classes if + * necessary. + */ + virtual ~Solver() {} + + /** + * Solver procedure. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b) = 0; + + /** + * Access to control object. + */ + SolverControl& control() const; + + protected: + /** + * Calculation of convergence + * criterion. To allow further + * flexibility in solvers, the + * convergence criterion can be + * implemented by the user. Each + * method has it's standard + * criterion (usually some kind of + * residual) already implemented. + */ + virtual double criterion() = 0; + + /** + * Additional entry point for + * examination. Allows access to + * internal variables of a solver + * at a given point dependent of + * the actual solver method. Does + * nothing by default and should + * only be implemented to test + * iterative methods. + */ + virtual void expert() {}; + + /** + * Control structure. + */ + SolverControl &cntrl; + + /** + * Memory for auxilliary vectors. + */ + VectorMemory &memory; +}; + + + + + +/** + * Base class for non-symmetric linear solvers. A second entry point + * allows the simultaneous computation of a dual problem. + */ +template +class SolverDual : public Solver +{ + public: + /** + * Constructor. + */ + SolverDual(SolverControl&, VectorMemory&); + + /** + * Solve the original problem + * $Ax=b$. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b) = 0; + + /** + * Solve the two problems + * $Ax=b1$ and $A^Tz=b2$ simultanously. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b1, + Vector &z, + const Vector &b2) = 0; +}; + + + + + + +/*-------------------------------- Inline functions ------------------------*/ + +template +inline +SolverControl & Solver::control() const { + return cntrl; +}; + + + +template +inline +Solver::Solver(SolverControl &cn, VectorMemory &mem) + : cntrl(cn), + memory(mem) +{}; + + + + + +/*---------------------------- solver.h ---------------------------*/ +/* end of #ifndef __solver_H */ +#endif +/*---------------------------- solver.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h new file mode 100644 index 0000000000..167aedf4a6 --- /dev/null +++ b/deal.II/lac/include/lac/solver_cg.h @@ -0,0 +1,160 @@ +/*---------------------------- solver_pcg.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_pcg_H +#define __solver_pcg_H +/*---------------------------- solver_pcg.h ---------------------------*/ + + + +#include +#include + + + + +/** + * Preconditioned cg method. + */ +template +class SolverPCG : public Solver { + public: + + /** + * Constructor. + */ + SolverPCG (SolverControl &cn, VectorMemory &mem) : + Solver(cn,mem) {}; + + /** + * Solver method. + */ + 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(); + + /** + * Temporary vectors, allocated through + * the #VectorMemory# object at the start + * of the actual solution process and + * deallocated at the end. + */ + Vector *Vr, *Vp, *Vz, *VAp; + + /** + * 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 ----------------------------*/ + + +template +double SolverPCG::criterion() +{ + return sqrt(res2); +}; + + + +template +Solver::ReturnState +SolverPCG::solve (const Matrix &A, + Vector &x, + const Vector &b) { + SolverControl::State conv=SolverControl::iterate; + + // Memory allocation + Vr = memory.alloc(); Vector& g = *Vr; + Vp = memory.alloc(); Vector& h = *Vp; + Vz = memory.alloc(); Vector& d = *Vz; + VAp = memory.alloc(); Vector& Ad = *VAp; + + // Implementation taken from the DEAL + // library + int it=0; + double res,gh,alpha,beta; + + res = A.residual(g,x,b); + conv = control().check(0,res); + if (conv) + { + memory.free(Vr); + memory.free(Vp); + memory.free(Vz); + memory.free(VAp); + + return success; + }; + + g.scale(-1.); + A.precondition(h,g); + + d.equ(-1.,h); + + gh = g*h; + + while (conv == SolverControl::iterate) + { + A.vmult(Ad,d); + + alpha = d*Ad; + alpha = gh/alpha; + + g.add(alpha,Ad); + x.add(alpha,d ); + res = sqrt(g*g); + + conv = control().check(it,res); + if (conv) break; + + A.precondition(h,g); + + beta = gh; + gh = g*h; + beta = gh/beta; + + d.sadd(beta,-1.,h); + it++; + }; + + + // Deallocate Memory + + memory.free(Vr); + memory.free(Vp); + memory.free(Vz); + memory.free(VAp); + + // Output + + if (conv == SolverControl::failure) + return exceeded; + else + return success; +}; + + + + +/*---------------------------- solver_pcg.h ---------------------------*/ +/* end of #ifndef __solver_pcg_H */ +#endif +/*---------------------------- solver_pcg.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/solver_control.h b/deal.II/lac/include/lac/solver_control.h new file mode 100644 index 0000000000..444f42b050 --- /dev/null +++ b/deal.II/lac/include/lac/solver_control.h @@ -0,0 +1,189 @@ +/*---------------------------- solver_control.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_control_H +#define __solver_control_H +/*---------------------------- solver_control.h ---------------------------*/ + + + +/** + * Control class for iterative solvers. + * + * Used by iterative methods to + * determine whether the iteration should be continued. To this respect, + * the virtual function #check()# is called in each iteration + * with the current iteration + * step and the value indicating convergence (usually the residual). + * + * After the iteration has terminated, the functions #last_value# and + * #last_step# can be used to obtain information about the final state + * of the iteration. + * + * #check()# can be replaced in derived classes to allow for more + * sophisticated tests. + */ +class SolverControl { + public: + /** + * Return states of the check + * function, which indicate the state the + * solver is in. + * + * The possible values of State are + *
    + *
  1. #iterate = 0#: continue + * the iteration. + *
  2. #success#: the goal is reached, + * the iterative method can terminate + * successfully. + *
  3. #failure#!: the iterative + * method should stop because + * convergence cannot be achieved or at + * least was not achieved within the given + * maximal number of iterations. + *
+ */ + enum State { + iterate = 0, success, failure + }; + + /** + * Constructor. The parameters + * #n# and #tol# are the + * maximum number of iteration + * steps before failure and the + * tolerance to determine success + * of the iteration. */ + SolverControl (const unsigned int n, const double tol); + + /** + * Decide about success or failure + * of an iteration. This function + * gets the current iteration step + * to determine, whether the + * allowed number of steps has + * been exceeded and returns + * #failure# in this case. If + * #check_value# is below the + * prescribed tolerance, it + * returns #success#. In all + * other cases #iterate# is + * returned to suggest + * continuation of the iterative + * procedure. + * + * #check()# additionally + * preserves #step# and + * #check_value#. These + * values are accessible by + * #last_value()# and + * #last_step()#. + * + * Derived classes may overload this + * function, e.g. to log the convergence + * indicators (#check_value#) or to do + * other computations. + */ + virtual State check (const unsigned int step, const double check_value); + + /** + * Return the convergence value of last + * iteration step for which #check# was + * called by the solver. + */ + double last_value() const; + + /** + * Number of last iteration step. + */ + unsigned int last_step() const; + + protected: + /** + * Maximum number of steps. + */ + const unsigned int maxsteps; + + /** + * Prescribed tolerance to be achieved. + */ + const double tol; + + /** + * Last value of the convergence criterion. + */ + double lvalue; + + /** + * Last step. + */ + unsigned int lstep; +}; + + + + +/** + * Specialization of #SolverControl# which returns #success# if either + * the specified tolerance is achieved or if the + * initial residual (or whatever criterion was chosen by the solver + * class) is reduced by a given factor. This is useful in cases where + * you don't want to solve exactly, but rather want to gain two digits. + */ +class ReductionControl : public SolverControl { + public: + /** + * Constructor. Provide the + * reduction factor additional to + * the arguments of the Control + * constructor. + */ + ReductionControl (const unsigned int maxiter, + const double tolerance, + const double reduce); + + /** + * Decide about success or failure + * of an iteration. This function + * calls the one in the base + * class, but sets the tolerance + * to #reduction * initial value# + * upon the first iteration. + */ + virtual State check (const unsigned int step, + const double check_value); + + /** + * Return the initial convergence + * criterion. + */ + double initial_value() const; + + + protected: + /** + * Desired reduction factor. + */ + const double reduce; + + /** + * Initial value. + */ + double initial_val; + + /** + * Reduced tolerance. Stop iterations + * if either this value is achieved + * or if the base class indicates + * success. + */ + double reduced_tol; +}; + + + + +/*---------------------------- solver_control.h ---------------------------*/ +/* end of #ifndef __solver_control_H */ +#endif +/*---------------------------- solver_control.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/vector_memory.h b/deal.II/lac/include/lac/vector_memory.h new file mode 100644 index 0000000000..7b670c2f17 --- /dev/null +++ b/deal.II/lac/include/lac/vector_memory.h @@ -0,0 +1,80 @@ +/*---------------------------- vector_memory.h ---------------------------*/ +/* $Id$ */ +#ifndef __vector_memory_H +#define __vector_memory_H +/*---------------------------- vector_memory.h ---------------------------*/ + + + + +/** + * Memory management for vectors. This class is used by all + * iterative methods to allocate space for auxilliary + * vectors. This class is used to avoid interaction with the + * operating system whenever a vector is needed. Especially, when + * an iterative method is invoked as part of an outer iteration, + * this would lead to runtime overhead and memory fragmentation. + * + * Classes derived from this class implement a more or less + * sophisticated management of vectors. One of these has to be + * applied by the user according to his needs. + */ +template +class VectorMemory { + public: + /** + * Return new vector from the pool. + */ + virtual Vector* alloc() = 0; + + /** + * Return a vector into the pool + * for later use. + */ + virtual void free(const Vector*) = 0; +}; + + + + + +/** + * Simple memory management. This memory class is just made for + * tests. It requires the vector to have a constructor with integer + * parameter indicating the size. It just allocates and deletes + * vectors as needed from the global heap, i.e. performs no + * specially adapted actions to the purpose of this class. + */ +template +class PrimitiveVectorMemory : public VectorMemory { + const unsigned int size; + public: + /** + * Constructor. + * #sz# is the length of vectors to + * be generated by alloc. + */ + PrimitiveVectorMemory (const unsigned int sz) : size(sz) {}; + + /** + * Allocate a vector of the given size + * from the global heap. + */ + virtual Vector* alloc() { + return new Vector(size); + }; + + /** + * Return a vector to the global heap. + */ + virtual void free (const Vector* v) { + delete v; + }; +}; + + + +/*---------------------------- vector_memory.h ---------------------------*/ +/* end of #ifndef __vector_memory_H */ +#endif +/*---------------------------- vector_memory.h ---------------------------*/ diff --git a/deal.II/lac/source/solver_control.cc b/deal.II/lac/source/solver_control.cc new file mode 100644 index 0000000000..26e9dba996 --- /dev/null +++ b/deal.II/lac/source/solver_control.cc @@ -0,0 +1,86 @@ +// $Id$ + + +#include + + + +/*----------------------- SolverControl ---------------------------------*/ + + +SolverControl::SolverControl (const unsigned int maxiter, + const double tolerance) : + maxsteps(maxiter), + tol(tolerance), + lvalue(1.e300), + lstep(0) +{}; + + + +SolverControl::State +SolverControl::check (const unsigned int step, + const double check_value) { + lstep = step; + lvalue = check_value; + if (step>=maxsteps) return failure; + if (check_value <= tol) return success; + return iterate; +}; + + + +double +SolverControl::last_value() const +{ + return lvalue; +}; + + + +unsigned int +SolverControl::last_step() const +{ + return lstep; +}; + + + + +/*----------------------- ReductionControl ---------------------------------*/ + + + +ReductionControl::ReductionControl(const unsigned int n, + const double tol, + const double red) : + SolverControl (n, tol), + reduce(red) +{}; + + + + +double +ReductionControl::initial_value() const { + return initial_val; +}; + + + +SolverControl::State +ReductionControl::check (const unsigned int step, + const double check_value) +{ + if (step==0) + { + initial_val = check_value; + reduced_tol = check_value * reduce; + }; + + if (check_value < reduced_tol) + return success; + else + return SolverControl::check(step, check_value); +}; + -- 2.39.5