solver_qmrs.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_solver_qmrs_h
#define dealii_solver_qmrs_h

#include <deal.II/base/config.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/base/logstream.h>
#include <cmath>
#include <deal.II/base/subscriptor.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


template <typename VectorType = Vector<double> >
class SolverQMRS : public Solver<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>
  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 <class VectorType>
SolverQMRS<VectorType>::IterationResult::IterationResult (const SolverControl::State state,
                                                          const double last_residual)
  :
  state (state),
  last_residual (last_residual)
{}



template <class VectorType>
SolverQMRS<VectorType>::SolverQMRS (SolverControl            &cn,
                                    VectorMemory<VectorType> &mem,
                                    const AdditionalData     &data)
  :
  Solver<VectorType> (cn, mem),
  additional_data (data),
  step (0)
{
}

template <class VectorType>
SolverQMRS<VectorType>::SolverQMRS (SolverControl        &cn,
                                    const AdditionalData &data)
  :
  Solver<VectorType> (cn),
  additional_data (data),
  step (0)
{
}

template <class VectorType>
void
SolverQMRS<VectorType>::print_vectors (const unsigned int,
                                       const VectorType &,
                                       const VectorType &,
                                       const VectorType &) const
{
}

template <class VectorType>
template <typename MatrixType, typename PreconditionerType>
void
SolverQMRS<VectorType>::solve (const MatrixType &A,
                               VectorType &x,
                               const VectorType &b,
                               const PreconditionerType &preconditioner)
{
  LogStream::Prefix prefix("SQMR");


// temporary vectors, allocated trough 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 <class 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, u);
    }
  return IterationResult (SolverControl::success, res);
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif