solver_qmrs.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2017 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(bool exact_residual = false,
                   double breakdown=1.e-16) :
      exact_residual(exact_residual),
      breakdown(breakdown)
    {}

    bool exact_residual;

    double breakdown;
  };

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

  virtual void print_vectors (const unsigned int step,
                              const VectorType   &x,
                              const VectorType   &r,
                              const VectorType   &d) const;
protected:
  virtual double criterion();

  VectorType *Vv;
  VectorType *Vp;
  VectorType *Vq;
  VectorType *Vt;
  VectorType *Vd;
  VectorType *Vx;
  const VectorType *Vb;

  double res2;

  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,
           const PreconditionerType &precondition);

  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)
{}



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



template<class VectorType>
double
SolverQMRS<VectorType>::criterion()
{
  return std::sqrt(res2);
}



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 &precondition)
{
  deallog.push("QMRS");

  // Memory allocation
  Vv  = this->memory.alloc();
  Vp  = this->memory.alloc();
  Vq  = this->memory.alloc();
  Vt  = this->memory.alloc();
  Vd  = this->memory.alloc();

  Vx = &x;
  Vb = &b;
  // resize the vectors, but do not set
  // the values since they'd be overwritten
  // soon anyway.
  Vv->reinit(x, true);
  Vp->reinit(x, true);
  Vq->reinit(x, true);
  Vt->reinit(x, true);

  step = 0;

  IterationResult state (SolverControl::failure,0);

  do
    {
      if (step > 0)
        deallog << "Restart step " << step << std::endl;
      state = iterate(A, precondition);
    }
  while (state.state == SolverControl::iterate);

  // Deallocate Memory
  this->memory.free(Vv);
  this->memory.free(Vp);
  this->memory.free(Vq);
  this->memory.free(Vt);
  this->memory.free(Vd);

  // Output
  deallog.pop();

  // 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,
                                const PreconditionerType &precondition)
{
  /* Remark: the matrix A in the article is the preconditioned matrix.
   * Therefore, we have to precondition x before we compute the first residual.
   * In step 1 we replace p by q to avoid one preconditioning step.
   * There are still two steps left, making this algorithm expensive.
   */

  SolverControl::State state = SolverControl::iterate;

  // define some aliases for simpler access
  VectorType &v  = *Vv;
  VectorType &p  = *Vp;
  VectorType &q  = *Vq;
  VectorType &t  = *Vt;
  VectorType &d  = *Vd;
  VectorType &x  = *Vx;
  const VectorType &b = *Vb;

  int  it=0;

  double tau, rho, theta=0;
  double res;

  d.reinit(x);

  // Apply right preconditioning to x
  precondition.vmult(q,x);
  // Preconditioned residual
  A.vmult(v,q);
  v.sadd(-1.,1.,b);
  res = v.l2_norm();

  if (this->iteration_status(step, res, x) == SolverControl::success)
    return IterationResult(SolverControl::success, res);

  p = v;

  precondition.vmult(q,p);

  tau = v.norm_sqr();
  rho = q*v;

  while (state == SolverControl::iterate)
    {
      step++;
      it++;
      // Step 1
      A.vmult(t,q);
      // Step 2
      const double sigma = q*t;

//TODO:[?] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence
      if (std::fabs(sigma/rho) < additional_data.breakdown)
        return IterationResult(SolverControl::iterate, std::fabs(sigma/rho));
      // Step 3
      const double alpha = rho/sigma;

      v.add(-alpha,t);
      // Step 4
      const double theta_old = theta;
      theta = v*v/tau;
      const double psi = 1./(1.+theta);
      tau *= theta*psi;

      d.sadd(psi*theta_old, psi*alpha, p);
      x.add(d);

      print_vectors(step,x,v,d);
      // Step 5
      if (additional_data.exact_residual)
        {
          A.vmult(q,x);
          q.sadd(-1.,1.,b);
          res = q.l2_norm();
        }
      else
        res = std::sqrt((it+1)*tau);
      state = this->iteration_status(step,res,x);
      if ((state == SolverControl::success)
          || (state == SolverControl::failure))
        return IterationResult(state, res);
      // Step 6
      if (std::fabs(rho) < additional_data.breakdown)
        return IterationResult(SolverControl::iterate, std::fabs(rho));
      // Step 7
      const double rho_old = rho;
      precondition.vmult(q,v);
      rho = q*v;

      const double beta = rho/rho_old;
      p.sadd(beta,v);
      precondition.vmult(q,p);
    }
  return IterationResult(SolverControl::success, std::fabs(rho));
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif