dune-istl  2.2.0
solvers.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set ts=4 sw=2 et sts=2:
00003 
00004 #ifndef DUNE_SOLVERS_HH
00005 #define DUNE_SOLVERS_HH
00006 
00007 #include<cmath>
00008 #include<complex>
00009 #include<iostream>
00010 #include<iomanip>
00011 #include<string>
00012 
00013 #include "istlexception.hh"
00014 #include "operators.hh"
00015 #include "preconditioners.hh"
00016 #include "scalarproducts.hh"
00017 #include <dune/common/timer.hh>
00018 #include <dune/common/ftraits.hh>
00019 #include <dune/common/static_assert.hh>
00020 
00021 namespace Dune {
00046   struct InverseOperatorResult
00047   {
00049     InverseOperatorResult ()
00050     {
00051       clear();
00052     }
00053 
00055     void clear ()
00056     {
00057       iterations = 0;
00058       reduction = 0;
00059       converged = false;
00060       conv_rate = 1;
00061       elapsed = 0;
00062     }
00063 
00065     int iterations;
00066 
00068     double reduction;
00069 
00071     bool converged;
00072 
00074     double conv_rate;
00075 
00077     double elapsed;
00078   };
00079 
00080 
00081   //=====================================================================
00093   template<class X, class Y>
00094   class InverseOperator {
00095   public:
00097     typedef X domain_type;
00098 
00100     typedef Y range_type;
00101 
00103     typedef typename X::field_type field_type;
00104 
00114     virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
00115 
00126     virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
00127 
00129     virtual ~InverseOperator () {}
00130 
00131   protected:
00132     // spacing values
00133     enum { iterationSpacing = 5 , normSpacing = 16 };
00134 
00136     void printHeader(std::ostream& s) const
00137     {
00138       s << std::setw(iterationSpacing)  << " Iter";
00139       s << std::setw(normSpacing) << "Defect";
00140       s << std::setw(normSpacing) << "Rate" << std::endl;
00141     }
00142 
00144     template <class DataType>
00145     void printOutput(std::ostream& s,
00146       const double iter,
00147       const DataType& norm,
00148       const DataType& norm_old) const
00149     {
00150       const DataType rate = norm/norm_old;
00151       s << std::setw(iterationSpacing)  << iter << " ";
00152       s << std::setw(normSpacing) << norm << " ";
00153       s << std::setw(normSpacing) << rate << std::endl;
00154     }
00155 
00157     template <class DataType>
00158     void printOutput(std::ostream& s,
00159       const double iter,
00160       const DataType& norm) const
00161     {
00162       s << std::setw(iterationSpacing)  << iter << " ";
00163       s << std::setw(normSpacing) << norm << std::endl;
00164     }
00165   };
00166 
00167 
00168   //=====================================================================
00169   // Implementation of this interface
00170   //=====================================================================
00171 
00180   template<class X>
00181   class LoopSolver : public InverseOperator<X,X> {
00182   public:
00184     typedef X domain_type;
00186     typedef X range_type;
00188     typedef typename X::field_type field_type;
00189 
00209     template<class L, class P>
00210     LoopSolver (L& op, P& prec,
00211       double reduction, int maxit, int verbose) :
00212       ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00213     {
00214       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
00215         "L and P have to have the same category!");
00216       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
00217         "L has to be sequential!");
00218     }
00219 
00240     template<class L, class S, class P>
00241     LoopSolver (L& op, S& sp, P& prec,
00242       double reduction, int maxit, int verbose) :
00243       _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00244     {
00245       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00246         "L and P must have the same category!");
00247       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
00248         "L and S must have the same category!");
00249     }
00250 
00251 
00253     virtual void apply (X& x, X& b, InverseOperatorResult& res)
00254     {
00255       // clear solver statistics
00256       res.clear();
00257 
00258       // start a timer
00259       Timer watch;
00260 
00261       // prepare preconditioner
00262       _prec.pre(x,b);
00263 
00264       // overwrite b with defect
00265       _op.applyscaleadd(-1,x,b);
00266 
00267       // compute norm, \todo parallelization
00268       double def0 = _sp.norm(b);
00269 
00270       // printing
00271       if (_verbose>0)
00272       {
00273         std::cout << "=== LoopSolver" << std::endl;
00274         if (_verbose>1)
00275         {
00276           this->printHeader(std::cout);
00277           this->printOutput(std::cout,0,def0);
00278         }
00279       }
00280 
00281       // allocate correction vector
00282       X v(x);
00283 
00284       // iteration loop
00285       int i=1; double def=def0;
00286       for ( ; i<=_maxit; i++ )
00287       {
00288         v = 0;                      // clear correction
00289         _prec.apply(v,b);           // apply preconditioner
00290         x += v;                     // update solution
00291         _op.applyscaleadd(-1,v,b);  // update defect
00292         double defnew=_sp.norm(b);  // comp defect norm
00293         if (_verbose>1)             // print
00294           this->printOutput(std::cout,i,defnew,def);
00295         //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
00296         def = defnew;               // update norm
00297         if (def<def0*_reduction || def<1E-30)    // convergence check
00298         {
00299           res.converged  = true;
00300           break;
00301         }
00302       }
00303 
00304       // print
00305       if (_verbose==1)
00306         this->printOutput(std::cout,i,def);
00307 
00308       // postprocess preconditioner
00309       _prec.post(x);
00310 
00311       // fill statistics
00312       res.iterations = i;
00313       res.reduction = def/def0;
00314       res.conv_rate  = pow(res.reduction,1.0/i);
00315       res.elapsed = watch.elapsed();
00316 
00317       // final print
00318       if (_verbose>0)
00319       {
00320         std::cout << "=== rate=" << res.conv_rate
00321                   << ", T=" << res.elapsed
00322                   << ", TIT=" << res.elapsed/i
00323                   << ", IT=" << i << std::endl;
00324       }
00325     }
00326 
00328     virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
00329     {
00330       std::swap(_reduction,reduction);
00331       (*this).apply(x,b,res);
00332       std::swap(_reduction,reduction);
00333     }
00334 
00335   private:
00336     SeqScalarProduct<X> ssp;
00337     LinearOperator<X,X>& _op;
00338     Preconditioner<X,X>& _prec;
00339     ScalarProduct<X>& _sp;
00340     double _reduction;
00341     int _maxit;
00342     int _verbose;
00343   };
00344 
00345 
00346   // all these solvers are taken from the SUMO library
00348   template<class X>
00349   class GradientSolver : public InverseOperator<X,X> {
00350   public:
00352     typedef X domain_type;
00354     typedef X range_type;
00356     typedef typename X::field_type field_type;
00357 
00363     template<class L, class P>
00364     GradientSolver (L& op, P& prec,
00365       double reduction, int maxit, int verbose) :
00366       ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00367     {
00368       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
00369         "L and P have to have the same category!");
00370       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
00371         "L has to be sequential!");
00372     }
00378     template<class L, class S, class P>
00379     GradientSolver (L& op, S& sp, P& prec,
00380       double reduction, int maxit, int verbose) :
00381       _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00382     {
00383       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
00384         "L and P have to have the same category!");
00385       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
00386         "L and S have to have the same category!");
00387     }
00388 
00394     virtual void apply (X& x, X& b, InverseOperatorResult& res)
00395     {
00396       res.clear();                  // clear solver statistics
00397       Timer watch;                // start a timer
00398       _prec.pre(x,b);             // prepare preconditioner
00399       _op.applyscaleadd(-1,x,b);  // overwrite b with defect
00400 
00401       X p(x);                     // create local vectors
00402       X q(b);
00403 
00404       double def0 = _sp.norm(b);// compute norm
00405 
00406       if (_verbose>0)             // printing
00407       {
00408         std::cout << "=== GradientSolver" << std::endl;
00409         if (_verbose>1)
00410         {
00411           this->printHeader(std::cout);
00412           this->printOutput(std::cout,0,def0);
00413         }
00414       }
00415 
00416       int i=1; double def=def0;   // loop variables
00417       field_type lambda;
00418       for ( ; i<=_maxit; i++ )
00419       {
00420         p = 0;                      // clear correction
00421         _prec.apply(p,b);           // apply preconditioner
00422         _op.apply(p,q);             // q=Ap
00423         lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization
00424         x.axpy(lambda,p);           // update solution
00425         b.axpy(-lambda,q);          // update defect
00426 
00427         double defnew=_sp.norm(b);// comp defect norm
00428         if (_verbose>1)             // print
00429           this->printOutput(std::cout,i,defnew,def);
00430 
00431         def = defnew;               // update norm
00432         if (def<def0*_reduction || def<1E-30)    // convergence check
00433         {
00434           res.converged  = true;
00435           break;
00436         }
00437       }
00438 
00439       if (_verbose==1)                // printing for non verbose
00440         this->printOutput(std::cout,i,def);
00441 
00442       _prec.post(x);                  // postprocess preconditioner
00443       res.iterations = i;               // fill statistics
00444       res.reduction = def/def0;
00445       res.conv_rate  = pow(res.reduction,1.0/i);
00446       res.elapsed = watch.elapsed();
00447       if (_verbose>0)                 // final print
00448         std::cout << "=== rate=" << res.conv_rate
00449                   << ", T=" << res.elapsed
00450                   << ", TIT=" << res.elapsed/i
00451                   << ", IT=" << i << std::endl;
00452     }
00453 
00459     virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
00460     {
00461       std::swap(_reduction,reduction);
00462       (*this).apply(x,b,res);
00463       std::swap(_reduction,reduction);
00464     }
00465 
00466   private:
00467     SeqScalarProduct<X> ssp;
00468     LinearOperator<X,X>& _op;
00469     Preconditioner<X,X>& _prec;
00470     ScalarProduct<X>& _sp;
00471     double _reduction;
00472     int _maxit;
00473     int _verbose;
00474   };
00475 
00476 
00477 
00479   template<class X>
00480   class CGSolver : public InverseOperator<X,X> {
00481   public:
00483     typedef X domain_type;
00485     typedef X range_type;
00487     typedef typename X::field_type field_type;
00488 
00494     template<class L, class P>
00495     CGSolver (L& op, P& prec, double reduction, int maxit, int verbose) :
00496       ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00497     {
00498       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00499         "L and P must have the same category!");
00500       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
00501         "L must be sequential!");
00502     }
00508     template<class L, class S, class P>
00509     CGSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) :
00510       _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00511     {
00512       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00513         "L and P must have the same category!");
00514       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
00515         "L and S must have the same category!");
00516     }
00517 
00523     virtual void apply (X& x, X& b, InverseOperatorResult& res)
00524     {
00525       res.clear();                  // clear solver statistics
00526       Timer watch;                // start a timer
00527       _prec.pre(x,b);             // prepare preconditioner
00528       _op.applyscaleadd(-1,x,b);  // overwrite b with defect
00529 
00530       X p(x);              // the search direction
00531       X q(x);              // a temporary vector
00532 
00533       double def0 = _sp.norm(b);// compute norm
00534       if (def0<1E-30)    // convergence check
00535       {
00536         res.converged  = true;
00537         res.iterations = 0;               // fill statistics
00538         res.reduction = 0;
00539         res.conv_rate  = 0;
00540         res.elapsed=0;
00541         if (_verbose>0)                 // final print
00542           std::cout << "=== rate=" << res.conv_rate
00543                     << ", T=" << res.elapsed << ", TIT=" << res.elapsed
00544                     << ", IT=0" << std::endl;
00545         return;
00546       }
00547 
00548       if (_verbose>0)             // printing
00549       {
00550         std::cout << "=== CGSolver" << std::endl;
00551         if (_verbose>1) {
00552           this->printHeader(std::cout);
00553           this->printOutput(std::cout,0,def0);
00554         }
00555       }
00556 
00557       // some local variables
00558       double def=def0;   // loop variables
00559       field_type rho,rholast,lambda,alpha,beta;
00560 
00561       // determine initial search direction
00562       p = 0;                          // clear correction
00563       _prec.apply(p,b);               // apply preconditioner
00564       rholast = _sp.dot(p,b);         // orthogonalization
00565 
00566       // the loop
00567       int i=1;
00568       for ( ; i<=_maxit; i++ )
00569       {
00570         // minimize in given search direction p
00571         _op.apply(p,q);             // q=Ap
00572         alpha = _sp.dot(p,q);       // scalar product
00573         lambda = rholast/alpha;     // minimization
00574         x.axpy(lambda,p);           // update solution
00575         b.axpy(-lambda,q);          // update defect
00576 
00577         // convergence test
00578         double defnew=_sp.norm(b);// comp defect norm
00579 
00580         if (_verbose>1)             // print
00581           this->printOutput(std::cout,i,defnew,def);
00582 
00583         def = defnew;               // update norm
00584         if (def<def0*_reduction || def<1E-30)    // convergence check
00585         {
00586           res.converged  = true;
00587           break;
00588         }
00589 
00590         // determine new search direction
00591         q = 0;                      // clear correction
00592         _prec.apply(q,b);           // apply preconditioner
00593         rho = _sp.dot(q,b);         // orthogonalization
00594         beta = rho/rholast;         // scaling factor
00595         p *= beta;                  // scale old search direction
00596         p += q;                     // orthogonalization with correction
00597         rholast = rho;              // remember rho for recurrence
00598       }
00599 
00600       if (_verbose==1)                // printing for non verbose
00601         this->printOutput(std::cout,i,def);
00602 
00603       _prec.post(x);                  // postprocess preconditioner
00604       res.iterations = i;               // fill statistics
00605       res.reduction = def/def0;
00606       res.conv_rate  = pow(res.reduction,1.0/i);
00607       res.elapsed = watch.elapsed();
00608 
00609       if (_verbose>0)                 // final print
00610       {
00611         std::cout << "=== rate=" << res.conv_rate
00612                   << ", T=" << res.elapsed
00613                   << ", TIT=" << res.elapsed/i
00614                   << ", IT=" << i << std::endl;
00615       }
00616     }
00617 
00623     virtual void apply (X& x, X& b, double reduction,
00624       InverseOperatorResult& res)
00625     {
00626       std::swap(_reduction,reduction);
00627       (*this).apply(x,b,res);
00628       std::swap(_reduction,reduction);
00629     }
00630 
00631   private:
00632     SeqScalarProduct<X> ssp;
00633     LinearOperator<X,X>& _op;
00634     Preconditioner<X,X>& _prec;
00635     ScalarProduct<X>& _sp;
00636     double _reduction;
00637     int _maxit;
00638     int _verbose;
00639   };
00640 
00641 
00642   // Ronald Kriemanns BiCG-STAB implementation from Sumo
00644   template<class X>
00645   class BiCGSTABSolver : public InverseOperator<X,X> {
00646   public:
00648     typedef X domain_type;
00650     typedef X range_type;
00652     typedef typename X::field_type field_type;
00653 
00659     template<class L, class P>
00660     BiCGSTABSolver (L& op, P& prec,
00661       double reduction, int maxit, int verbose) :
00662       ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00663     {
00664       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), "L and P must be of the same category!");
00665       dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), "L must be sequential!");
00666     }
00672     template<class L, class S, class P>
00673     BiCGSTABSolver (L& op, S& sp, P& prec,
00674       double reduction, int maxit, int verbose) :
00675       _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00676     {
00677       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00678         "L and P must have the same category!");
00679       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
00680         "L and S must have the same category!");
00681     }
00682 
00688     virtual void apply (X& x, X& b, InverseOperatorResult& res)
00689     {
00690       const double EPSILON=1e-80;
00691 
00692       double               it;
00693       field_type           rho, rho_new, alpha, beta, h, omega;
00694       field_type           norm, norm_old, norm_0;
00695 
00696       //
00697       // get vectors and matrix
00698       //
00699       X& r=b;
00700       X p(x);
00701       X v(x);
00702       X t(x);
00703       X y(x);
00704       X rt(x);
00705 
00706       //
00707       // begin iteration
00708       //
00709 
00710       // r = r - Ax; rt = r
00711       res.clear();                // clear solver statistics
00712       Timer watch;                // start a timer
00713       _prec.pre(x,r);             // prepare preconditioner
00714       _op.applyscaleadd(-1,x,r);  // overwrite b with defect
00715 
00716       rt=r;
00717 
00718       norm = norm_old = norm_0 = _sp.norm(r);
00719 
00720       p=0;
00721       v=0;
00722 
00723       rho   = 1;
00724       alpha = 1;
00725       omega = 1;
00726 
00727       if (_verbose>0)             // printing
00728       {
00729         std::cout << "=== BiCGSTABSolver" << std::endl;
00730         if (_verbose>1)
00731         {
00732           this->printHeader(std::cout);
00733           this->printOutput(std::cout,0,norm_0);
00734           //std::cout << " Iter       Defect         Rate" << std::endl;
00735           //std::cout << "    0" << std::setw(14) << norm_0 << std::endl;
00736         }
00737       }
00738 
00739       if ( norm < (_reduction * norm_0)  || norm<1E-30)
00740       {
00741         res.converged = 1;
00742         _prec.post(x);                  // postprocess preconditioner
00743         res.iterations = 0;             // fill statistics
00744         res.reduction = 0;
00745         res.conv_rate  = 0;
00746         res.elapsed = watch.elapsed();
00747         return;
00748       }
00749 
00750       //
00751       // iteration
00752       //
00753 
00754       for (it = 0.5; it < _maxit; it+=.5)
00755       {
00756         //
00757         // preprocess, set vecsizes etc.
00758         //
00759 
00760         // rho_new = < rt , r >
00761         rho_new = _sp.dot(rt,r);
00762 
00763         // look if breakdown occured
00764         if (std::abs(rho) <= EPSILON)
00765           DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
00766             << rho << " <= EPSILON " << EPSILON
00767             << " after " << it << " iterations");
00768         if (std::abs(omega) <= EPSILON)
00769           DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
00770             << omega << " <= EPSILON " << EPSILON
00771             << " after " << it << " iterations");
00772 
00773 
00774         if (it<1)
00775           p = r;
00776         else
00777         {
00778           beta = ( rho_new / rho ) * ( alpha / omega );
00779           p.axpy(-omega,v); // p = r + beta (p - omega*v)
00780           p *= beta;
00781           p += r;
00782         }
00783 
00784         // y = W^-1 * p
00785         y = 0;
00786         _prec.apply(y,p);           // apply preconditioner
00787 
00788         // v = A * y
00789         _op.apply(y,v);
00790 
00791         // alpha = rho_new / < rt, v >
00792         h = _sp.dot(rt,v);
00793 
00794         if ( std::abs(h) < EPSILON )
00795           DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
00796 
00797         alpha = rho_new / h;
00798 
00799         // apply first correction to x
00800         // x <- x + alpha y
00801         x.axpy(alpha,y);
00802 
00803         // r = r - alpha*v
00804         r.axpy(-alpha,v);
00805 
00806         //
00807         // test stop criteria
00808         //
00809 
00810         norm = _sp.norm(r);
00811 
00812         if (_verbose>1) // print
00813         {
00814           this->printOutput(std::cout,it,norm,norm_old);
00815         }
00816 
00817         if ( norm < (_reduction * norm_0) )
00818         {
00819           res.converged = 1;
00820           break;
00821         }
00822         it+=.5;
00823 
00824         norm_old = norm;
00825 
00826         // y = W^-1 * r
00827         y = 0;
00828         _prec.apply(y,r);
00829 
00830         // t = A * y
00831         _op.apply(y,t);
00832 
00833         // omega = < t, r > / < t, t >
00834         omega = _sp.dot(t,r)/_sp.dot(t,t);
00835 
00836         // apply second correction to x
00837         // x <- x + omega y
00838         x.axpy(omega,y);
00839 
00840         // r = s - omega*t (remember : r = s)
00841         r.axpy(-omega,t);
00842 
00843         rho = rho_new;
00844 
00845         //
00846         // test stop criteria
00847         //
00848 
00849         norm = _sp.norm(r);
00850 
00851         if (_verbose > 1)             // print
00852         {
00853           this->printOutput(std::cout,it,norm,norm_old);
00854         }
00855 
00856         if ( norm < (_reduction * norm_0)  || norm<1E-30)
00857         {
00858           res.converged = 1;
00859           break;
00860         }
00861 
00862         norm_old = norm;
00863       } // end for
00864 
00865       if (_verbose==1)                // printing for non verbose
00866         this->printOutput(std::cout,it,norm);
00867 
00868       _prec.post(x);                  // postprocess preconditioner
00869       res.iterations = static_cast<int>(std::ceil(it));              // fill statistics
00870       res.reduction = norm/norm_0;
00871       res.conv_rate  = pow(res.reduction,1.0/it);
00872       res.elapsed = watch.elapsed();
00873       if (_verbose>0)                 // final print
00874         std::cout << "=== rate=" << res.conv_rate
00875                   << ", T=" << res.elapsed
00876                   << ", TIT=" << res.elapsed/it
00877                   << ", IT=" << it << std::endl;
00878     }
00879 
00885     virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
00886     {
00887       std::swap(_reduction,reduction);
00888       (*this).apply(x,b,res);
00889       std::swap(_reduction,reduction);
00890     }
00891 
00892   private:
00893     SeqScalarProduct<X> ssp;
00894     LinearOperator<X,X>& _op;
00895     Preconditioner<X,X>& _prec;
00896     ScalarProduct<X>& _sp;
00897     double _reduction;
00898     int _maxit;
00899     int _verbose;
00900   };
00901 
00908   template<class X>
00909   class MINRESSolver : public InverseOperator<X,X> {
00910   public:
00912     typedef X domain_type;
00914     typedef X range_type;
00916     typedef typename X::field_type field_type;
00918     typedef typename FieldTraits<field_type>::real_type real_type;
00919 
00925     template<class L, class P>
00926     MINRESSolver (L& op, P& prec, double reduction, int maxit, int verbose) :
00927       ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00928     {
00929       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00930         "L and P must have the same category!");
00931       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
00932         "L must be sequential!");
00933     }
00939     template<class L, class S, class P>
00940     MINRESSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) :
00941       _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
00942     {
00943       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
00944         "L and P must have the same category!");
00945       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
00946         "L and S must have the same category!");
00947     }
00948 
00954     virtual void apply (X& x, X& b, InverseOperatorResult& res)
00955     {
00956       res.clear();                // clear solver statistics
00957       Timer watch;                // start a timer
00958       _prec.pre(x,b);             // prepare preconditioner
00959       _op.applyscaleadd(-1,x,b);  // overwrite b with defect/residual
00960 
00961       real_type def0 = _sp.norm(b);   // compute residual norm
00962 
00963       if (def0<1E-30)    // convergence check
00964       {
00965         res.converged  = true;
00966         res.iterations = 0;               // fill statistics
00967         res.reduction = 0;
00968         res.conv_rate  = 0;
00969         res.elapsed=0;
00970         if (_verbose>0)                 // final print
00971           std::cout << "=== rate=" << res.conv_rate << ", T=" << res.elapsed << ", TIT=" << res.elapsed << ", IT=0" << std::endl;
00972         return;
00973       }
00974 
00975       if (_verbose>0)             // printing
00976       {
00977         std::cout << "=== MINRESSolver" << std::endl;
00978         if (_verbose>1) {
00979           this->printHeader(std::cout);
00980           this->printOutput(std::cout,0,def0);
00981         }
00982       }
00983 
00984       // some local variables
00985       real_type def=def0;                 // the defect/residual norm
00986       field_type alpha,                   // recurrence coefficients as computed in the Lanczos alg making up the matrix T
00987         c[2]={0.0, 0.0},         // diagonal entry of Givens rotation
00988         s[2]={0.0, 0.0};         // off-diagonal entries of Givens rotation
00989         real_type beta;
00990 
00991         field_type T[3]={0.0, 0.0, 0.0};    // recurrence coefficients (column k of Matrix T)
00992 
00993         X   z(b.size()),    // some temporary vectors
00994           dummy(b.size());
00995 
00996         field_type xi[2]={1.0, 0.0};
00997 
00998         // initialize
00999         z = 0.0;                // clear correction
01000 
01001         _prec.apply(z,b);       // apply preconditioner z=M^-1*b
01002 
01003         beta = std::sqrt(std::abs(_sp.dot(z,b)));
01004         real_type beta0 = beta;
01005 
01006         X p[3];     // the search directions
01007         X q[3];     // Orthonormal basis vectors (in unpreconditioned case)
01008 
01009         q[0].resize(b.size());
01010         q[1].resize(b.size());
01011         q[2].resize(b.size());
01012         q[0] = 0.0;
01013         q[1] = b;
01014         q[1] /= beta;
01015         q[2] = 0.0;
01016 
01017         p[0].resize(b.size());
01018         p[1].resize(b.size());
01019         p[2].resize(b.size());
01020         p[0] = 0.0;
01021         p[1] = 0.0;
01022         p[2] = 0.0;
01023 
01024 
01025         z /= beta;      // this is w_current
01026 
01027         // the loop
01028         int i=1;
01029         for ( ; i<=_maxit; i++)
01030         {
01031           dummy = z; // remember z_old for the computation of the search direction p in the next iteration
01032 
01033           int i1 = i%3,
01034             i0 = (i1+2)%3,
01035             i2 = (i1+1)%3;
01036 
01037           // Symmetrically Preconditioned Lanczos (Greenbaum p.121)
01038           _op.apply(z,q[i2]);             // q[i2] = Az
01039           q[i2].axpy(-beta, q[i0]);
01040           alpha = _sp.dot(q[i2],z);
01041           q[i2].axpy(-alpha, q[i1]);
01042 
01043           z=0.0;
01044           _prec.apply(z,q[i2]);
01045 
01046           beta = std::sqrt(std::abs(_sp.dot(q[i2],z)));
01047 
01048           q[i2] /= beta;
01049           z /= beta;
01050 
01051           // QR Factorization of recurrence coefficient matrix
01052           // apply previous Givens rotations to last column of T
01053           T[1] = T[2];
01054           if (i>2)
01055           {
01056             T[0] = s[i%2]*T[1];
01057             T[1] = c[i%2]*T[1];
01058           }
01059           if (i>1)
01060           {
01061             T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
01062             T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
01063           }
01064           else
01065             T[2] = alpha;
01066 
01067           // recompute c, s -> current Givens rotation \TODO use BLAS-routine drotg instead for greater robustness
01068 //          cblas_drotg (a, b, c, s);
01069           c[i%2] = 1.0/std::sqrt(T[2]*T[2] + beta*beta);
01070           s[i%2] = beta*c[i%2];
01071           c[i%2] *= T[2];
01072 
01073           // apply current Givens rotation to T eliminating the last entry...
01074           T[2] = c[i%2]*T[2] + s[i%2]*beta;
01075 
01076           // ...and to xi, the right hand side of the least squares problem min_y||beta*xi-T*y||
01077           xi[i%2] = -s[i%2]*xi[(i+1)%2];
01078           xi[(i+1)%2] *= c[i%2];
01079 
01080           // compute correction direction
01081           p[i2] = dummy;
01082           p[i2].axpy(-T[1],p[i1]);
01083           p[i2].axpy(-T[0],p[i0]);
01084           p[i2] /= T[2];
01085 
01086           // apply correction/update solution
01087           x.axpy(beta0*xi[(i+1)%2], p[i2]);
01088 
01089           // remember beta_old
01090           T[2] = beta;
01091 
01092           // update residual - not necessary if in the preconditioned case we are content with the residual norm of the
01093           // preconditioned system as convergence test
01094 //          _op.apply(p[i2],dummy);
01095 //          b.axpy(-beta0*xi[(i+1)%2],dummy);
01096 
01097 //          convergence test
01098           real_type defnew = std::abs(beta0*xi[i%2]); // the last entry the QR-transformed least squares RHS is the new residual norm
01099 
01100           if (_verbose>1)             // print
01101             this->printOutput(std::cout,i,defnew,def);
01102 
01103           def = defnew;               // update norm
01104           if (def<def0*_reduction || def<1E-30 || i==_maxit)    // convergence check
01105           {
01106             res.converged  = true;
01107             break;
01108           }
01109         }
01110 
01111         if (_verbose==1)                // printing for non verbose
01112           this->printOutput(std::cout,i,def);
01113 
01114         _prec.post(x);                  // postprocess preconditioner
01115         res.iterations = i;               // fill statistics
01116         res.reduction = def/def0;
01117         res.conv_rate  = pow(res.reduction,1.0/i);
01118         res.elapsed = watch.elapsed();
01119 
01120         if (_verbose>0)                 // final print
01121         {
01122           std::cout << "=== rate=" << res.conv_rate
01123                     << ", T=" << res.elapsed
01124                     << ", TIT=" << res.elapsed/i
01125                     << ", IT=" << i << std::endl;
01126         }
01127 
01128     }
01129 
01135     virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
01136     {
01137       std::swap(_reduction,reduction);
01138       (*this).apply(x,b,res);
01139       std::swap(_reduction,reduction);
01140     }
01141 
01142   private:
01143     SeqScalarProduct<X> ssp;
01144     LinearOperator<X,X>& _op;
01145     Preconditioner<X,X>& _prec;
01146     ScalarProduct<X>& _sp;
01147     double _reduction;
01148     int _maxit;
01149     int _verbose;
01150   };
01151 
01163   template<class X, class Y=X, class F = Y>
01164   class RestartedGMResSolver : public InverseOperator<X,Y>
01165   {
01166   public:
01168     typedef X domain_type;
01170     typedef Y range_type;
01172     typedef typename X::field_type field_type;
01174     typedef F basis_type;
01175 
01183     template<class L, class P>
01184     RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
01185       _A_(op), _M(prec),
01186       ssp(), _sp(ssp), _restart(restart),
01187       _reduction(reduction), _maxit(maxit), _verbose(verbose),
01188       _recalc_defect(recalc_defect)
01189     {
01190       dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
01191         "P and L must be the same category!");
01192       dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
01193         "L must be sequential!");
01194     }
01195 
01203     template<class L, class S, class P>
01204     RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
01205       _A_(op), _M(prec),
01206       _sp(sp), _restart(restart),
01207       _reduction(reduction), _maxit(maxit), _verbose(verbose),
01208       _recalc_defect(recalc_defect)
01209     {
01210       dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
01211         "P and L must have the same category!");
01212       dune_static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
01213         "P and S must have the same category!");
01214     }
01215 
01217     virtual void apply (X& x, X& b, InverseOperatorResult& res)
01218     {
01219       apply(x,b,_reduction,res);
01220     };
01221 
01227     virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
01228     {
01229       int m = _restart;
01230       field_type norm;
01231       field_type norm_old = 0.0;
01232       field_type norm_0;
01233       field_type beta;
01234       int i, j = 1, k;
01235       std::vector<field_type> s(m+1), cs(m), sn(m);
01236       // helper vector
01237       X w(b);
01238       std::vector< std::vector<field_type> > H(m+1,s);
01239       std::vector<F> v(m+1,b);
01240 
01241       // start timer
01242       Timer watch;                // start a timer
01243 
01244       // clear solver statistics
01245       res.clear();
01246       _M.pre(x,b);
01247       if (_recalc_defect)
01248       {
01249         // norm_0 = norm(M^-1 b)
01250         w = 0.0; _M.apply(w,b); // w = M^-1 b
01251         norm_0 = _sp.norm(w);
01252         // r = _M.solve(b - A * x);
01253         w = b;
01254         _A_.applyscaleadd(-1,x, /* => */ w); // w = b - Ax;
01255         v[0] = 0.0; _M.apply(v[0],w); // r = M^-1 w
01256         beta = _sp.norm(v[0]);
01257       }
01258       else
01259       {
01260         // norm_0 = norm(M^-1 b)
01261         w = 0.0; _M.apply(w,b); // w = M^-1 b
01262         // r = _M.solve(b - A * x);
01263         _A_.applyscaleadd(-1,x, /* => */ b); // b = b - Ax;
01264         norm_0 = _sp.norm(b);
01265         v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
01266         beta = _sp.norm(v[0]);
01267       }
01268 
01269       // avoid division by zero
01270       if (norm_0 == 0.0)
01271         norm_0 = 1.0;
01272       norm = norm_old = _sp.norm(v[0]);
01273 
01274       // print header
01275       if (_verbose > 0)
01276       {
01277         std::cout << "=== RestartedGMResSolver" << std::endl;
01278         if (_verbose > 1)
01279         {
01280           this->printHeader(std::cout);
01281           this->printOutput(std::cout,0,norm_0);
01282           this->printOutput(std::cout,0,norm, norm_0);
01283         }
01284       }
01285 
01286       // check convergence
01287       if (norm <= reduction * norm_0) {
01288         _M.post(x);                  // postprocess preconditioner
01289         res.converged  = true;
01290         if (_verbose > 0)                 // final print
01291           print_result(res);
01292         return;
01293       }
01294 
01295       while (j <= _maxit && res.converged != true) {
01296         v[0] *= (1.0 / beta);
01297         for (i=1; i<=m; i++) s[i] = 0.0;
01298         s[0] = beta;
01299 
01300         for (i = 0; i < m && j <= _maxit && res.converged != true; i++, j++) {
01301           w = 0.0;
01302           v[i+1] = 0.0; // use v[i+1] as temporary vector
01303           _A_.apply(v[i], /* => */ v[i+1]);
01304           _M.apply(w, v[i+1]);
01305           for (k = 0; k <= i; k++) {
01306             H[k][i] = _sp.dot(w, v[k]);
01307             // w -= H[k][i] * v[k];
01308             w.axpy(-H[k][i], v[k]);
01309           }
01310           H[i+1][i] = _sp.norm(w);
01311           if (H[i+1][i] == 0.0)
01312             DUNE_THROW(ISTLError,"breakdown in GMRes - |w| "
01313               << w << " == 0.0 after " << j << " iterations");
01314           // v[i+1] = w * (1.0 / H[i+1][i]);
01315           v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]);
01316 
01317           for (k = 0; k < i; k++)
01318             applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]);
01319 
01320           generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
01321           applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
01322           applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]);
01323 
01324           norm = std::abs(s[i+1]);
01325 
01326           if (_verbose > 1)             // print
01327           {
01328             this->printOutput(std::cout,j,norm,norm_old);
01329           }
01330 
01331           norm_old = norm;
01332 
01333           if (norm < reduction * norm_0) {
01334             res.converged = true;
01335           }
01336         }
01337 
01338         if (_recalc_defect)
01339         {
01340           // update x
01341           update(x, i - 1, H, s, v);
01342 
01343           // update residuum
01344           // r = M^-1 (b - A * x);
01345           w = b; _A_.applyscaleadd(-1,x, /* => */ w);
01346           _M.apply(v[0], w);
01347           beta = _sp.norm(v[0]);
01348           norm = beta;
01349         }
01350         else
01351         {
01352           // calc update vector
01353           w = 0;
01354           update(w, i - 1, H, s, v);
01355 
01356           // update x
01357           x += w;
01358 
01359           // r = M^-1 (b - A * x);
01360           // update defect
01361           _A_.applyscaleadd(-1,w, /* => */ b);
01362           // r = M^-1 (b - A * x);
01363           v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
01364           beta = _sp.norm(v[0]);
01365           norm = beta;
01366 
01367           res.converged = false;
01368         }
01369 
01370         if (_verbose > 1)             // print
01371         {
01372           this->printOutput(std::cout,j,norm,norm_old);
01373         }
01374 
01375         norm_old = norm;
01376 
01377         if (norm < reduction * norm_0) {
01378           // fill statistics
01379           res.converged = true;
01380         }
01381 
01382         if (res.converged != true && _verbose > 0)
01383           std::cout << "=== GMRes::restart\n";
01384       }
01385 
01386       _M.post(x);                  // postprocess preconditioner
01387 
01388       res.iterations = j;
01389       res.reduction = norm / norm_0;
01390       res.conv_rate  = pow(res.reduction,1.0/j);
01391       res.elapsed = watch.elapsed();
01392 
01393       if (_verbose>0)
01394         print_result(res);
01395     }
01396   private:
01397 
01398     void
01399     print_result (const InverseOperatorResult & res) const
01400     {
01401       int j = res.iterations>0?res.iterations:1;
01402       std::cout << "=== rate=" << res.conv_rate
01403                 << ", T=" << res.elapsed
01404                 << ", TIT=" << res.elapsed/j
01405                 << ", IT=" << res.iterations
01406                 << std::endl;
01407     }
01408 
01409     static void
01410     update(X &x, int k,
01411       std::vector< std::vector<field_type> > & h,
01412       std::vector<field_type> & s, std::vector<F> v)
01413     {
01414       std::vector<field_type> y(s);
01415 
01416       // Backsolve:
01417       for (int i = k; i >= 0; i--) {
01418         y[i] /= h[i][i];
01419         for (int j = i - 1; j >= 0; j--)
01420           y[j] -= h[j][i] * y[i];
01421       }
01422 
01423       for (int j = 0; j <= k; j++)
01424         // x += v[j] * y[j];
01425         x.axpy(y[j],v[j]);
01426     }
01427 
01428     void
01429     generatePlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
01430     {
01431       if (dy == 0.0) {
01432         cs = 1.0;
01433         sn = 0.0;
01434       } else if (std::abs(dy) > std::abs(dx)) {
01435         field_type temp = dx / dy;
01436         sn = 1.0 / std::sqrt( 1.0 + temp*temp );
01437         cs = temp * sn;
01438       } else {
01439         field_type temp = dy / dx;
01440         cs = 1.0 / std::sqrt( 1.0 + temp*temp );
01441         sn = temp * cs;
01442       }
01443     }
01444 
01445 
01446     void
01447     applyPlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
01448     {
01449       field_type temp  =  cs * dx + sn * dy;
01450       dy = -sn * dx + cs * dy;
01451       dx = temp;
01452     }
01453 
01454     LinearOperator<X,X>& _A_;
01455     Preconditioner<X,X>& _M;
01456     SeqScalarProduct<X> ssp;
01457     ScalarProduct<X>& _sp;
01458     int _restart;
01459     double _reduction;
01460     int _maxit;
01461     int _verbose;
01462     bool _recalc_defect;
01463   };
01464 
01467 } // end namespace
01468 
01469 #endif