Und hier ist die Krylovmethode für heute: BiCGStab

Ist noch nicht ausführlich getestet (wie alles hier), liefert
aber erst mal ähnliche Ergebnisse wie CG

[[Imported from SVN: r961]]

istl/solvers.hh

@@ -193,25 +193,25 @@ namespace Dune {
   };
 
 
-
+  // all these solvers are taken from the SUMO library
 
   //! conjugate gradient method
-  template<class X, class Y>
-  class CGSolver : public InverseOperator<X,Y> {
+  template<class X>
+  class CGSolver : public InverseOperator<X,X> {
   public:
     //! export types, usually they come from the derived class
     typedef X domain_type;
-    typedef Y range_type;
+    typedef X range_type;
     typedef typename X::field_type field_type;
 
     //! set up Loop solver
-    CGSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec,
+    CGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& 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)
+    virtual void apply (X& x, X& b, InverseOperatorResult& r)
     {
       r.clear();                      // clear solver statistics
       Timer watch;                    // start a timer
@@ -219,7 +219,7 @@ namespace Dune {
       _op.applyscaleadd(-1,x,b);      // overwrite b with defect
 
       X p(x);                         // create local vectors
-      Y q(b);
+      X q(b);
 
       double def0 = _prec.norm(b);    // compute norm
 
@@ -267,15 +267,219 @@ namespace Dune {
     }
 
     //! apply inverse operator, with given reduction factor
-    virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& r)
+    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& r)
     {
       _reduction = reduction;
       (*this).apply(x,b,r);
     }
 
   private:
-    LinearOperator<X,Y>& _op;
-    Preconditioner<X,Y>& _prec;
+    LinearOperator<X,X>& _op;
+    Preconditioner<X,X>& _prec;
+    double _reduction;
+    int _maxit;
+    int _verbose;
+  };
+
+
+  // Ronald Kriemanns BiCG-STAB implementation from Sumo
+  //! Bi-conjugate Gradient Stabilized (BiCG-STAB)
+  template<class X>
+  class BiCGSTABSolver : public InverseOperator<X,X> {
+  public:
+    //! export types, usually they come from the derived class
+    typedef X domain_type;
+    typedef X range_type;
+    typedef typename X::field_type field_type;
+
+    //! set up Loop solver
+    BiCGSTABSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& 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, X& b, InverseOperatorResult& res)
+    {
+      const double EPSILON=1e-20;
+
+      int it, me;
+      field_type rho, rho_new, alpha, beta, h, omega;
+      field_type norm, norm_old, norm_0;
+
+      //
+      // get vectors and matrix
+      //
+      X& r=b;
+      X p(x);
+      X v(x);
+      X t(x);
+      X y(x);
+      X rt(x);
+
+      //
+      // begin iteration
+      //
+
+      // r = r - Ax; rt = r
+      res.clear();                      // clear solver statistics
+      Timer watch;                    // start a timer
+      _op.applyscaleadd(-1,x,r);      // overwrite b with defect
+      _prec.pre(x,r);                 // prepare preconditioner
+
+      rt=r;
+
+      norm = norm_old = norm_0 = r.two_norm();
+
+      p=0;
+      v=0;
+
+      rho   = 1;
+      alpha = 1;
+      omega = 1;
+
+      if (_verbose>0)                 // printing
+      {
+        printf("=== BiCGSTABSolver\n");
+        if (_verbose>1) printf(" Iter       Defect         Rate\n");
+        if (_verbose>1) printf("%5d %12.4E\n",0,norm_0);
+      }
+
+      //
+      // iteration
+      //
+
+      it = 0;
+
+      while ( true )
+      {
+        //
+        // preprocess, set vecsizes etc.
+        //
+
+        // rho_new = < rt , r >
+        rho_new = _prec.dot(rt,r);
+
+        // look if breakdown occured
+        if ((std::abs(rho) < EPSILON) || (std::abs(omega) < EPSILON))
+          DUNE_THROW(ISTLError,"breakdown in BiCGSTAB");
+
+        if (it==0)
+          p = r;
+        else
+        {
+          beta = ( rho_new / rho ) * ( alpha / omega );
+          p.axpy(-omega,v);                 // p = r + beta (p - omega*v)
+          p *= beta;
+          p += r;
+        }
+
+        // y = W^-1 * p
+        y = 0;
+        _prec.apply(y,p);                     // apply preconditioner
+
+        // v = A * y
+        _op.apply(y,v);
+
+        // alpha = rho_new / < rt, v >
+        h = _prec.dot(rt,v);
+
+        if ( std::abs(h) < EPSILON )
+          DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
+
+        alpha = rho_new / h;
+
+        // apply first correction to x
+        // x <- x + alpha y
+        x.axpy(alpha,y);
+
+        // r = r - alpha*v
+        r.axpy(-alpha,v);
+
+        //
+        // test stop criteria
+        //
+
+        it++;
+        norm = r.two_norm();
+
+        if (_verbose>1)                       // print
+          printf("%5d %12.4E %12.4g\n",it,norm,norm/norm_old);
+
+        if ( norm < (_reduction * norm_0) )
+        {
+          res.converged = 1;
+          break;
+        }
+
+        if (it >= _maxit)
+          break;
+
+        norm_old = norm;
+
+        // y = W^-1 * r
+        y = 0;
+        _prec.apply(y,r);
+
+        // t = A * y
+        _op.apply(y,t);
+
+        // omega = < t, r > / < t, t >
+        omega = _prec.dot(t,r)/_prec.dot(t,t);
+
+        // apply second correction to x
+        // x <- x + omega y
+        x.axpy(omega,y);
+
+        // r = s - omega*t (remember : r = s)
+        r.axpy(-omega,t);
+
+        rho = rho_new;
+
+        //
+        // test stop criteria
+        //
+
+        it++;
+
+        norm = r.two_norm();
+
+        if (_verbose>1)                       // print
+          printf("%5d %12.4E %12.4g\n",it,norm,norm/norm_old);
+
+        if ( norm < (_reduction * norm_0) )
+        {
+          res.converged = 1;
+          break;
+        }
+
+        if (it >= _maxit)
+          break;
+
+        norm_old = norm;
+      }          // while
+
+      if (_verbose==1)                    // printing for non verbose
+        printf("%5d %12.4E\n",it,norm);
+      _prec.post(x);                      // postprocess preconditioner
+      res.iterations = it;                  // fill statistics
+      res.reduction = norm/norm_0;
+      res.conv_rate  = pow(res.reduction,1.0/it);
+      res.elapsed = watch.elapsed();
+      if (_verbose>0)                     // final print
+        printf("=== rate=%g, T=%g, TIT=%g\n",res.conv_rate,res.elapsed,res.elapsed/it);
+    }
+
+    //! apply inverse operator, with given reduction factor
+    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& r)
+    {
+      _reduction = reduction;
+      (*this).apply(x,b,r);
+    }
+
+  private:
+    LinearOperator<X,X>& _op;
+    Preconditioner<X,X>& _prec;
     double _reduction;
     int _maxit;
     int _verbose;

istl/tutorial/example.cc

@@ -452,7 +452,7 @@ void test_Interface ()
     E[i][i] = -1;
 
   // make a block compressed row matrix with five point stencil
-  const int BW2=512, BW1=1, N=BW2*BW2;
+  const int BW2=64, 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)
   {
@@ -479,11 +479,12 @@ void test_Interface ()
   Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);        // make linear operator from A
   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
+  Dune::CGSolver<Vector> cg(op,ilu0,1E-8,8000,2);         // an inverse operator
+  Dune::BiCGSTABSolver<Vector> bcgs(op,ilu0,1E-8,8000,2); // an inverse operator
 
   // call the solver
   Dune::InverseOperatorResult r;
-  solver.apply(x,b,r);
+  cg.apply(x,b,r);
 }