solver_bicgstab.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1998 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_solver_bicgstab_h
#define dealii_solver_bicgstab_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/template_constraints.h>
 
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
 
#include <cmath>
#include <limits>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename VectorType = Vector<double>>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
class SolverBicgstab : public SolverBase<VectorType>
{
public:
  struct AdditionalData
  {
    explicit AdditionalData(
      const bool   exact_residual = true,
      const double breakdown =
        std::numeric_limits<typename VectorType::value_type>::min())
      : exact_residual(exact_residual)
      , breakdown(breakdown)
    {}
    bool exact_residual;
    double breakdown;
  };
 
  SolverBicgstab(SolverControl            &cn,
                 VectorMemory<VectorType> &mem,
                 const AdditionalData     &data = AdditionalData());
 
  SolverBicgstab(SolverControl        &cn,
                 const AdditionalData &data = AdditionalData());
 
  virtual ~SolverBicgstab() override = default;
 
  template <typename MatrixType, typename PreconditionerType>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_linear_operator_on<MatrixType, VectorType> &&
     concepts::is_linear_operator_on<PreconditionerType, VectorType>))
  void solve(const MatrixType         &A,
             VectorType               &x,
             const VectorType         &b,
             const PreconditionerType &preconditioner);
 
protected:
  template <typename MatrixType>
  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;
 
  AdditionalData additional_data;
 
private:
  struct IterationResult
  {
    bool                 breakdown;
    SolverControl::State state;
    unsigned int         last_step;
    double               last_residual;
 
    IterationResult(const bool                 breakdown,
                    const SolverControl::State state,
                    const unsigned int         last_step,
                    const double               last_residual);
  };
 
  template <typename MatrixType, typename PreconditionerType>
  IterationResult
  iterate(const MatrixType         &A,
          VectorType               &x,
          const VectorType         &b,
          const PreconditionerType &preconditioner,
          const unsigned int        step);
};
 
 
/*-------------------------Inline functions -------------------------------*/
 
#ifndef DOXYGEN
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverBicgstab<VectorType>::IterationResult::IterationResult(
  const bool                 breakdown,
  const SolverControl::State state,
  const unsigned int         last_step,
  const double               last_residual)
  : breakdown(breakdown)
  , state(state)
  , last_step(last_step)
  , last_residual(last_residual)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverBicgstab<VectorType>::SolverBicgstab(SolverControl            &cn,
                                           VectorMemory<VectorType> &mem,
                                           const AdditionalData     &data)
  : SolverBase<VectorType>(cn, mem)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverBicgstab<VectorType>::SolverBicgstab(SolverControl        &cn,
                                           const AdditionalData &data)
  : SolverBase<VectorType>(cn)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType>
double SolverBicgstab<VectorType>::criterion(const MatrixType &A,
                                             const VectorType &x,
                                             const VectorType &b,
                                             VectorType       &t)
{
  A.vmult(t, x);
  return std::sqrt(t.add_and_dot(-1.0, b, t));
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
void SolverBicgstab<VectorType>::print_vectors(const unsigned int,
                                               const VectorType &,
                                               const VectorType &,
                                               const VectorType &) const
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType, typename PreconditionerType>
typename SolverBicgstab<VectorType>::IterationResult
  SolverBicgstab<VectorType>::iterate(const MatrixType         &A,
                                      VectorType               &x,
                                      const VectorType         &b,
                                      const PreconditionerType &preconditioner,
                                      const unsigned int        last_step)
{
  // Allocate temporary memory.
  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
  typename VectorMemory<VectorType>::Pointer Vrbar(this->memory);
  typename VectorMemory<VectorType>::Pointer Vp(this->memory);
  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
  typename VectorMemory<VectorType>::Pointer Vz(this->memory);
  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
  typename VectorMemory<VectorType>::Pointer Vv(this->memory);
 
  // Define a few aliases for simpler use of the vectors
  VectorType &r    = *Vr;
  VectorType &rbar = *Vrbar;
  VectorType &p    = *Vp;
  VectorType &y    = *Vy;
  VectorType &z    = *Vz;
  VectorType &t    = *Vt;
  VectorType &v    = *Vv;
 
  r.reinit(x, true);
  rbar.reinit(x, true);
  p.reinit(x, true);
  y.reinit(x, true);
  z.reinit(x, true);
  t.reinit(x, true);
  v.reinit(x, true);
 
  using value_type = typename VectorType::value_type;
  using real_type  = typename numbers::NumberTraits<value_type>::real_type;
 
  A.vmult(r, x);
  r.sadd(-1., 1., b);
  value_type res = r.l2_norm();
 
  unsigned int step = last_step;
 
  SolverControl::State state = this->iteration_status(step, res, x);
  if (state == SolverControl::State::success)
    return IterationResult(false, state, step, res);
 
  rbar = r;
 
  value_type alpha = 1.;
  value_type rho   = 1.;
  value_type omega = 1.;
 
  do
    {
      ++step;
 
      const value_type rhobar = (step == 1 + last_step) ? res * res : r * rbar;
 
      if (std::fabs(rhobar) < additional_data.breakdown)
        {
          return IterationResult(true, state, step, res);
        }
 
      const value_type beta = rhobar * alpha / (rho * omega);
      rho                   = rhobar;
      if (step == last_step + 1)
        {
          p = r;
        }
      else
        {
          p.sadd(beta, 1., r);
          p.add(-beta * omega, v);
        }
 
      preconditioner.vmult(y, p);
      A.vmult(v, y);
      const value_type rbar_dot_v = rbar * v;
      if (std::fabs(rbar_dot_v) < additional_data.breakdown)
        {
          return IterationResult(true, state, step, res);
        }
 
      alpha = rho / rbar_dot_v;
 
      res = std::sqrt(real_type(r.add_and_dot(-alpha, v, r)));
 
      // check for early success, see the lac/bicgstab_early testcase as to
      // why this is necessary
      //
      // note: the vector *Vx we pass to the iteration_status signal here is
      // only the current approximation, not the one we will return with, which
      // will be x=*Vx + alpha*y
      if (this->iteration_status(step, res, x) == SolverControl::success)
        {
          x.add(alpha, y);
          print_vectors(step, x, r, y);
          return IterationResult(false, SolverControl::success, step, res);
        }
 
      preconditioner.vmult(z, r);
      A.vmult(t, z);
      const value_type t_dot_r   = t * r;
      const real_type  t_squared = t * t;
      if (t_squared < additional_data.breakdown)
        {
          return IterationResult(true, state, step, res);
        }
      omega = t_dot_r / t_squared;
      x.add(alpha, y, omega, z);
 
      if (additional_data.exact_residual)
        {
          r.add(-omega, t);
          res = criterion(A, x, b, t);
        }
      else
        res = std::sqrt(real_type(r.add_and_dot(-omega, t, r)));
 
      state = this->iteration_status(step, res, x);
      print_vectors(step, x, r, y);
    }
  while (state == SolverControl::iterate);
 
  return IterationResult(false, state, step, res);
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType, typename PreconditionerType>
DEAL_II_CXX20_REQUIRES(
  (concepts::is_linear_operator_on<MatrixType, VectorType> &&
   concepts::is_linear_operator_on<PreconditionerType, VectorType>))
void SolverBicgstab<VectorType>::solve(const MatrixType         &A,
                                       VectorType               &x,
                                       const VectorType         &b,
                                       const PreconditionerType &preconditioner)
{
  LogStream::Prefix prefix("Bicgstab");
 
  IterationResult state(false, SolverControl::failure, 0, 0);
  do
    {
      state = iterate(A, x, b, preconditioner, state.last_step);
    }
  while (state.state == SolverControl::iterate);
 
  // In case of failure: throw exception
  AssertThrow(state.state == SolverControl::success,
              SolverControl::NoConvergence(state.last_step,
                                           state.last_residual));
  // Otherwise exit as normal
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif