cg implemented

ilu0 preconditioner can now be used

[[Imported from SVN: r951]]

istl/istlexception.hh

@@ -5,7 +5,7 @@
 
 #include <stdlib.h>
 
-#include "exceptions.hh"
+#include "dune/common/exceptions.hh"
 
 namespace Dune {
 

istl/preconditioners.hh

@@ -158,6 +158,56 @@ namespace Dune {
 
 
 
+  /*! Wraps the naked ISTL generic ILU0 preconditioner into the
+      solver framework.
+   */
+  template<class M, class X, class Y>
+  class SeqILU0 : 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.
+    SeqILU0 (const M& A)
+      : ILU(A)     // copy A
+    {
+      bilu0_decomposition(ILU);
+    }
+
+    //! 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)
+    {
+      bilu_backsolve(ILU,v,d);
+    }
+
+    //! 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:
+    M ILU;
+  };
+
+
+
   /** @} end documentation */
 
 } // end namespace

istl/solvers.hh

@@ -13,7 +13,7 @@
 #include "istlexception.hh"
 #include "operators.hh"
 #include "preconditioners.hh"
-#include "timer.hh"
+#include "dune/common/timer.hh"
 
 /*! \file __FILE__
 
@@ -193,6 +193,95 @@ namespace Dune {
   };
 
 
+
+
+  //! conjugate gradient method
+  template<class X, class Y>
+  class CGSolver : 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
+    CGSolver (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)
+    {
+      r.clear();                      // clear solver statistics
+      Timer watch;                    // start a timer
+      _prec.pre(x,b);                 // prepare preconditioner
+      _op.applyscaleadd(-1,x,b);      // overwrite b with defect
+
+      X p(x);                         // create local vectors
+      Y q(b);
+
+      double def0 = _prec.norm(b);    // compute norm
+
+      if (_verbose>0)                 // printing
+      {
+        printf("=== CGSolver\n");
+        if (_verbose>1) printf(" Iter       Defect         Rate\n");
+        if (_verbose>1) printf("%5d %12.4E\n",0,def0);
+      }
+
+      int i=1; double def=def0;       // loop variables
+      field_type rho,rholast,lambda;
+      for ( ; i<=_maxit; i++ )
+      {
+        p = 0;                                // clear correction
+        _prec.apply(p,b);                     // apply preconditioner
+        rho = _prec.dot(p,b);                 // orthogonalization
+        if (i>1) p.axpy(rho/rholast,q);
+        rholast = rho;                        // remember rho for recurrence
+        _op.apply(p,q);                       // q=Ap
+        lambda = rho/_prec.dot(q,p);          // minimization
+        x.axpy(lambda,p);                     // update solution
+        b.axpy(-lambda,q);                    // update defect
+        q=p;                                  // remember search direction
+        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;
+        }
+      }
+
+      if (_verbose==1)                    // printing for non verbose
+        printf("%5d %12.4E\n",i,def);
+      _prec.post(x);                      // postprocess preconditioner
+      r.iterations = i;                   // fill statistics
+      r.reduction = def/def0;
+      r.conv_rate  = pow(r.reduction,1.0/i);
+      r.elapsed = watch.elapsed();
+      if (_verbose>0)                     // final print
+        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

istl/tutorial/example.cc

@@ -435,7 +435,7 @@ void test_Iter ()
 void test_Interface ()
 {
   // define Types
-  const int BlockSize = 6;
+  const int BlockSize = 1;
   typedef Dune::FieldVector<double,BlockSize> VB;
   typedef Dune::FieldMatrix<double,BlockSize,BlockSize> MB;
   typedef Dune::BlockVector<VB> Vector;
@@ -452,7 +452,7 @@ void test_Interface ()
     E[i][i] = -1;
 
   // make a block compressed row matrix with five point stencil
-  const int N=10000, BW1=1, BW2=100;
+  const int BW2=512, BW1=1, N=BW2*BW2;
   Matrix A(N,N,5*N,Dune::BCRSMatrix<MB>::row_wise);
   for (Matrix::CreateIterator i=A.createbegin(); i!=A.createend(); ++i)
   {
@@ -477,8 +477,9 @@ void test_Interface ()
 
   // set up the high-level solver objects
   Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);        // make linear operator from A
-  Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,2,1.0);      // preconditioner object
-  Dune::LoopSolver<Vector,Vector> solver(op,ssor,1E-4,50,2); // an inverse operator
+  Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,1.78);     // SSOR preconditioner
+  Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A);            // preconditioner object
+  Dune::CGSolver<Vector,Vector> solver(op,ssor,1E-8,8000,2); // an inverse operator
 
   // call the solver
   Dune::InverseOperatorResult r;