The ISTL high-level solver interface

- supports operator concept, i.e. you can use e.g. CG with on-the-fly operators
  and preconditioners
- supports parallel preconditioners

[[Imported from SVN: r923]]

istl/operators.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_ISTLOPERATORS_HH__
+#define __DUNE_ISTLOPERATORS_HH__
+
+#include <math.h>
+#include <complex>
+#include <iostream>
+#include <iomanip>
+#include <string>
+
+
+/*! \file __FILE__
+
+    Define general, extensible interface for operators.
+        The available implementation wraps a matrix.
+ */
+
+namespace Dune {
+
+  /** @addtogroup ISTL
+          @{
+   */
+
+
+  //=====================================================================
+  // Abstract operator interface
+  //=====================================================================
+
+  /*! Abstract base class defining an operator \f$ A : X\to Y\f$. The
+     simplest solvers just need the application  \f$ A(x)\f$ of
+     the operator. The operator might even be nonlinear (but is not in
+     our application here).
+
+        - enables on the fly computation through operator concept. If explicit
+        representation of the operator is required use AssembledLinearOperator
+
+        - Some inverters may need an explicit formation of the operator
+        as a matrix, e.g. BiCGStab, ILU, AMG, etc. In that case use the
+        derived class
+   */
+  template<class X, class Y>
+  class Operator {
+  public:
+    //! export types, they come from the derived class
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! apply operator to x:  \f$ y = A(x) \f$
+    virtual void apply (const X& x, Y& y) const = 0;
+
+    //! every abstract base class has a virtual destructor
+    virtual ~Operator () {}
+  };
+
+
+  /*! This type ensures that the operator \f$ A \f$ is a
+      linear operator, i.e. \f$ A(\alpha x) = \alpha A(x) \f$ and
+      \f$ A(x+y) = A(x)+A(y)\f$. The additional interface function
+      reflects this fact.
+   */
+  template<class X, class Y>
+  class LinearOperator : public Operator<X,Y> {
+  public:
+    //! export types, they come from the derived class
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
+    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const = 0;
+
+    //! every abstract base class has a virtual destructor
+    virtual ~LinearOperator () {}
+  };
+
+
+  /*! Linear Operator that exports the operator in
+     matrix form. This is needed for certain solvers, such as
+     LU decomposition, ILU preconditioners or BiCG-Stab (because
+     of multiplication with A^T).
+   */
+  template<class M, class X, class Y>
+  class AssembledLinearOperator : public LinearOperator<X,Y> {
+  public:
+    //! export types, usually they come from the derived class
+    typedef M matrix_type;
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! get matrix via *
+    virtual const M& getmat () const = 0;
+
+    //! every abstract base class has a virtual destructor
+    virtual ~AssembledLinearOperator () {}
+  };
+
+
+
+  //=====================================================================
+  // Implementation for ISTL-matrix based operator
+  //=====================================================================
+
+  //! Adapts a matrix to the assembled linear operator interface
+  template<class M, class X, class Y>
+  class MatrixAdapter : public AssembledLinearOperator<M,X,Y>
+  {
+  public:
+    //! export types
+    typedef M matrix_type;
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! constructor: just store a reference to a matrix
+    MatrixAdapter (const M& A) : _A(A) {}
+
+    //! apply operator to x:  \f$ y = A(x) \f$
+    virtual void apply (const X& x, Y& y) const
+    {
+      y = 0;
+      _A.umv(x,y);
+    }
+
+    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
+    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+    {
+      _A.usmv(alpha,x,y);
+    }
+
+    //! get matrix via *
+    virtual const M& getmat () const
+    {
+      return _A;
+    }
+
+  private:
+    const M& _A;
+  };
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif

istl/preconditioners.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_PRECONDITIONERS_HH__
+#define __DUNE_PRECONDITIONERS_HH__
+
+#include <math.h>
+#include <complex>
+#include <iostream>
+#include <iomanip>
+#include <string>
+
+#include "istlexception.hh"
+#include "io.hh"
+#include "gsetc.hh"
+#include "ilu.hh"
+
+/*! \file __FILE__
+
+    Define general preconditioner interface and wrap the
+    methods implemented by ISTL in this interface.
+        However, the interface is extensible such that new preconditioners
+    can be implemented and used with the solvers.
+ */
+
+namespace Dune {
+
+  /** @addtogroup ISTL
+          @{
+   */
+
+
+  //=====================================================================
+  /*! Base class for matrix free definition of preconditioners.
+          Note that the operator, which is the basis for the preconditioning,
+      is supplied to the preconditioner from the outside in the
+      constructor or some other method.
+
+          This interface allows the encapsulation of all parallelization
+          aspects into the preconditioners.
+   */
+  //=====================================================================
+  template<class X, class Y>
+  class Preconditioner {
+  public:
+    //! export types, they come from the derived class
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    /*! A solver solves a linear operator equation A(x)=b by applying
+       one or several steps of the preconditioner. The method pre()
+       is called before before the first apply operation.
+       x and b are right hand side and solution vector of the linear
+       system. It may. e.g., scale the system, allocate memory or
+       compute a (I)LU decomposition.
+       Note: The ILU decomposition could also be computed in the constructor
+       or with a seperate method of the derived method if several
+       linear systems with the same matrix are to be solved.
+     */
+    virtual void pre (X& x, Y& b) = 0;
+
+    /*! Apply one step of the preconditioner to the system A(v)=d.
+       On entry v=0 and d=b-A(x) (although this might not be
+       computed in that way. On exit v contains the update, i.e
+       one step computes \f$ v = M^{-1} d \f$ where \f$ M \f$ is the
+       approximate inverse of the operator \f$ A \f$ characterizing
+       the preconditioner.
+     */
+    virtual void apply (X& v, const Y& d) = 0;
+
+    /*! Dot product of two right-hand side vectors. This
+       method is in the interface for parallel implementations.
+            It will require at least a global communication and also
+            some local communication depending on the consistency model.
+            It allows the solver to be independent of parallelization
+            issues.
+     */
+    virtual field_type dot (const Y& y, const Y& z) = 0;
+
+    /*! Norm of a right-hand side vector. This
+       method is in the interface for parallel implementations.
+            It will require at least a global communication and also
+            some local communication depending on the consistency model.
+            It allows the solver to be independent of parallelization
+            issues.
+            You may also subclass to overload your favourite norm :-)
+     */
+    virtual double norm (const Y& y) = 0;
+
+
+    /*! This method is called after the last apply call for the
+       linear system to be solved. Memory may be deallocated safely
+       here. x is the solution of the linear equation.
+     */
+    virtual void post (X& x) = 0;
+
+    //! every abstract base class has a virtual destructor
+    virtual ~Preconditioner () {}
+  };
+
+
+  //=====================================================================
+  // Implementation of this interface for sequential ISTL-preconditioners
+  //=====================================================================
+
+
+  /*! Wraps the naked ISTL generic SOR preconditioner into the
+      solver framework.
+   */
+  template<class M, class X, class Y>
+  class SeqSSOR : public Preconditioner<X,Y> {
+  public:
+    //! export types, they come from the derived class
+    typedef M matrix_type;
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! constructor gets all parameters to operate the prec.
+    SeqSSOR (const M& A, int n, field_type w)
+      : _A(A), _n(n), _w(w)
+    {       }
+
+    //! prepare: nothing to do here
+    virtual void pre (X& x, Y& b) {}
+
+    //! just calls the istl functions
+    virtual void apply (X& v, const Y& d)
+    {
+      for (int i=0; i<_n; i++)
+      {
+        bsorf(_A,v,d,_w);
+        bsorb(_A,v,d,_w);
+      }
+    }
+
+    //! sequential case: just call vector function
+    virtual field_type dot (const Y& y, const Y& z)
+    {
+      return y*z;
+    }
+
+    //! sequential case: just call vector function
+    virtual double norm (const Y& y)
+    {
+      return y.two_norm();     // my favourite norm
+    }
+
+
+    // nothing to do here
+    virtual void post (X& x) {}
+
+  private:
+    const M& _A;       // my matrix to operate on
+    int _n;            // number of iterations
+    field_type _w;     // relaxation factor
+  };
+
+
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif

istl/solvers.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_SOLVERS_HH__
+#define __DUNE_SOLVERS_HH__
+
+#include <math.h>
+#include <complex>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <stdio.h> // there is nothing better than printf
+
+#include "istlexception.hh"
+#include "operators.hh"
+#include "preconditioners.hh"
+#include "timer.hh"
+
+/*! \file __FILE__
+
+    Define general, extensible interface for inverse operators.
+        Implementation here coverse only inversion of linear operators,
+        but the implementation might be used for nonlinear operators
+    as well.
+ */
+
+namespace Dune {
+
+  /** @addtogroup ISTL
+          @{
+   */
+
+
+  /*! The return value of an application of the inverse
+      operator delivers some important information about
+          the iteration.
+   */
+  struct InverseOperatorResult
+  {
+    InverseOperatorResult ()
+    {
+      clear();
+    }
+
+    void clear ()
+    {
+      iterations = 0;
+      reduction = 0;
+      converged = false;
+      conv_rate = 1;
+      elapsed = 0;
+    }
+
+    int iterations;       // number of iterations
+    double reduction;     // reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$
+    bool converged;       // true if convergence criterion has been met
+    double conv_rate;     // convergence rate (average reduction per step)
+    double elapsed;       // elapsed time in seconds
+  };
+
+
+  //=====================================================================
+  /*! An InverseOperator computes the solution of \f$ A(x)=b\f$ where
+      \f$ A : X \to Y \f$ is an operator.
+      Note that the solver "knows" which operator
+      to invert and which preconditioner to apply (if any). The
+      user is only interested in inverting the operator.
+          InverseOperator might be a Newton scheme, a Krylov subspace method,
+      or a direct solver or just anything.
+   */
+  //=====================================================================
+  template<class X, class Y>
+  class InverseOperator {
+  public:
+    //! export types, usually they come from the derived class
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! apply inverse operator, Note: right hand side b may be overwritten!
+    virtual void apply (X& x, Y& b, InverseOperatorResult& r) = 0;
+
+    //! apply inverse operator, with given convergence criteria, right hand side b may be overwritten!
+    virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& r) = 0;
+
+    //! the usual thing
+    virtual ~InverseOperator () {}
+  };
+
+
+  //=====================================================================
+  // Implementation of this interface
+  //=====================================================================
+
+  /*! Implements a preconditioned loop.
+          Verbose levels are:
+          0 : print nothing
+          1 : print initial and final defect and statistics
+          2 : print line for each iteration
+   */
+  template<class X, class Y>
+  class LoopSolver : public InverseOperator<X,Y> {
+  public:
+    //! export types, usually they come from the derived class
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! set up Loop solver
+    LoopSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec,
+                double reduction, int maxit, int verbose) :
+      _op(op), _prec(prec), _reduction(reduction), _maxit(maxit), _verbose(verbose)
+    {}
+
+    //! apply inverse operator
+    virtual void apply (X& x, Y& b, InverseOperatorResult& r)
+    {
+      // clear solver statistics
+      r.clear();
+
+      // start a timer
+      Timer watch;
+
+      // prepare preconditioner
+      _prec.pre(x,b);
+
+      // overwrite b with defect
+      _op.applyscaleadd(-1,x,b);
+
+      // compute norm, \todo parallelization
+      double def0 = _prec.norm(b);
+
+      // printing
+      if (_verbose>0)
+      {
+        printf("=== LoopSolver\n");
+        if (_verbose>1) printf(" Iter       Defect         Rate\n");
+        if (_verbose>1) printf("%5d %12.4E\n",0,def0);
+      }
+
+      // allocate correction vector
+      X v(x);
+
+      // iteration loop
+      int i=1; double def=def0;
+      for ( ; i<=_maxit; i++ )
+      {
+        v = 0;                                // clear correction
+        _prec.apply(v,b);                     // apply preconditioner
+        x += v;                               // update solution
+        _op.applyscaleadd(-1,v,b);            // update defect
+        double defnew=_prec.norm(b);          // comp defect norm
+        if (_verbose>1)                       // print
+          printf("%5d %12.4E %12.4g\n",i,defnew,defnew/def);
+        def = defnew;                         // update norm
+        if (def<def0*_reduction)              // convergence check
+        {
+          r.converged  = true;
+          break;
+        }
+      }
+
+      // print
+      if (_verbose==1)
+        printf("%5d %12.4E\n",i,def);
+
+      // postprocess preconditioner
+      _prec.post(x);
+
+      // fill statistics
+      r.iterations = i;
+      r.reduction = def/def0;
+      r.conv_rate  = pow(r.reduction,1.0/i);
+      r.elapsed = watch.elapsed();
+
+      // final print
+      if (_verbose>0)
+        printf("=== rate=%g, T=%g, TIT=%g\n",r.conv_rate,r.elapsed,r.elapsed/i);
+    }
+
+    //! apply inverse operator, with given reduction factor
+    virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& r)
+    {
+      _reduction = reduction;
+      (*this).apply(x,b,r);
+    }
+
+  private:
+    LinearOperator<X,Y>& _op;
+    Preconditioner<X,Y>& _prec;
+    double _reduction;
+    int _maxit;
+    int _verbose;
+  };
+
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif