Reference documentation for deal.II version 9.4.1
\(\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\}}\)
Loading...
Searching...
No Matches
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
25
26#include <deal.II/lac/solver.h>
28
29#include <cmath>
30
32
35
36namespace 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
123template <typename VectorType = Vector<double>>
124class SolverBicgstab : public SolverBase<VectorType>,
126{
127public:
139 {
148 const bool exact_residual = true,
149 const double breakdown =
150 std::numeric_limits<typename VectorType::value_type>::min())
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
193protected:
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
262private:
268 {
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
295template <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
309template <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
321template <typename VectorType>
323 const AdditionalData &data)
324 : SolverBase<VectorType>(cn)
325 , Vx(nullptr)
326 , Vb(nullptr)
327 , additional_data(data)
328{}
329
330
331
332template <typename VectorType>
333template <typename MatrixType>
334double
335SolverBicgstab<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
348template <typename VectorType>
349void
351 const VectorType &,
352 const VectorType &,
353 const VectorType &) const
354{}
355
356
357
358template <typename VectorType>
359template <typename MatrixType, typename PreconditionerType>
361SolverBicgstab<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);
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
461template <typename VectorType>
462template <typename MatrixType, typename PreconditionerType>
463void
464SolverBicgstab<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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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)