15#ifndef dealii_solver_bicgstab_h
16#define dealii_solver_bicgstab_h
77template <
typename VectorType = Vector<
double>>
103 std::numeric_limits<typename VectorType::value_type>::min())
139 template <
typename MatrixType,
typename PreconditionerType>
144 const PreconditionerType &preconditioner);
150 template <
typename MatrixType>
166 const VectorType &d)
const;
195 template <
typename MatrixType,
typename PreconditionerType>
200 const PreconditionerType &preconditioner,
201 const unsigned int step);
211template <
typename VectorType>
213 const bool breakdown,
215 const unsigned int last_step,
216 const double last_residual)
217 : breakdown(breakdown)
219 , last_step(last_step)
220 , last_residual(last_residual)
225template <
typename VectorType>
228 const AdditionalData &data)
230 , additional_data(data)
235template <
typename VectorType>
237 const AdditionalData &data)
239 , additional_data(data)
244template <
typename VectorType>
245template <
typename MatrixType>
253 return std::sqrt(t.add_and_dot(-1.0, b, t));
258template <
typename VectorType>
263 const VectorType &)
const
268template <
typename VectorType>
269template <
typename MatrixType,
typename PreconditionerType>
274 const PreconditionerType &preconditioner,
275 const unsigned int last_step)
288 VectorType &rbar = *Vrbar;
296 rbar.reinit(x,
true);
303 using value_type =
typename VectorType::value_type;
308 value_type res = r.l2_norm();
310 unsigned int step = last_step;
314 return IterationResult(
false, state, step, res);
318 value_type alpha = 1.;
320 value_type omega = 1.;
326 const value_type rhobar = (step == 1 + last_step) ? res * res : r * rbar;
328 if (std::fabs(rhobar) < additional_data.breakdown)
330 return IterationResult(
true, state, step, res);
333 const value_type beta = rhobar * alpha / (rho * omega);
335 if (step == last_step + 1)
342 p.add(-beta * omega, v);
345 preconditioner.vmult(y, p);
347 const value_type rbar_dot_v = rbar * v;
348 if (std::fabs(rbar_dot_v) < additional_data.breakdown)
350 return IterationResult(
true, state, step, res);
353 alpha = rho / rbar_dot_v;
355 res =
std::sqrt(real_type(r.add_and_dot(-alpha, v, r)));
366 print_vectors(step, x, r, y);
370 preconditioner.vmult(z, r);
372 const value_type t_dot_r = t * r;
373 const real_type t_squared = t * t;
374 if (t_squared < additional_data.breakdown)
376 return IterationResult(
true, state, step, res);
378 omega = t_dot_r / t_squared;
379 x.add(alpha, y, omega, z);
381 if (additional_data.exact_residual)
384 res = criterion(A, x, b, t);
387 res =
std::sqrt(real_type(r.add_and_dot(-omega, t, r)));
389 state = this->iteration_status(step, res, x);
390 print_vectors(step, x, r, y);
394 return IterationResult(
false, state, step, res);
399template <
typename VectorType>
400template <
typename MatrixType,
typename PreconditionerType>
405 const PreconditionerType &preconditioner)
412 state = iterate(A, x, b, preconditioner, state.last_step);
419 state.last_residual));
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
#define DEAL_II_NAMESPACE_CLOSE
#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())
SolverControl::State state
IterationResult(const bool breakdown, const SolverControl::State state, const unsigned int last_step, const double last_residual)