Reference documentation for deal.II version GIT f3d24b0aad 2022-12-01 05:00:02+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 - 2022 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 #include <limits>
31 
33 
78 template <typename VectorType = Vector<double>>
79 class SolverBicgstab : public SolverBase<VectorType>
80 {
81 public:
93  {
101  explicit AdditionalData(
102  const bool exact_residual = true,
103  const double breakdown =
107  {}
115  double breakdown;
116  };
117 
123  const AdditionalData & data = AdditionalData());
124 
130  const AdditionalData &data = AdditionalData());
131 
135  virtual ~SolverBicgstab() override = default;
136 
140  template <typename MatrixType, typename PreconditionerType>
141  void
142  solve(const MatrixType & A,
143  VectorType & x,
144  const VectorType & b,
145  const PreconditionerType &preconditioner);
146 
147 protected:
151  template <typename MatrixType>
152  double
153  criterion(const MatrixType &A,
154  const VectorType &x,
155  const VectorType &b,
156  VectorType & t);
157 
163  virtual void
164  print_vectors(const unsigned int step,
165  const VectorType & x,
166  const VectorType & r,
167  const VectorType & d) const;
168 
173 
174 private:
180  {
181  bool breakdown;
183  unsigned int last_step;
185 
188  const unsigned int last_step,
189  const double last_residual);
190  };
191 
196  template <typename MatrixType, typename PreconditionerType>
198  iterate(const MatrixType & A,
199  VectorType & x,
200  const VectorType & b,
201  const PreconditionerType &preconditioner,
202  const unsigned int step);
203 };
204 
205 
207 /*-------------------------Inline functions -------------------------------*/
208 
209 #ifndef DOXYGEN
210 
211 
212 template <typename VectorType>
214  const bool breakdown,
215  const SolverControl::State state,
216  const unsigned int last_step,
217  const double last_residual)
218  : breakdown(breakdown)
219  , state(state)
220  , last_step(last_step)
221  , last_residual(last_residual)
222 {}
223 
224 
225 
226 template <typename VectorType>
229  const AdditionalData & data)
230  : SolverBase<VectorType>(cn, mem)
231  , additional_data(data)
232 {}
233 
234 
235 
236 template <typename VectorType>
238  const AdditionalData &data)
239  : SolverBase<VectorType>(cn)
240  , additional_data(data)
241 {}
242 
243 
244 
245 template <typename VectorType>
246 template <typename MatrixType>
247 double
248 SolverBicgstab<VectorType>::criterion(const MatrixType &A,
249  const VectorType &x,
250  const VectorType &b,
251  VectorType & t)
252 {
253  A.vmult(t, x);
254  return std::sqrt(t.add_and_dot(-1.0, b, t));
255 }
256 
257 
258 
259 template <typename VectorType>
260 void
262  const VectorType &,
263  const VectorType &,
264  const VectorType &) const
265 {}
266 
267 
268 
269 template <typename VectorType>
270 template <typename MatrixType, typename PreconditionerType>
272 SolverBicgstab<VectorType>::iterate(const MatrixType & A,
273  VectorType & x,
274  const VectorType & b,
275  const PreconditionerType &preconditioner,
276  const unsigned int last_step)
277 {
278  // Allocate temporary memory.
279  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
280  typename VectorMemory<VectorType>::Pointer Vrbar(this->memory);
281  typename VectorMemory<VectorType>::Pointer Vp(this->memory);
282  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
283  typename VectorMemory<VectorType>::Pointer Vz(this->memory);
284  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
285  typename VectorMemory<VectorType>::Pointer Vv(this->memory);
286 
287  // Define a few aliases for simpler use of the vectors
288  VectorType &r = *Vr;
289  VectorType &rbar = *Vrbar;
290  VectorType &p = *Vp;
291  VectorType &y = *Vy;
292  VectorType &z = *Vz;
293  VectorType &t = *Vt;
294  VectorType &v = *Vv;
295 
296  r.reinit(x, true);
297  rbar.reinit(x, true);
298  p.reinit(x, true);
299  y.reinit(x, true);
300  z.reinit(x, true);
301  t.reinit(x, true);
302  v.reinit(x, true);
303 
304  using value_type = typename VectorType::value_type;
305  using real_type = typename numbers::NumberTraits<value_type>::real_type;
306 
307  A.vmult(r, x);
308  r.sadd(-1., 1., b);
309  value_type res = r.l2_norm();
310 
311  unsigned int step = last_step;
312 
313  SolverControl::State state = this->iteration_status(step, res, x);
314  if (state == SolverControl::State::success)
315  return IterationResult(false, state, step, res);
316 
317  rbar = r;
318 
319  value_type alpha = 1.;
320  value_type rho = 1.;
321  value_type omega = 1.;
322 
323  do
324  {
325  ++step;
326 
327  const value_type rhobar = (step == 1 + last_step) ? res * res : r * rbar;
328 
329  if (std::fabs(rhobar) < additional_data.breakdown)
330  {
331  return IterationResult(true, state, step, res);
332  }
333 
334  const value_type beta = rhobar * alpha / (rho * omega);
335  rho = rhobar;
336  if (step == last_step + 1)
337  {
338  p = r;
339  }
340  else
341  {
342  p.sadd(beta, 1., r);
343  p.add(-beta * omega, v);
344  }
345 
346  preconditioner.vmult(y, p);
347  A.vmult(v, y);
348  const value_type rbar_dot_v = rbar * v;
349  if (std::fabs(rbar_dot_v) < additional_data.breakdown)
350  {
351  return IterationResult(true, state, step, res);
352  }
353 
354  alpha = rho / rbar_dot_v;
355 
356  res = std::sqrt(real_type(r.add_and_dot(-alpha, v, r)));
357 
358  // check for early success, see the lac/bicgstab_early testcase as to
359  // why this is necessary
360  //
361  // note: the vector *Vx we pass to the iteration_status signal here is
362  // only the current approximation, not the one we will return with, which
363  // will be x=*Vx + alpha*y
364  if (this->iteration_status(step, res, x) == SolverControl::success)
365  {
366  x.add(alpha, y);
367  print_vectors(step, x, r, y);
368  return IterationResult(false, SolverControl::success, step, res);
369  }
370 
371  preconditioner.vmult(z, r);
372  A.vmult(t, z);
373  const value_type t_dot_r = t * r;
374  const real_type t_squared = t * t;
375  if (t_squared < additional_data.breakdown)
376  {
377  return IterationResult(true, state, step, res);
378  }
379  omega = t_dot_r / t_squared;
380  x.add(alpha, y, omega, z);
381 
382  if (additional_data.exact_residual)
383  {
384  r.add(-omega, t);
385  res = criterion(A, x, b, t);
386  }
387  else
388  res = std::sqrt(real_type(r.add_and_dot(-omega, t, r)));
389 
390  state = this->iteration_status(step, res, x);
391  print_vectors(step, x, r, y);
392  }
393  while (state == SolverControl::iterate);
394 
395  return IterationResult(false, state, step, res);
396 }
397 
398 
399 
400 template <typename VectorType>
401 template <typename MatrixType, typename PreconditionerType>
402 void
403 SolverBicgstab<VectorType>::solve(const MatrixType & A,
404  VectorType & x,
405  const VectorType & b,
406  const PreconditionerType &preconditioner)
407 {
408  LogStream::Prefix prefix("Bicgstab");
409 
410  IterationResult state(false, SolverControl::failure, 0, 0);
411  do
412  {
413  state = iterate(A, x, b, preconditioner, state.last_step);
414  }
415  while (state.state == SolverControl::iterate);
416 
417  // In case of failure: throw exception
418  AssertThrow(state.state == SolverControl::success,
419  SolverControl::NoConvergence(state.last_step,
420  state.last_residual));
421  // Otherwise exit as normal
422 }
423 
424 #endif // DOXYGEN
425 
427 
428 #endif
double criterion(const MatrixType &A, const VectorType &x, const VectorType &b, VectorType &t)
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
SolverBicgstab(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
IterationResult iterate(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner, const unsigned int step)
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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)