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