Reference documentation for deal.II version 9.5.0
\(\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#include <limits>
31
33
78template <typename VectorType = Vector<double>>
79class SolverBicgstab : public SolverBase<VectorType>
80{
81public:
93 {
102 const bool exact_residual = true,
103 const double breakdown =
104 std::numeric_limits<typename VectorType::value_type>::min())
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
147protected:
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
174private:
180 {
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
212template <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
226template <typename VectorType>
229 const AdditionalData & data)
230 : SolverBase<VectorType>(cn, mem)
231 , additional_data(data)
232{}
233
234
235
236template <typename VectorType>
238 const AdditionalData &data)
239 : SolverBase<VectorType>(cn)
240 , additional_data(data)
241{}
242
243
244
245template <typename VectorType>
246template <typename MatrixType>
247double
248SolverBicgstab<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
259template <typename VectorType>
260void
262 const VectorType &,
263 const VectorType &,
264 const VectorType &) const
265{}
266
267
268
269template <typename VectorType>
270template <typename MatrixType, typename PreconditionerType>
272SolverBicgstab<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);
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
400template <typename VectorType>
401template <typename MatrixType, typename PreconditionerType>
402void
403SolverBicgstab<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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define AssertThrow(cond, exc)
::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)