]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the iterative solver classes and helpers from the indepenant MIA library to...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 3 Aug 1998 11:29:35 +0000 (11:29 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 3 Aug 1998 11:29:35 +0000 (11:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@469 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver.h [new file with mode: 0644]
deal.II/lac/include/lac/solver_cg.h [new file with mode: 0644]
deal.II/lac/include/lac/solver_control.h [new file with mode: 0644]
deal.II/lac/include/lac/vector_memory.h [new file with mode: 0644]
deal.II/lac/source/solver_control.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/solver.h b/deal.II/lac/include/lac/solver.h
new file mode 100644 (file)
index 0000000..687ad83
--- /dev/null
@@ -0,0 +1,235 @@
+/*----------------------------   solver.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __solver_H
+#define __solver_H
+/*----------------------------   solver.h     ---------------------------*/
+
+
+
+// forward declaration
+class SolverControl;
+template <class Vector> 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 Matrix, class Vector>
+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<Vector> &);
+    
+                                    /**
+                                     * 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<Vector> &memory;
+};
+
+
+
+
+
+/**
+ * Base class for non-symmetric linear solvers.  A second entry point
+ * allows the simultaneous computation of a dual problem.
+ */
+template<class Matrix, class Vector>
+class SolverDual : public Solver<Matrix, Vector>
+{
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverDual(SolverControl&, VectorMemory<Vector>&);
+    
+                                    /**
+                                     * 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 <class Matrix, class Vector>
+inline
+SolverControl & Solver<Matrix,Vector>::control() const {
+  return cntrl;
+};
+
+
+
+template<class Matrix, class Vector>
+inline
+Solver<Matrix, Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &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 (file)
index 0000000..167aedf
--- /dev/null
@@ -0,0 +1,160 @@
+/*----------------------------   solver_pcg.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __solver_pcg_H
+#define __solver_pcg_H
+/*----------------------------   solver_pcg.h     ---------------------------*/
+
+
+
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+
+
+
+
+/**
+ * Preconditioned cg method.
+ */
+template<class Matrix, class Vector>
+class SolverPCG : public Solver<Matrix,Vector> {
+  public:
+
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverPCG (SolverControl &cn, VectorMemory<Vector> &mem) :
+                   Solver<Matrix,Vector>(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<class Matrix, class Vector>
+double SolverPCG<Matrix,Vector>::criterion()
+{
+  return sqrt(res2);
+};
+
+
+
+template<class Matrix, class Vector>
+Solver<Matrix,Vector>::ReturnState 
+SolverPCG<Matrix,Vector>::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 (file)
index 0000000..444f42b
--- /dev/null
@@ -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
+                                     * <OL>
+                                     * <LI> #iterate = 0#: continue
+                                     * the iteration.
+                                     * <LI> #success#: the goal is reached,
+                                     * the iterative method can terminate
+                                     * successfully.
+                                     * <LI> #failure#!: the iterative
+                                     * method should stop because
+                                     * convergence cannot be achieved or at
+                                     * least was not achieved within the given
+                                     * maximal number of iterations.
+                                     * </OL>
+                                     */
+    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 (file)
index 0000000..7b670c2
--- /dev/null
@@ -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 Vector>
+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 Vector>
+class PrimitiveVectorMemory : public VectorMemory<Vector> {
+    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 (file)
index 0000000..26e9dba
--- /dev/null
@@ -0,0 +1,86 @@
+// $Id$
+
+
+#include <lac/solver_control.h>
+
+
+
+/*----------------------- 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);
+};
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.