Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
solver_bicgstab.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_bicgstab.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 1998 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__solver_bicgstab_h
18 #define __deal2__solver_bicgstab_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/logstream.h>
23 #include <deal.II/lac/solver.h>
24 #include <deal.II/lac/solver_control.h>
25 #include <cmath>
26 #include <deal.II/base/subscriptor.h>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
32 
66 template <class VECTOR = Vector<double> >
67 class SolverBicgstab : public Solver<VECTOR>
68 {
69 public:
81  {
88  AdditionalData(const bool exact_residual = true,
89  const double breakdown = 1.e-10) :
92  {}
100  double breakdown;
101  };
102 
108  const AdditionalData &data=AdditionalData());
109 
115  const AdditionalData &data=AdditionalData());
116 
120  virtual ~SolverBicgstab ();
121 
125  template<class MATRIX, class PRECONDITIONER>
126  void
127  solve (const MATRIX &A,
128  VECTOR &x,
129  const VECTOR &b,
130  const PRECONDITIONER &precondition);
131 
132 protected:
136  template <class MATRIX>
137  double criterion (const MATRIX &A, const VECTOR &x, const VECTOR &b);
138 
144  virtual void print_vectors(const unsigned int step,
145  const VECTOR &x,
146  const VECTOR &r,
147  const VECTOR &d) const;
148 
152  VECTOR *Vx;
156  VECTOR *Vr;
160  VECTOR *Vrbar;
164  VECTOR *Vp;
168  VECTOR *Vy;
172  VECTOR *Vz;
176  VECTOR *Vt;
180  VECTOR *Vv;
184  const VECTOR *Vb;
185 
189  double alpha;
193  double beta;
197  double omega;
201  double rho;
205  double rhobar;
206 
210  unsigned int step;
211 
215  double res;
216 
221 
222 private:
226  template <class MATRIX>
228 
232  template<class MATRIX, class PRECONDITIONER>
233  bool
234  iterate(const MATRIX &A, const PRECONDITIONER &precondition);
235 
236 };
237 
239 /*-------------------------Inline functions -------------------------------*/
240 
241 #ifndef DOXYGEN
242 
243 template<class VECTOR>
246  const AdditionalData &data)
247  :
248  Solver<VECTOR>(cn,mem),
249  additional_data(data)
250 {}
251 
252 
253 
254 template<class VECTOR>
256  const AdditionalData &data)
257  :
258  Solver<VECTOR>(cn),
259  additional_data(data)
260 {}
261 
262 
263 
264 template<class VECTOR>
266 {}
267 
268 
269 
270 template <class VECTOR>
271 template <class MATRIX>
272 double
273 SolverBicgstab<VECTOR>::criterion (const MATRIX &A, const VECTOR &x, const VECTOR &b)
274 {
275  A.vmult(*Vt, x);
276  Vt->add(-1.,b);
277  res = Vt->l2_norm();
278 
279  return res;
280 }
281 
282 
283 
284 template <class VECTOR >
285 template <class MATRIX>
288 {
289  A.vmult(*Vr, *Vx);
290  Vr->sadd(-1.,1.,*Vb);
291  res = Vr->l2_norm();
292 
293  return this->control().check(step, res);
294 }
295 
296 
297 
298 template<class VECTOR>
299 void
300 SolverBicgstab<VECTOR>::print_vectors(const unsigned int,
301  const VECTOR &,
302  const VECTOR &,
303  const VECTOR &) const
304 {}
305 
306 
307 
308 template<class VECTOR>
309 template<class MATRIX, class PRECONDITIONER>
310 bool
312  const PRECONDITIONER &precondition)
313 {
314 //TODO:[GK] Implement "use the length of the computed orthogonal residual" in the BiCGStab method.
316  alpha = omega = rho = 1.;
317 
318  VECTOR &r = *Vr;
319  VECTOR &rbar = *Vrbar;
320  VECTOR &p = *Vp;
321  VECTOR &y = *Vy;
322  VECTOR &z = *Vz;
323  VECTOR &t = *Vt;
324  VECTOR &v = *Vv;
325 
326  rbar = r;
327  bool startup = true;
328 
329  do
330  {
331  ++step;
332 
333  rhobar = r*rbar;
334  beta = rhobar * alpha / (rho * omega);
335  rho = rhobar;
336  if (startup == true)
337  {
338  p = r;
339  startup = false;
340  }
341  else
342  p.sadd(beta, 1., r, -beta*omega, v);
343 
344  precondition.vmult(y,p);
345  A.vmult(v,y);
346  rhobar = rbar * v;
347 
348  alpha = rho/rhobar;
349 
350 //TODO:[?] Find better breakdown criterion
351 
352  if (std::fabs(alpha) > 1.e10)
353  return true;
354 
355  r.add(-alpha, v);
356 
357  // check for early success, see the lac/bicgstab_early testcase as to
358  // why this is necessary
359  if (this->control().check(step, r.l2_norm()) == SolverControl::success)
360  {
361  Vx->add(alpha, y);
362  print_vectors(step, *Vx, r, y);
363  return false;
364  }
365 
366  precondition.vmult(z,r);
367  A.vmult(t,z);
368  rhobar = t*r;
369  omega = rhobar/(t*t);
370  Vx->add(alpha, y, omega, z);
371  r.add(-omega, t);
372 
373  if (additional_data.exact_residual)
374  res = criterion(A, *Vx, *Vb);
375  else
376  res = r.l2_norm();
377 
378  state = this->control().check(step, res);
379  print_vectors(step, *Vx, r, y);
380  }
381  while (state == SolverControl::iterate);
382  return false;
383 }
384 
385 
386 template<class VECTOR>
387 template<class MATRIX, class PRECONDITIONER>
388 void
390  VECTOR &x,
391  const VECTOR &b,
392  const PRECONDITIONER &precondition)
393 {
394  deallog.push("Bicgstab");
395  Vr = this->memory.alloc();
396  Vr->reinit(x, true);
397  Vrbar = this->memory.alloc();
398  Vrbar->reinit(x, true);
399  Vp = this->memory.alloc();
400  Vp->reinit(x, true);
401  Vy = this->memory.alloc();
402  Vy->reinit(x, true);
403  Vz = this->memory.alloc();
404  Vz->reinit(x, true);
405  Vt = this->memory.alloc();
406  Vt->reinit(x, true);
407  Vv = this->memory.alloc();
408  Vv->reinit(x, true);
409 
410  Vx = &x;
411  Vb = &b;
412 
413  step = 0;
414 
415  bool state;
416 
417  do
418  {
419  if (step != 0)
420  deallog << "Restart step " << step << std::endl;
421  if (start(A) == SolverControl::success)
422  break;
423  state = iterate(A, precondition);
424  }
425  while (state);
426 
427  this->memory.free(Vr);
428  this->memory.free(Vrbar);
429  this->memory.free(Vp);
430  this->memory.free(Vy);
431  this->memory.free(Vz);
432  this->memory.free(Vt);
433  this->memory.free(Vv);
434 
435  deallog.pop();
436 
437  // in case of failure: throw exception
438  if (this->control().last_check() != SolverControl::success)
439  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
440  this->control().last_value()));
441  // otherwise exit as normal
442 }
443 
444 #endif // DOXYGEN
445 
446 DEAL_II_NAMESPACE_CLOSE
447 
448 #endif
AdditionalData(const bool exact_residual=true, const double breakdown=1.e-10)
void pop()
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
SolverBicgstab(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
double criterion(const MATRIX &A, const VECTOR &x, const VECTOR &b)
Stop iteration, goal reached.
unsigned int step
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
Definition: solver.h:147
bool iterate(const MATRIX &A, const PRECONDITIONER &precondition)
void push(const std::string &text)
AdditionalData additional_data
const VECTOR * Vb
SolverControl::State start(const MATRIX &A)
virtual ~SolverBicgstab()