solver_minres.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_solver_minres_h
#define dealii_solver_minres_h


#include <deal.II/base/config.h>

#include <deal.II/base/logstream.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


template <class VectorType = Vector<double>>
class SolverMinRes : public Solver<VectorType>
{
public:
  struct AdditionalData
  {};

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

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

  virtual ~SolverMinRes() override = default;

  template <typename MatrixType, typename PreconditionerType>
  void
  solve(const MatrixType &        A,
        VectorType &              x,
        const VectorType &        b,
        const PreconditionerType &preconditioner);

  DeclException0(ExcPreconditionerNotDefinite);

protected:
  virtual double
  criterion();

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

  double res2;
};

/*------------------------- Implementation ----------------------------*/

#ifndef DOXYGEN

template <class VectorType>
SolverMinRes<VectorType>::SolverMinRes(SolverControl &           cn,
                                       VectorMemory<VectorType> &mem,
                                       const AdditionalData &)
  : Solver<VectorType>(cn, mem)
  , res2(numbers::signaling_nan<double>())
{}



template <class VectorType>
SolverMinRes<VectorType>::SolverMinRes(SolverControl &cn,
                                       const AdditionalData &)
  : Solver<VectorType>(cn)
  , res2(numbers::signaling_nan<double>())
{}



template <class VectorType>
double
SolverMinRes<VectorType>::criterion()
{
  return res2;
}


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



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

  // Memory allocation
  typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
  typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
  typename VectorMemory<VectorType>::Pointer Vu2(this->memory);

  typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
  typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
  typename VectorMemory<VectorType>::Pointer Vm2(this->memory);

  typename VectorMemory<VectorType>::Pointer Vv(this->memory);

  // define some aliases for simpler access
  using vecptr     = VectorType *;
  vecptr      u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
  vecptr      m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
  VectorType &v    = *Vv;

  // resize the vectors, but do not set the values since they'd be overwritten
  // soon anyway.
  u[0]->reinit(b, true);
  u[1]->reinit(b, true);
  u[2]->reinit(b, true);
  m[0]->reinit(b, true);
  m[1]->reinit(b, true);
  m[2]->reinit(b, true);
  v.reinit(b, true);

  // some values needed
  double delta[3] = {0, 0, 0};
  double f[2]     = {0, 0};
  double e[2]     = {0, 0};

  double r_l2 = 0;
  double r0   = 0;
  double tau  = 0;
  double c    = 0;
  double s    = 0;
  double d_   = 0;

  // The iteration step.
  unsigned int j = 1;


  // Start of the solution process
  A.vmult(*m[0], x);
  *u[1] = b;
  *u[1] -= *m[0];
  // Precondition is applied.
  // The preconditioner has to be
  // positive definite and symmetric

  // M v = u[1]
  preconditioner.vmult(v, *u[1]);

  delta[1] = v * (*u[1]);
  // Preconditioner positive
  Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());

  r0   = std::sqrt(delta[1]);
  r_l2 = r0;


  u[0]->reinit(b);
  delta[0] = 1.;
  m[0]->reinit(b);
  m[1]->reinit(b);
  m[2]->reinit(b);

  SolverControl::State conv = this->iteration_status(0, r_l2, x);
  while (conv == SolverControl::iterate)
    {
      if (delta[1] != 0)
        v *= 1. / std::sqrt(delta[1]);
      else
        v.reinit(b);

      A.vmult(*u[2], v);
      u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);

      const double gamma = *u[2] * v;
      u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
      *m[0] = v;

      // precondition: solve M v = u[2]
      // Preconditioner has to be positive
      // definite and symmetric.
      preconditioner.vmult(v, *u[2]);

      delta[2] = v * (*u[2]);

      Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());

      if (j == 1)
        {
          d_   = gamma;
          e[1] = std::sqrt(delta[2]);
        }
      if (j > 1)
        {
          d_   = s * e[0] - c * gamma;
          e[0] = c * e[0] + s * gamma;
          f[1] = s * std::sqrt(delta[2]);
          e[1] = -c * std::sqrt(delta[2]);
        }

      const double d = std::sqrt(d_ * d_ + delta[2]);

      if (j > 1)
        tau *= s / c;
      c = d_ / d;
      tau *= c;

      s = std::sqrt(delta[2]) / d;

      if (j == 1)
        tau = r0 * c;

      m[0]->add(-e[0], *m[1]);
      if (j > 1)
        m[0]->add(-f[0], *m[2]);
      *m[0] *= 1. / d;
      x.add(tau, *m[0]);
      r_l2 *= std::fabs(s);

      conv = this->iteration_status(j, r_l2, x);

      // next iteration step
      ++j;
      // All vectors have to be shifted
      // one iteration step.
      // This should be changed one time.
      //
      // the previous code was like this:
      //   m[2] = m[1];
      //   m[1] = m[0];
      // but it can be made more efficient,
      // since the value of m[0] is no more
      // needed in the next iteration
      swap(*m[2], *m[1]);
      swap(*m[1], *m[0]);

      // likewise, but reverse direction:
      //   u[0] = u[1];
      //   u[1] = u[2];
      swap(*u[0], *u[1]);
      swap(*u[1], *u[2]);

      // these are scalars, so need
      // to bother
      f[0]     = f[1];
      e[0]     = e[1];
      delta[0] = delta[1];
      delta[1] = delta[2];
    }

  // in case of failure: throw exception
  AssertThrow(conv == SolverControl::success,
              SolverControl::NoConvergence(j, r_l2));

  // otherwise exit as normal
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif