solver_bicgstab.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 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_bicgstab_h
#define dealii__solver_bicgstab_h


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

DEAL_II_NAMESPACE_OPEN


template <typename VectorType = Vector<double> >
class SolverBicgstab : public Solver<VectorType>
{
public:
  struct AdditionalData
  {
    explicit
    AdditionalData(const bool   exact_residual = true,
                   const double breakdown      = 1.e-10)
      : exact_residual(exact_residual),
        breakdown(breakdown)
    {}
    bool exact_residual;
    double breakdown;
  };

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

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

  virtual ~SolverBicgstab ();

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

protected:
  template <typename MatrixType>
  double criterion (const MatrixType &A, const VectorType &x, const VectorType &b);

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

  VectorType *Vx;
  VectorType *Vr;
  VectorType *Vrbar;
  VectorType *Vp;
  VectorType *Vy;
  VectorType *Vz;
  VectorType *Vt;
  VectorType *Vv;
  const VectorType *Vb;

  double alpha;
  double beta;
  double omega;
  double rho;
  double rhobar;

  unsigned int step;

  double res;

  AdditionalData additional_data;

private:
  template <typename MatrixType>
  SolverControl::State start(const MatrixType &A);

  struct IterationResult
  {
    bool                 breakdown;
    SolverControl::State state;
    unsigned int         last_step;
    double               last_residual;

    IterationResult (const bool breakdown,
                     const SolverControl::State state,
                     const unsigned int         last_step,
                     const double               last_residual);
  };

  template<typename MatrixType, typename PreconditionerType>
  IterationResult
  iterate(const MatrixType         &A,
          const PreconditionerType &precondition);
};

/*-------------------------Inline functions -------------------------------*/

#ifndef DOXYGEN


template<typename VectorType>
SolverBicgstab<VectorType>::IterationResult::IterationResult
(const bool                 breakdown,
 const SolverControl::State state,
 const unsigned int         last_step,
 const double               last_residual)
  :
  breakdown (breakdown),
  state (state),
  last_step (last_step),
  last_residual (last_residual)
{}


template<typename VectorType>
SolverBicgstab<VectorType>::SolverBicgstab (SolverControl            &cn,
                                            VectorMemory<VectorType> &mem,
                                            const AdditionalData     &data)
  :
  Solver<VectorType>(cn,mem),
  additional_data(data)
{}



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



template<typename VectorType>
SolverBicgstab<VectorType>::~SolverBicgstab ()
{}



template <typename VectorType>
template <typename MatrixType>
double
SolverBicgstab<VectorType>::criterion (const MatrixType &A,
                                       const VectorType &x,
                                       const VectorType &b)
{
  A.vmult(*Vt, x);
  Vt->add(-1.,b);
  res = Vt->l2_norm();

  return res;
}



template <typename VectorType >
template <typename MatrixType>
SolverControl::State
SolverBicgstab<VectorType>::start (const MatrixType &A)
{
  A.vmult(*Vr, *Vx);
  Vr->sadd(-1.,1.,*Vb);
  res = Vr->l2_norm();

  return this->iteration_status(step, res, *Vx);
}



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



template<typename VectorType>
template<typename MatrixType, typename PreconditionerType>
typename SolverBicgstab<VectorType>::IterationResult
SolverBicgstab<VectorType>::iterate(const MatrixType         &A,
                                    const PreconditionerType &precondition)
{
//TODO:[GK] Implement "use the length of the computed orthogonal residual" in the BiCGStab method.
  SolverControl::State state = SolverControl::iterate;
  alpha = omega = rho = 1.;

  VectorType &r = *Vr;
  VectorType &rbar = *Vrbar;
  VectorType &p = *Vp;
  VectorType &y = *Vy;
  VectorType &z = *Vz;
  VectorType &t = *Vt;
  VectorType &v = *Vv;

  rbar = r;
  bool startup = true;

  do
    {
      ++step;

      rhobar = r*rbar;
      beta   = rhobar * alpha / (rho * omega);
      rho    = rhobar;
      if (startup == true)
        {
          p = r;
          startup = false;
        }
      else
        {
          p.sadd(beta, 1., r);
          p.add(-beta*omega, v);
        }

      precondition.vmult(y,p);
      A.vmult(v,y);
      rhobar = rbar * v;

      alpha = rho/rhobar;

//TODO:[?] Find better breakdown criterion

      if (std::fabs(alpha) > 1.e10)
        return IterationResult(true, state, step, res);

      res = std::sqrt(r.add_and_dot(-alpha, v, r));

      // check for early success, see the lac/bicgstab_early testcase as to
      // why this is necessary
      //
      // note: the vector *Vx we pass to the iteration_status signal here is only
      // the current approximation, not the one we will return with,
      // which will be x=*Vx + alpha*y
      if (this->iteration_status(step, res, *Vx) == SolverControl::success)
        {
          Vx->add(alpha, y);
          print_vectors(step, *Vx, r, y);
          return IterationResult(false, SolverControl::success, step, res);
        }

      precondition.vmult(z,r);
      A.vmult(t,z);
      rhobar = t*r;
      omega = rhobar/(t*t);
      Vx->add(alpha, y, omega, z);

      if (additional_data.exact_residual)
        {
          r.add(-omega, t);
          res = criterion(A, *Vx, *Vb);
        }
      else
        res = std::sqrt(r.add_and_dot(-omega, t, r));

      state = this->iteration_status(step, res, *Vx);
      print_vectors(step, *Vx, r, y);
    }
  while (state == SolverControl::iterate);
  return IterationResult(false, state, step, res);
}


template<typename VectorType>
template<typename MatrixType, typename PreconditionerType>
void
SolverBicgstab<VectorType>::solve(const MatrixType         &A,
                                  VectorType               &x,
                                  const VectorType         &b,
                                  const PreconditionerType &precondition)
{
  deallog.push("Bicgstab");
  Vr    = this->memory.alloc();
  Vr->reinit(x, true);
  Vrbar = this->memory.alloc();
  Vrbar->reinit(x, true);
  Vp    = this->memory.alloc();
  Vp->reinit(x, true);
  Vy    = this->memory.alloc();
  Vy->reinit(x, true);
  Vz    = this->memory.alloc();
  Vz->reinit(x, true);
  Vt    = this->memory.alloc();
  Vt->reinit(x, true);
  Vv    = this->memory.alloc();
  Vv->reinit(x, true);

  Vx = &x;
  Vb = &b;

  step = 0;

  IterationResult state(false,SolverControl::failure,0,0);

  // iterate while the inner iteration returns a breakdown
  do
    {
      if (step != 0)
        deallog << "Restart step " << step << std::endl;
      if (start(A) == SolverControl::success)
        {
          state.state = SolverControl::success;
          break;
        }
      state = iterate(A, precondition);
    }
  while (state.breakdown == true);

  this->memory.free(Vr);
  this->memory.free(Vrbar);
  this->memory.free(Vp);
  this->memory.free(Vy);
  this->memory.free(Vz);
  this->memory.free(Vt);
  this->memory.free(Vv);

  deallog.pop();

  // in case of failure: throw exception
  AssertThrow(state.state == SolverControl::success,
              SolverControl::NoConvergence (state.last_step,
                                            state.last_residual));
  // otherwise exit as normal
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif