Reference documentation for deal.II version GIT relicensing-1054-gedd46abbd6 2024-07-12 01:20: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\}}\)
Loading...
Searching...
No Matches
solver_bicgstab.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_bicgstab_h
16#define dealii_solver_bicgstab_h
17
18
19#include <deal.II/base/config.h>
20
25
26#include <deal.II/lac/solver.h>
28
29#include <cmath>
30#include <limits>
31
33
78template <typename VectorType = Vector<double>>
80class SolverBicgstab : public SolverBase<VectorType>
81{
82public:
94 {
103 const bool exact_residual = true,
104 const double breakdown =
105 std::numeric_limits<typename VectorType::value_type>::min())
106 : exact_residual(exact_residual)
107 , breakdown(breakdown)
108 {}
116 double breakdown;
117 };
118
124 const AdditionalData &data = AdditionalData());
125
131 const AdditionalData &data = AdditionalData());
132
136 virtual ~SolverBicgstab() override = default;
137
141 template <typename MatrixType, typename PreconditionerType>
145 void solve(const MatrixType &A,
146 VectorType &x,
147 const VectorType &b,
148 const PreconditionerType &preconditioner);
149
150protected:
154 template <typename MatrixType>
155 double
156 criterion(const MatrixType &A,
157 const VectorType &x,
158 const VectorType &b,
159 VectorType &t);
160
166 virtual void
167 print_vectors(const unsigned int step,
168 const VectorType &x,
169 const VectorType &r,
170 const VectorType &d) const;
171
175 AdditionalData additional_data;
176
177private:
183 {
186 unsigned int last_step;
188
189 IterationResult(const bool breakdown,
190 const SolverControl::State state,
191 const unsigned int last_step,
192 const double last_residual);
193 };
194
199 template <typename MatrixType, typename PreconditionerType>
201 iterate(const MatrixType &A,
202 VectorType &x,
203 const VectorType &b,
204 const PreconditionerType &preconditioner,
205 const unsigned int step);
206};
207
208
210/*-------------------------Inline functions -------------------------------*/
211
212#ifndef DOXYGEN
213
214
215template <typename VectorType>
218 const bool breakdown,
219 const SolverControl::State state,
220 const unsigned int last_step,
221 const double last_residual)
222 : breakdown(breakdown)
223 , state(state)
224 , last_step(last_step)
225 , last_residual(last_residual)
226{}
227
228
229
230template <typename VectorType>
234 const AdditionalData &data)
235 : SolverBase<VectorType>(cn, mem)
236 , additional_data(data)
237{}
238
239
240
241template <typename VectorType>
244 const AdditionalData &data)
246 , additional_data(data)
247{}
248
249
250
251template <typename VectorType>
253template <typename MatrixType>
254double SolverBicgstab<VectorType>::criterion(const MatrixType &A,
255 const VectorType &x,
256 const VectorType &b,
257 VectorType &t)
258{
259 A.vmult(t, x);
260 return std::sqrt(t.add_and_dot(-1.0, b, t));
261}
262
263
264
265template <typename VectorType>
267void SolverBicgstab<VectorType>::print_vectors(const unsigned int,
268 const VectorType &,
269 const VectorType &,
270 const VectorType &) const
271{}
272
273
274
275template <typename VectorType>
277template <typename MatrixType, typename PreconditionerType>
279 SolverBicgstab<VectorType>::iterate(const MatrixType &A,
280 VectorType &x,
281 const VectorType &b,
282 const PreconditionerType &preconditioner,
283 const unsigned int last_step)
284{
285 // Allocate temporary memory.
286 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
287 typename VectorMemory<VectorType>::Pointer Vrbar(this->memory);
288 typename VectorMemory<VectorType>::Pointer Vp(this->memory);
289 typename VectorMemory<VectorType>::Pointer Vy(this->memory);
290 typename VectorMemory<VectorType>::Pointer Vz(this->memory);
291 typename VectorMemory<VectorType>::Pointer Vt(this->memory);
292 typename VectorMemory<VectorType>::Pointer Vv(this->memory);
293
294 // Define a few aliases for simpler use of the vectors
295 VectorType &r = *Vr;
296 VectorType &rbar = *Vrbar;
297 VectorType &p = *Vp;
298 VectorType &y = *Vy;
299 VectorType &z = *Vz;
300 VectorType &t = *Vt;
301 VectorType &v = *Vv;
302
303 r.reinit(x, true);
304 rbar.reinit(x, true);
305 p.reinit(x, true);
306 y.reinit(x, true);
307 z.reinit(x, true);
308 t.reinit(x, true);
309 v.reinit(x, true);
310
311 using value_type = typename VectorType::value_type;
312 using real_type = typename numbers::NumberTraits<value_type>::real_type;
313
314 A.vmult(r, x);
315 r.sadd(-1., 1., b);
316 value_type res = r.l2_norm();
317
318 unsigned int step = last_step;
319
320 SolverControl::State state = this->iteration_status(step, res, x);
322 return IterationResult(false, state, step, res);
323
324 rbar = r;
325
326 value_type alpha = 1.;
327 value_type rho = 1.;
328 value_type omega = 1.;
329
330 do
331 {
332 ++step;
333
334 const value_type rhobar = (step == 1 + last_step) ? res * res : r * rbar;
335
336 if (std::fabs(rhobar) < additional_data.breakdown)
337 {
338 return IterationResult(true, state, step, res);
339 }
340
341 const value_type beta = rhobar * alpha / (rho * omega);
342 rho = rhobar;
343 if (step == last_step + 1)
344 {
345 p = r;
346 }
347 else
348 {
349 p.sadd(beta, 1., r);
350 p.add(-beta * omega, v);
351 }
352
353 preconditioner.vmult(y, p);
354 A.vmult(v, y);
355 const value_type rbar_dot_v = rbar * v;
356 if (std::fabs(rbar_dot_v) < additional_data.breakdown)
357 {
358 return IterationResult(true, state, step, res);
359 }
360
361 alpha = rho / rbar_dot_v;
362
363 res = std::sqrt(real_type(r.add_and_dot(-alpha, v, r)));
364
365 // check for early success, see the lac/bicgstab_early testcase as to
366 // why this is necessary
367 //
368 // note: the vector *Vx we pass to the iteration_status signal here is
369 // only the current approximation, not the one we will return with, which
370 // will be x=*Vx + alpha*y
371 if (this->iteration_status(step, res, x) == SolverControl::success)
372 {
373 x.add(alpha, y);
374 print_vectors(step, x, r, y);
375 return IterationResult(false, SolverControl::success, step, res);
376 }
377
378 preconditioner.vmult(z, r);
379 A.vmult(t, z);
380 const value_type t_dot_r = t * r;
381 const real_type t_squared = t * t;
382 if (t_squared < additional_data.breakdown)
383 {
384 return IterationResult(true, state, step, res);
385 }
386 omega = t_dot_r / t_squared;
387 x.add(alpha, y, omega, z);
388
389 if (additional_data.exact_residual)
390 {
391 r.add(-omega, t);
392 res = criterion(A, x, b, t);
393 }
394 else
395 res = std::sqrt(real_type(r.add_and_dot(-omega, t, r)));
396
397 state = this->iteration_status(step, res, x);
398 print_vectors(step, x, r, y);
399 }
400 while (state == SolverControl::iterate);
401
402 return IterationResult(false, state, step, res);
403}
404
405
406
407template <typename VectorType>
409template <typename MatrixType, typename PreconditionerType>
413void SolverBicgstab<VectorType>::solve(const MatrixType &A,
414 VectorType &x,
415 const VectorType &b,
416 const PreconditionerType &preconditioner)
417{
418 LogStream::Prefix prefix("Bicgstab");
419
420 IterationResult state(false, SolverControl::failure, 0, 0);
421 do
422 {
423 state = iterate(A, x, b, preconditioner, state.last_step);
424 }
425 while (state.state == SolverControl::iterate);
426
427 // In case of failure: throw exception
428 AssertThrow(state.state == SolverControl::success,
429 SolverControl::NoConvergence(state.last_step,
430 state.last_residual));
431 // Otherwise exit as normal
432}
433
434#endif // DOXYGEN
435
437
438#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
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:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define AssertThrow(cond, exc)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
::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)