Reference documentation for deal.II version GIT 3d869ba6cd 2022-05-16 18:15:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_bicgstab.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 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>
25 
26 #include <deal.II/lac/solver.h>
28 
29 #include <cmath>
30 
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  {
147  explicit AdditionalData(
148  const bool exact_residual = true,
149  const double breakdown =
153  {}
161  double breakdown;
162  };
163 
169  const AdditionalData & data = AdditionalData());
170 
176  const AdditionalData &data = AdditionalData());
177 
181  virtual ~SolverBicgstab() override = default;
182 
186  template <typename MatrixType, typename PreconditionerType>
187  void
188  solve(const MatrixType & A,
189  VectorType & x,
190  const VectorType & b,
191  const PreconditionerType &preconditioner);
192 
193 protected:
197  VectorType *Vx;
198 
203 
208 
213 
218 
223 
228 
233 
237  const VectorType *Vb;
238 
242  template <typename MatrixType>
243  double
244  criterion(const MatrixType &A, const VectorType &x, const VectorType &b);
245 
251  virtual void
252  print_vectors(const unsigned int step,
253  const VectorType & x,
254  const VectorType & r,
255  const VectorType & d) const;
256 
261 
262 private:
268  {
269  bool breakdown;
271  unsigned int last_step;
273 
276  const unsigned int last_step,
277  const double last_residual);
278  };
279 
284  template <typename MatrixType, typename PreconditionerType>
286  iterate(const MatrixType &A, const PreconditionerType &preconditioner);
287 };
288 
290 /*-------------------------Inline functions -------------------------------*/
291 
292 #ifndef DOXYGEN
293 
294 
295 template <typename VectorType>
297  const bool breakdown,
298  const SolverControl::State state,
299  const unsigned int last_step,
300  const double last_residual)
301  : breakdown(breakdown)
302  , state(state)
303  , last_step(last_step)
304  , last_residual(last_residual)
305 {}
306 
307 
308 
309 template <typename VectorType>
312  const AdditionalData & data)
313  : SolverBase<VectorType>(cn, mem)
314  , Vx(nullptr)
315  , Vb(nullptr)
316  , additional_data(data)
317 {}
318 
319 
320 
321 template <typename VectorType>
323  const AdditionalData &data)
324  : SolverBase<VectorType>(cn)
325  , Vx(nullptr)
326  , Vb(nullptr)
327  , additional_data(data)
328 {}
329 
330 
331 
332 template <typename VectorType>
333 template <typename MatrixType>
334 double
335 SolverBicgstab<VectorType>::criterion(const MatrixType &A,
336  const VectorType &x,
337  const VectorType &b)
338 {
339  A.vmult(*Vt, x);
340  Vt->add(-1., b);
341  res = Vt->l2_norm();
342 
343  return res;
344 }
345 
346 
347 
348 template <typename VectorType>
349 void
351  const VectorType &,
352  const VectorType &,
353  const VectorType &) const
354 {}
355 
356 
357 
358 template <typename VectorType>
359 template <typename MatrixType, typename PreconditionerType>
361 SolverBicgstab<VectorType>::iterate(const MatrixType & A,
362  const PreconditionerType &preconditioner)
363 {
364  A.vmult(*Vr, *Vx);
365  Vr->sadd(-1., 1., *Vb);
366  res = Vr->l2_norm();
367 
368  SolverControl::State state = this->iteration_status(step, res, *Vx);
369  if (state == SolverControl::State::success)
370  return IterationResult(false, state, step, res);
371 
372  alpha = omega = rho = 1.;
373 
374  VectorType &r = *Vr;
375  VectorType &rbar = *Vrbar;
376  VectorType &p = *Vp;
377  VectorType &y = *Vy;
378  VectorType &z = *Vz;
379  VectorType &t = *Vt;
380  VectorType &v = *Vv;
381 
382  rbar = r;
383  bool startup = true;
384 
385  do
386  {
387  ++step;
388 
389  rhobar = r * rbar;
390  if (std::fabs(rhobar) < additional_data.breakdown)
391  {
392  return IterationResult(true, state, step, res);
393  }
394  beta = rhobar * alpha / (rho * omega);
395  rho = rhobar;
396  if (startup == true)
397  {
398  p = r;
399  startup = false;
400  }
401  else
402  {
403  p.sadd(beta, 1., r);
404  p.add(-beta * omega, v);
405  }
406 
407  preconditioner.vmult(y, p);
408  A.vmult(v, y);
409  rhobar = rbar * v;
410  if (std::fabs(rhobar) < additional_data.breakdown)
411  {
412  return IterationResult(true, state, step, res);
413  }
414 
415  alpha = rho / rhobar;
416 
417  res = std::sqrt(r.add_and_dot(-alpha, v, r));
418 
419  // check for early success, see the lac/bicgstab_early testcase as to
420  // why this is necessary
421  //
422  // note: the vector *Vx we pass to the iteration_status signal here is
423  // only the current approximation, not the one we will return with, which
424  // will be x=*Vx + alpha*y
425  if (this->iteration_status(step, res, *Vx) == SolverControl::success)
426  {
427  Vx->add(alpha, y);
428  print_vectors(step, *Vx, r, y);
429  return IterationResult(false, SolverControl::success, step, res);
430  }
431 
432  preconditioner.vmult(z, r);
433  A.vmult(t, z);
434  rhobar = t * r;
435  auto t_squared = t * t;
436  if (t_squared < additional_data.breakdown)
437  {
438  return IterationResult(true, state, step, res);
439  }
440  omega = rhobar / (t * t);
441  Vx->add(alpha, y, omega, z);
442 
443  if (additional_data.exact_residual)
444  {
445  r.add(-omega, t);
446  res = criterion(A, *Vx, *Vb);
447  }
448  else
449  res = std::sqrt(r.add_and_dot(-omega, t, r));
450 
451  state = this->iteration_status(step, res, *Vx);
452  print_vectors(step, *Vx, r, y);
453  }
454  while (state == SolverControl::iterate);
455 
456  return IterationResult(false, state, step, res);
457 }
458 
459 
460 
461 template <typename VectorType>
462 template <typename MatrixType, typename PreconditionerType>
463 void
464 SolverBicgstab<VectorType>::solve(const MatrixType & A,
465  VectorType & x,
466  const VectorType & b,
467  const PreconditionerType &preconditioner)
468 {
469  LogStream::Prefix prefix("Bicgstab");
470 
471  // Allocate temporary memory.
472  Vr = typename VectorMemory<VectorType>::Pointer(this->memory);
473  Vrbar = typename VectorMemory<VectorType>::Pointer(this->memory);
474  Vp = typename VectorMemory<VectorType>::Pointer(this->memory);
475  Vy = typename VectorMemory<VectorType>::Pointer(this->memory);
476  Vz = typename VectorMemory<VectorType>::Pointer(this->memory);
477  Vt = typename VectorMemory<VectorType>::Pointer(this->memory);
478  Vv = typename VectorMemory<VectorType>::Pointer(this->memory);
479 
480  Vr->reinit(x, true);
481  Vrbar->reinit(x, true);
482  Vp->reinit(x, true);
483  Vy->reinit(x, true);
484  Vz->reinit(x, true);
485  Vt->reinit(x, true);
486  Vv->reinit(x, true);
487 
488  Vx = &x;
489  Vb = &b;
490 
491  step = 0;
492 
493  IterationResult state(false, SolverControl::failure, 0, 0);
494  do
495  {
496  state = iterate(A, preconditioner);
497  }
498  while (state.state == SolverControl::iterate);
499 
500 
501  // Release the temporary memory again.
502  Vr.reset();
503  Vrbar.reset();
504  Vp.reset();
505  Vy.reset();
506  Vz.reset();
507  Vt.reset();
508  Vv.reset();
509 
510  // In case of failure: throw exception
511  AssertThrow(state.state == SolverControl::success,
512  SolverControl::NoConvergence(state.last_step,
513  state.last_residual));
514  // Otherwise exit as normal
515 }
516 
517 #endif // DOXYGEN
518 
520 
521 #endif
const VectorType * Vb
VectorMemory< VectorType >::Pointer Vr
VectorMemory< VectorType >::Pointer Vp
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
virtual ~SolverBicgstab() override=default
AdditionalData additional_data
VectorMemory< VectorType >::Pointer Vv
VectorMemory< VectorType >::Pointer Vy
double criterion(const MatrixType &A, const VectorType &x, const VectorType &b)
VectorMemory< VectorType >::Pointer Vrbar
VectorMemory< VectorType >::Pointer Vt
SolverBicgstab(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
VectorType * Vx
IterationResult iterate(const MatrixType &A, const PreconditionerType &preconditioner)
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
VectorMemory< VectorType >::Pointer Vz
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
@ failure
Stop iteration, goal not reached.
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Expression fabs(const Expression &x)
static const char A
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const bool exact_residual=true, const double breakdown=std::numeric_limits< typename VectorType::value_type >::min())
IterationResult(const bool breakdown, const SolverControl::State state, const unsigned int last_step, const double last_residual)