solver_minres.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2015 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_minres_h
#define dealii__solver_minres_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>

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

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

  DeclException0 (ExcPreconditionerNotDefinite);

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

  VectorType *Vu0, *Vu1, *Vu2;
  VectorType *Vm0, *Vm1, *Vm2;
  VectorType *Vv;

  double res2;
};

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

#ifndef DOXYGEN

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



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


template<class VectorType>
SolverMinRes<VectorType>::~SolverMinRes ()
{}



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 &precondition)
{
  SolverControl::State conv=SolverControl::iterate;

  deallog.push("minres");

  // Memory allocation
  Vu0  = this->memory.alloc();
  Vu1  = this->memory.alloc();
  Vu2  = this->memory.alloc();
  Vv   = this->memory.alloc();
  Vm0  = this->memory.alloc();
  Vm1  = this->memory.alloc();
  Vm2  = this->memory.alloc();
  // define some aliases for simpler access
  typedef VectorType *vecptr;
  vecptr u[3] = {Vu0, Vu1, Vu2};
  vecptr m[3] = {Vm0, Vm1, Vm2};
  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 gamma = 0;
  double s = 0;
  double d_ = 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]
  precondition.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);

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

      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.
      precondition.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]);
        }

      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];
    }

  // Deallocation of Memory
  this->memory.free(Vu0);
  this->memory.free(Vu1);
  this->memory.free(Vu2);
  this->memory.free(Vv);
  this->memory.free(Vm0);
  this->memory.free(Vm1);
  this->memory.free(Vm2);
  // Output
  deallog.pop ();

  // 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