solver_qmrs.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 1999 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_solver_qmrs_h
#define dealii_solver_qmrs_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/template_constraints.h>
 
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
 
#include <cmath>
 
DEAL_II_NAMESPACE_OPEN


template <typename VectorType = Vector<double>>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
class SolverQMRS : public SolverBase<VectorType>
{
public:
  struct AdditionalData
  {
    explicit AdditionalData(const bool   left_preconditioning = false,
                            const double solver_tolerance     = 1.e-9,
                            const bool   breakdown_testing    = true,
                            const double breakdown_threshold  = 1.e-16)
      : left_preconditioning(left_preconditioning)
      , solver_tolerance(solver_tolerance)
      , breakdown_testing(breakdown_testing)
      , breakdown_threshold(breakdown_threshold)
    {}

    bool left_preconditioning;

    double solver_tolerance;

    bool breakdown_testing;

    double breakdown_threshold;
  };

  SolverQMRS(SolverControl            &cn,
             VectorMemory<VectorType> &mem,
             const AdditionalData     &data = AdditionalData());

  SolverQMRS(SolverControl &cn, const AdditionalData &data = AdditionalData());

  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);

  virtual void
  print_vectors(const unsigned int step,
                const VectorType  &x,
                const VectorType  &r,
                const VectorType  &d) const;
 
protected:
  AdditionalData additional_data;
 
private:
  struct IterationResult
  {
    SolverControl::State state;
    double               last_residual;
 
    IterationResult(const SolverControl::State state,
                    const double               last_residual);
  };

  template <typename MatrixType, typename PreconditionerType>
  IterationResult
  iterate(const MatrixType         &A,
          VectorType               &x,
          const VectorType         &b,
          const PreconditionerType &preconditioner,
          VectorType               &r,
          VectorType               &u,
          VectorType               &q,
          VectorType               &t,
          VectorType               &d);

  unsigned int step;
};

/*------------------------- Implementation ----------------------------*/
 
#ifndef DOXYGEN
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverQMRS<VectorType>::IterationResult::IterationResult(
  const SolverControl::State state,
  const double               last_residual)
  : state(state)
  , last_residual(last_residual)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverQMRS<VectorType>::SolverQMRS(SolverControl            &cn,
                                   VectorMemory<VectorType> &mem,
                                   const AdditionalData     &data)
  : SolverBase<VectorType>(cn, mem)
  , additional_data(data)
  , step(0)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverQMRS<VectorType>::SolverQMRS(SolverControl        &cn,
                                   const AdditionalData &data)
  : SolverBase<VectorType>(cn)
  , additional_data(data)
  , step(0)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
void SolverQMRS<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>
DEAL_II_CXX20_REQUIRES(
  (concepts::is_linear_operator_on<MatrixType, VectorType> &&
   concepts::is_linear_operator_on<PreconditionerType, VectorType>))
void SolverQMRS<VectorType>::solve(const MatrixType         &A,
                                   VectorType               &x,
                                   const VectorType         &b,
                                   const PreconditionerType &preconditioner)
{
  LogStream::Prefix prefix("SQMR");
 
 
  // temporary vectors, allocated through the @p VectorMemory object at the
  // start of the actual solution process and deallocated at the end.
  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
  typename VectorMemory<VectorType>::Pointer Vu(this->memory);
  typename VectorMemory<VectorType>::Pointer Vq(this->memory);
  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
 
 
  // resize the vectors, but do not set
  // the values since they'd be overwritten
  // soon anyway.
  Vr->reinit(x, true);
  Vu->reinit(x, true);
  Vq->reinit(x, true);
  Vt->reinit(x, true);
  Vd->reinit(x, true);
 
  step = 0;
 
  IterationResult state(SolverControl::failure, 0);
 
  do
    {
      if (step > 0)
        deallog << "Restart step " << step << std::endl;
      state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
    }
  while (state.state == SolverControl::iterate);
 
 
  // in case of failure: throw exception
  AssertThrow(state.state == SolverControl::success,
              SolverControl::NoConvergence(step, state.last_residual));
  // otherwise exit as normal
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType, typename PreconditionerType>
typename SolverQMRS<VectorType>::IterationResult
  SolverQMRS<VectorType>::iterate(const MatrixType         &A,
                                  VectorType               &x,
                                  const VectorType         &b,
                                  const PreconditionerType &preconditioner,
                                  VectorType               &r,
                                  VectorType               &u,
                                  VectorType               &q,
                                  VectorType               &t,
                                  VectorType               &d)
{
  SolverControl::State state = SolverControl::iterate;
 
  int it = 0;
 
  double tau, rho, theta = 0;
  double res;
 
  // Compute the start residual
  A.vmult(r, x);
  r.sadd(-1., 1., b);
 
  // Doing the initial preconditioning
  if (additional_data.left_preconditioning)
    {
      // Left preconditioning
      preconditioner.vmult(t, r);
      q = t;
    }
  else
    {
      // Right preconditioning
      t = r;
      preconditioner.vmult(q, t);
    }
 
  tau = t.norm_sqr();
  res = std::sqrt(tau);
 
  if (this->iteration_status(step, res, x) == SolverControl::success)
    return IterationResult(SolverControl::success, res);
 
  rho = q * r;
 
  while (state == SolverControl::iterate)
    {
      ++step;
      ++it;
      //--------------------------------------------------------------
      // Step 1: apply the system matrix and compute one inner product
      //--------------------------------------------------------------
      A.vmult(t, q);
      const double sigma = q * t;
 
      // Check the breakdown criterion
      if (additional_data.breakdown_testing == true &&
          std::fabs(sigma) < additional_data.breakdown_threshold)
        return IterationResult(SolverControl::iterate, res);
      // Update the residual
      const double alpha = rho / sigma;
      r.add(-alpha, t);
 
      //--------------------------------------------------------------
      // Step 2: update the solution vector
      //--------------------------------------------------------------
      const double theta_old = theta;
 
      // Apply the preconditioner
      if (additional_data.left_preconditioning)
        {
          // Left Preconditioning
          preconditioner.vmult(t, r);
        }
      else
        {
          // Right Preconditioning
          t = r;
        }
 
      // Double updates
      theta            = t * t / tau;
      const double psi = 1. / (1. + theta);
      tau *= theta * psi;
 
      // Actual update of the solution vector
      d.sadd(psi * theta_old, psi * alpha, q);
      x += d;
 
      print_vectors(step, x, r, d);
 
      // Check for convergence
      // Compute a simple and cheap upper bound of the norm of the residual
      // vector b-Ax
      res = std::sqrt((it + 1) * tau);
      // If res lies close enough, within the desired tolerance, calculate the
      // exact residual
      if (res < additional_data.solver_tolerance)
        {
          A.vmult(u, x);
          u.sadd(-1., 1., b);
          res = u.l2_norm();
        }
      state = this->iteration_status(step, res, x);
      if ((state == SolverControl::success) ||
          (state == SolverControl::failure))
        return IterationResult(state, res);
 
      //--------------------------------------------------------------
      // Step 3: check breakdown criterion and update the vectors
      //--------------------------------------------------------------
      if (additional_data.breakdown_testing == true &&
          std::fabs(sigma) < additional_data.breakdown_threshold)
        return IterationResult(SolverControl::iterate, res);
 
      const double rho_old = rho;
 
      // Applying the preconditioner
      if (additional_data.left_preconditioning)
        {
          // Left preconditioning
          u = t;
        }
      else
        {
          // Right preconditioning
          preconditioner.vmult(u, t);
        }
 
      // Double and vector updates
      rho               = u * r;
      const double beta = rho / rho_old;
      q.sadd(beta, 1., u);
    }
  return IterationResult(SolverControl::success, res);
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif