Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_bicgstab.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_solver_bicgstab_h
17 #define dealii_solver_bicgstab_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/logstream.h>
23 #include <deal.II/base/signaling_nan.h>
24 #include <deal.II/base/subscriptor.h>
25 
26 #include <deal.II/lac/solver.h>
27 #include <deal.II/lac/solver_control.h>
28 
29 #include <cmath>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
35 
36 namespace internal
37 {
43  {
44  protected:
48  double alpha;
52  double beta;
56  double omega;
60  double rho;
64  double rhobar;
65 
69  unsigned int step;
70 
74  double res;
75 
81  };
82 } // namespace internal
83 
123 template <typename VectorType = Vector<double>>
124 class SolverBicgstab : public SolverBase<VectorType>,
126 {
127 public:
139  {
146  explicit AdditionalData(const bool exact_residual = true,
147  const double breakdown = 1.e-10)
148  : exact_residual(exact_residual)
149  , breakdown(breakdown)
150  {}
158  double breakdown;
159  };
160 
166  const AdditionalData & data = AdditionalData());
167 
173  const AdditionalData &data = AdditionalData());
174 
178  virtual ~SolverBicgstab() override = default;
179 
183  template <typename MatrixType, typename PreconditionerType>
184  void
185  solve(const MatrixType & A,
186  VectorType & x,
187  const VectorType & b,
188  const PreconditionerType &preconditioner);
189 
190 protected:
194  VectorType *Vx;
195 
200 
205 
210 
215 
220 
225 
230 
234  const VectorType *Vb;
235 
239  template <typename MatrixType>
240  double
241  criterion(const MatrixType &A, const VectorType &x, const VectorType &b);
242 
248  virtual void
249  print_vectors(const unsigned int step,
250  const VectorType & x,
251  const VectorType & r,
252  const VectorType & d) const;
253 
258 
259 private:
263  template <typename MatrixType>
265  start(const MatrixType &A);
266 
272  {
273  bool breakdown;
274  SolverControl::State state;
275  unsigned int last_step;
276  double last_residual;
277 
278  IterationResult(const bool breakdown,
279  const SolverControl::State state,
280  const unsigned int last_step,
281  const double last_residual);
282  };
283 
288  template <typename MatrixType, typename PreconditionerType>
290  iterate(const MatrixType &A, const PreconditionerType &preconditioner);
291 };
292 
294 /*-------------------------Inline functions -------------------------------*/
295 
296 #ifndef DOXYGEN
297 
298 
299 template <typename VectorType>
301  const bool breakdown,
302  const SolverControl::State state,
303  const unsigned int last_step,
304  const double last_residual)
305  : breakdown(breakdown)
306  , state(state)
307  , last_step(last_step)
308  , last_residual(last_residual)
309 {}
310 
311 
312 
313 template <typename VectorType>
316  const AdditionalData & data)
317  : SolverBase<VectorType>(cn, mem)
318  , Vx(nullptr)
319  , Vb(nullptr)
320  , additional_data(data)
321 {}
322 
323 
324 
325 template <typename VectorType>
327  const AdditionalData &data)
329  , Vx(nullptr)
330  , Vb(nullptr)
331  , additional_data(data)
332 {}
333 
334 
335 
336 template <typename VectorType>
337 template <typename MatrixType>
338 double
339 SolverBicgstab<VectorType>::criterion(const MatrixType &A,
340  const VectorType &x,
341  const VectorType &b)
342 {
343  A.vmult(*Vt, x);
344  Vt->add(-1., b);
345  res = Vt->l2_norm();
346 
347  return res;
348 }
349 
350 
351 
352 template <typename VectorType>
353 template <typename MatrixType>
355 SolverBicgstab<VectorType>::start(const MatrixType &A)
356 {
357  A.vmult(*Vr, *Vx);
358  Vr->sadd(-1., 1., *Vb);
359  res = Vr->l2_norm();
360 
361  return this->iteration_status(step, res, *Vx);
362 }
363 
364 
365 
366 template <typename VectorType>
367 void
369  const VectorType &,
370  const VectorType &,
371  const VectorType &) const
372 {}
373 
374 
375 
376 template <typename VectorType>
377 template <typename MatrixType, typename PreconditionerType>
379 SolverBicgstab<VectorType>::iterate(const MatrixType & A,
380  const PreconditionerType &preconditioner)
381 {
382  // TODO:[GK] Implement "use the length of the computed orthogonal residual" in
383  // the BiCGStab method.
385  alpha = omega = rho = 1.;
386 
387  VectorType &r = *Vr;
388  VectorType &rbar = *Vrbar;
389  VectorType &p = *Vp;
390  VectorType &y = *Vy;
391  VectorType &z = *Vz;
392  VectorType &t = *Vt;
393  VectorType &v = *Vv;
394 
395  rbar = r;
396  bool startup = true;
397 
398  do
399  {
400  ++step;
401 
402  rhobar = r * rbar;
403  beta = rhobar * alpha / (rho * omega);
404  rho = rhobar;
405  if (startup == true)
406  {
407  p = r;
408  startup = false;
409  }
410  else
411  {
412  p.sadd(beta, 1., r);
413  p.add(-beta * omega, v);
414  }
415 
416  preconditioner.vmult(y, p);
417  A.vmult(v, y);
418  rhobar = rbar * v;
419 
420  alpha = rho / rhobar;
421 
422  // TODO:[?] Find better breakdown criterion
423 
424  if (std::fabs(alpha) > 1.e10)
425  return IterationResult(true, state, step, res);
426 
427  res = std::sqrt(r.add_and_dot(-alpha, v, r));
428 
429  // check for early success, see the lac/bicgstab_early testcase as to
430  // why this is necessary
431  //
432  // note: the vector *Vx we pass to the iteration_status signal here is
433  // only the current approximation, not the one we will return with, which
434  // will be x=*Vx + alpha*y
435  if (this->iteration_status(step, res, *Vx) == SolverControl::success)
436  {
437  Vx->add(alpha, y);
438  print_vectors(step, *Vx, r, y);
440  }
441 
442  preconditioner.vmult(z, r);
443  A.vmult(t, z);
444  rhobar = t * r;
445  omega = rhobar / (t * t);
446  Vx->add(alpha, y, omega, z);
447 
448  if (additional_data.exact_residual)
449  {
450  r.add(-omega, t);
451  res = criterion(A, *Vx, *Vb);
452  }
453  else
454  res = std::sqrt(r.add_and_dot(-omega, t, r));
455 
456  state = this->iteration_status(step, res, *Vx);
457  print_vectors(step, *Vx, r, y);
458  }
459  while (state == SolverControl::iterate);
460  return IterationResult(false, state, step, res);
461 }
462 
463 
464 template <typename VectorType>
465 template <typename MatrixType, typename PreconditionerType>
466 void
467 SolverBicgstab<VectorType>::solve(const MatrixType & A,
468  VectorType & x,
469  const VectorType & b,
470  const PreconditionerType &preconditioner)
471 {
472  LogStream::Prefix prefix("Bicgstab");
473  Vr = typename VectorMemory<VectorType>::Pointer(this->memory);
474  Vrbar = typename VectorMemory<VectorType>::Pointer(this->memory);
475  Vp = typename VectorMemory<VectorType>::Pointer(this->memory);
476  Vy = typename VectorMemory<VectorType>::Pointer(this->memory);
477  Vz = typename VectorMemory<VectorType>::Pointer(this->memory);
478  Vt = typename VectorMemory<VectorType>::Pointer(this->memory);
479  Vv = typename VectorMemory<VectorType>::Pointer(this->memory);
480 
481  Vr->reinit(x, true);
482  Vrbar->reinit(x, true);
483  Vp->reinit(x, true);
484  Vy->reinit(x, true);
485  Vz->reinit(x, true);
486  Vt->reinit(x, true);
487  Vv->reinit(x, true);
488 
489  Vx = &x;
490  Vb = &b;
491 
492  step = 0;
493 
494  IterationResult state(false, SolverControl::failure, 0, 0);
495 
496  // iterate while the inner iteration returns a breakdown
497  do
498  {
499  if (step != 0)
500  deallog << "Restart step " << step << std::endl;
501  if (start(A) == SolverControl::success)
502  {
503  state.state = SolverControl::success;
504  break;
505  }
506  state = iterate(A, preconditioner);
507  ++step;
508  }
509  while (state.breakdown == true);
510 
511  // in case of failure: throw exception
512  AssertThrow(state.state == SolverControl::success,
513  SolverControl::NoConvergence(state.last_step,
514  state.last_residual));
515  // otherwise exit as normal
516 }
517 
518 #endif // DOXYGEN
519 
520 DEAL_II_NAMESPACE_CLOSE
521 
522 #endif
Stop iteration, goal not reached.
SolverBicgstab(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Continue iteration.
VectorMemory< VectorType >::Pointer Vrbar
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
VectorMemory< VectorType >::Pointer Vv
VectorMemory< VectorType >::Pointer Vp
VectorMemory< VectorType >::Pointer Vr
Stop iteration, goal reached.
AdditionalData additional_data
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
VectorMemory< VectorType >::Pointer Vt
AdditionalData(const bool exact_residual=true, const double breakdown=1.e-10)
const VectorType * Vb
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
IterationResult iterate(const MatrixType &A, const PreconditionerType &preconditioner)
VectorMemory< VectorType >::Pointer Vy
VectorMemory< VectorType >::Pointer Vz
SolverControl::State start(const MatrixType &A)
VectorType * Vx
double criterion(const MatrixType &A, const VectorType &x, const VectorType &b)