16#ifndef dealii_solver_bicgstab_h
17#define dealii_solver_bicgstab_h
123template <
typename VectorType = Vector<
double>>
150 std::numeric_limits<typename VectorType::value_type>::min())
186 template <
typename MatrixType,
typename PreconditionerType>
190 const VectorType & b,
191 const PreconditionerType &preconditioner);
237 const VectorType *
Vb;
242 template <
typename MatrixType>
244 criterion(
const MatrixType &A,
const VectorType &x,
const VectorType &b);
253 const VectorType & x,
254 const VectorType & r,
255 const VectorType & d)
const;
284 template <
typename MatrixType,
typename PreconditionerType>
286 iterate(
const MatrixType &A,
const PreconditionerType &preconditioner);
295template <
typename VectorType>
297 const bool breakdown,
299 const unsigned int last_step,
300 const double last_residual)
301 : breakdown(breakdown)
303 , last_step(last_step)
304 , last_residual(last_residual)
309template <
typename VectorType>
312 const AdditionalData & data)
316 , additional_data(data)
321template <
typename VectorType>
323 const AdditionalData &data)
327 , additional_data(data)
332template <
typename VectorType>
333template <
typename MatrixType>
348template <
typename VectorType>
353 const VectorType &)
const
358template <
typename VectorType>
359template <
typename MatrixType,
typename PreconditionerType>
362 const PreconditionerType &preconditioner)
365 Vr->sadd(-1., 1., *Vb);
370 return IterationResult(
false, state, step, res);
372 alpha = omega = rho = 1.;
375 VectorType &rbar = *Vrbar;
390 if (std::fabs(rhobar) < additional_data.breakdown)
392 return IterationResult(
true, state, step, res);
394 beta = rhobar * alpha / (rho * omega);
404 p.add(-beta * omega, v);
407 preconditioner.vmult(y, p);
410 if (std::fabs(rhobar) < additional_data.breakdown)
412 return IterationResult(
true, state, step, res);
415 alpha = rho / rhobar;
417 res =
std::sqrt(r.add_and_dot(-alpha, v, r));
428 print_vectors(step, *Vx, r, y);
432 preconditioner.vmult(z, r);
435 auto t_squared = t * t;
436 if (t_squared < additional_data.breakdown)
438 return IterationResult(
true, state, step, res);
440 omega = rhobar / (t * t);
441 Vx->add(alpha, y, omega, z);
443 if (additional_data.exact_residual)
446 res = criterion(A, *Vx, *Vb);
449 res =
std::sqrt(r.add_and_dot(-omega, t, r));
451 state = this->iteration_status(step, res, *Vx);
452 print_vectors(step, *Vx, r, y);
456 return IterationResult(
false, state, step, res);
461template <
typename VectorType>
462template <
typename MatrixType,
typename PreconditionerType>
466 const VectorType & b,
467 const PreconditionerType &preconditioner)
481 Vrbar->reinit(x,
true);
496 state = iterate(A, preconditioner);
513 state.last_residual));
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())
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
#define DEAL_II_NAMESPACE_CLOSE
#define AssertThrow(cond, exc)
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())
SolverControl::State state
IterationResult(const bool breakdown, const SolverControl::State state, const unsigned int last_step, const double last_residual)