00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
00355
00356 if (std::fabs(alpha) > 1.e10)
00357 return true;
00358
00359 s.equ(1., r, -alpha, v);
00360
00361
00362
00363
00364
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
00440
00441 if (this->control().last_check() != SolverControl::success)
00442 throw SolverControl::NoConvergence (this->control().last_step(),
00443 this->control().last_value());
00444
00445 }
00446
00447 #endif // DOXYGEN
00448
00449 DEAL_II_NAMESPACE_CLOSE
00450
00451 #endif
00452