include/deal.II/lac/solver_bicgstab.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: solver_bicgstab.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__solver_bicgstab_h
00013 #define __deal2__solver_bicgstab_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/logstream.h>
00018 #include <deal.II/lac/solver.h>
00019 #include <deal.II/lac/solver_control.h>
00020 #include <cmath>
00021 #include <deal.II/base/subscriptor.h>
00022 
00023 DEAL_II_NAMESPACE_OPEN
00024 
00027 
00061 template <class VECTOR = Vector<double> >
00062 class SolverBicgstab : public Solver<VECTOR>
00063 {
00064   public:
00080     struct AdditionalData
00081     {
00089         AdditionalData(const bool   exact_residual = true,
00090                        const double breakdown      = 1.e-10) :
00091                         exact_residual(exact_residual),
00092                         breakdown(breakdown)
00093           {}
00097         bool exact_residual;
00101         double breakdown;
00102     };
00103 
00107     SolverBicgstab (SolverControl        &cn,
00108                     VectorMemory<VECTOR> &mem,
00109                     const AdditionalData &data=AdditionalData());
00110 
00116     SolverBicgstab (SolverControl        &cn,
00117                     const AdditionalData &data=AdditionalData());
00118 
00122     virtual ~SolverBicgstab ();
00123 
00127     template<class MATRIX, class PRECONDITIONER>
00128     void
00129     solve (const MATRIX &A,
00130            VECTOR       &x,
00131            const VECTOR &b,
00132            const PRECONDITIONER& precondition);
00133 
00134   protected:
00138     template <class MATRIX>
00139     double criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b);
00140 
00150     virtual void print_vectors(const unsigned int step,
00151                                const VECTOR& x,
00152                                const VECTOR& r,
00153                                const VECTOR& d) const;
00154 
00158     VECTOR *Vx;
00162     VECTOR *Vr;
00166     VECTOR *Vrbar;
00170     VECTOR *Vp;
00174     VECTOR *Vy;
00178     VECTOR *Vz;
00182     VECTOR *Vs;
00186     VECTOR *Vt;
00190     VECTOR *Vv;
00194     const VECTOR *Vb;
00195 
00199     double alpha;
00203     double beta;
00207     double omega;
00211     double rho;
00215     double rhobar;
00216 
00220     unsigned int step;
00221 
00225     double res;
00226 
00230     AdditionalData additional_data;
00231 
00232   private:
00236     template <class MATRIX>
00237     SolverControl::State start(const MATRIX& A);
00238 
00242     template<class MATRIX, class PRECONDITIONER>
00243     bool
00244     iterate(const MATRIX& A, const PRECONDITIONER& precondition);
00245 
00246 };
00247 
00249 /*-------------------------Inline functions -------------------------------*/
00250 
00251 #ifndef DOXYGEN
00252 
00253 template<class VECTOR>
00254 SolverBicgstab<VECTOR>::SolverBicgstab (SolverControl &cn,
00255                                         VectorMemory<VECTOR> &mem,
00256                                         const AdditionalData &data)
00257                 :
00258                 Solver<VECTOR>(cn,mem),
00259                 additional_data(data)
00260 {}
00261 
00262 
00263 
00264 template<class VECTOR>
00265 SolverBicgstab<VECTOR>::SolverBicgstab (SolverControl &cn,
00266                                         const AdditionalData &data)
00267                 :
00268                 Solver<VECTOR>(cn),
00269                 additional_data(data)
00270 {}
00271 
00272 
00273 
00274 template<class VECTOR>
00275 SolverBicgstab<VECTOR>::~SolverBicgstab ()
00276 {}
00277 
00278 
00279 
00280 template <class VECTOR>
00281 template <class MATRIX>
00282 double
00283 SolverBicgstab<VECTOR>::criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b)
00284 {
00285   A.vmult(*Vt, x);
00286   Vt->sadd(-1.,1.,b);
00287   res = Vt->l2_norm();
00288 
00289   return res;
00290 }
00291 
00292 
00293 
00294 template <class VECTOR >
00295 template <class MATRIX>
00296 SolverControl::State
00297 SolverBicgstab<VECTOR>::start(const MATRIX& A)
00298 {
00299   A.vmult(*Vr, *Vx);
00300   Vr->sadd(-1.,1.,*Vb);
00301   res = Vr->l2_norm();
00302 
00303   Vp->reinit(*Vx);
00304   Vv->reinit(*Vx);
00305   *Vrbar = *Vr;
00306   return this->control().check(step, res);
00307 }
00308 
00309 
00310 
00311 template<class VECTOR>
00312 void
00313 SolverBicgstab<VECTOR>::print_vectors(const unsigned int,
00314                                       const VECTOR&,
00315                                       const VECTOR&,
00316                                       const VECTOR&) const
00317 {}
00318 
00319 
00320 
00321 template<class VECTOR>
00322 template<class MATRIX, class PRECONDITIONER>
00323 bool
00324 SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
00325                                 const PRECONDITIONER& precondition)
00326 {
00327 //TODO:[GK] Implement "use the length of the computed orthogonal residual" in the BiCGStab method.
00328   SolverControl::State state = SolverControl::iterate;
00329   alpha = omega = rho = 1.;
00330 
00331   VECTOR& r = *Vr;
00332   VECTOR& rbar = *Vrbar;
00333   VECTOR& p = *Vp;
00334   VECTOR& y = *Vy;
00335   VECTOR& z = *Vz;
00336   VECTOR& s = *Vs;
00337   VECTOR& t = *Vt;
00338   VECTOR& v = *Vv;
00339 
00340   do
00341     {
00342       ++step;
00343 
00344       rhobar = r*rbar;
00345       beta   = rhobar * alpha / (rho * omega);
00346       rho    = rhobar;
00347       p.sadd(beta, 1., r, -beta*omega, v);
00348       precondition.vmult(y,p);
00349       A.vmult(v,y);
00350       rhobar = rbar * v;
00351 
00352       alpha = rho/rhobar;
00353 
00354 //TODO:[?] Find better breakdown criterion
00355 
00356       if (std::fabs(alpha) > 1.e10)
00357         return true;
00358 
00359       s.equ(1., r, -alpha, v);
00360 
00361                                        // check for early success, see
00362                                        // the lac/bicgstab_early
00363                                        // testcase as to why this is
00364                                        // necessary
00365       if (this->control().check(step, s.l2_norm()/Vb->l2_norm())
00366           == SolverControl::success)
00367         {
00368           Vx->add(alpha, y);
00369           print_vectors(step, *Vx, r, y);
00370           return false;
00371         }
00372 
00373       precondition.vmult(z,s);
00374       A.vmult(t,z);
00375       rhobar = t*s;
00376       omega = rhobar/(t*t);
00377       Vx->add(alpha, y, omega, z);
00378       r.equ(1., s, -omega, t);
00379 
00380       if (additional_data.exact_residual)
00381         res = criterion(A, *Vx, *Vb);
00382       else
00383         res = r.l2_norm();
00384 
00385       state = this->control().check(step, res);
00386       print_vectors(step, *Vx, r, y);
00387     }
00388   while (state == SolverControl::iterate);
00389   return false;
00390 }
00391 
00392 
00393 template<class VECTOR>
00394 template<class MATRIX, class PRECONDITIONER>
00395 void
00396 SolverBicgstab<VECTOR>::solve(const MATRIX &A,
00397                               VECTOR       &x,
00398                               const VECTOR &b,
00399                               const PRECONDITIONER& precondition)
00400 {
00401   deallog.push("Bicgstab");
00402   Vr    = this->memory.alloc(); Vr->reinit(x);
00403   Vrbar = this->memory.alloc(); Vrbar->reinit(x);
00404   Vp    = this->memory.alloc();
00405   Vy    = this->memory.alloc(); Vy->reinit(x);
00406   Vz    = this->memory.alloc(); Vz->reinit(x);
00407   Vs    = this->memory.alloc(); Vs->reinit(x);
00408   Vt    = this->memory.alloc(); Vt->reinit(x);
00409   Vv    = this->memory.alloc();
00410 
00411   Vx = &x;
00412   Vb = &b;
00413 
00414   step = 0;
00415 
00416   bool state;
00417 
00418   do
00419     {
00420       if (step != 0)
00421         deallog << "Restart step " << step << std::endl;
00422       if (start(A) == SolverControl::success)
00423         break;
00424       state = iterate(A, precondition);
00425     }
00426   while (state);
00427 
00428   this->memory.free(Vr);
00429   this->memory.free(Vrbar);
00430   this->memory.free(Vp);
00431   this->memory.free(Vy);
00432   this->memory.free(Vz);
00433   this->memory.free(Vs);
00434   this->memory.free(Vt);
00435   this->memory.free(Vv);
00436 
00437   deallog.pop();
00438 
00439                                    // in case of failure: throw
00440                                    // exception
00441   if (this->control().last_check() != SolverControl::success)
00442     throw SolverControl::NoConvergence (this->control().last_step(),
00443                                         this->control().last_value());
00444                                    // otherwise exit as normal
00445 }
00446 
00447 #endif // DOXYGEN
00448 
00449 DEAL_II_NAMESPACE_CLOSE
00450 
00451 #endif
00452 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3