solver_richardson.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_richardson_h
#define dealii__solver_richardson_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 <deal.II/base/subscriptor.h>

DEAL_II_NAMESPACE_OPEN


template <class VectorType = Vector<double> >
class SolverRichardson : public Solver<VectorType>
{
public:
  struct AdditionalData
  {
    explicit
    AdditionalData (const double omega                       = 1,
                    const bool   use_preconditioned_residual = false);

    double omega;

    bool use_preconditioned_residual;

  };

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

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

  virtual ~SolverRichardson ();

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

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

  void set_omega (const double om=1.);

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

protected:
  virtual typename VectorType::value_type criterion();

  VectorType *Vr;
  VectorType *Vd;

  AdditionalData additional_data;

  typename VectorType::value_type res;
};

/*----------------- Implementation of the Richardson Method ------------------*/

#ifndef DOXYGEN

template <class VectorType>
inline
SolverRichardson<VectorType>::AdditionalData::
AdditionalData (const double omega,
                const bool   use_preconditioned_residual)
  :
  omega(omega),
  use_preconditioned_residual(use_preconditioned_residual)
{}


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



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



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


template <class VectorType>
template <typename MatrixType, typename PreconditionerType>
void
SolverRichardson<VectorType>::solve (const MatrixType         &A,
                                     VectorType               &x,
                                     const VectorType         &b,
                                     const PreconditionerType &precondition)
{
  SolverControl::State conv=SolverControl::iterate;

  double last_criterion = -std::numeric_limits<double>::max();

  unsigned int iter = 0;

  // Memory allocation
  Vr  = this->memory.alloc();
  VectorType &r  = *Vr;
  r.reinit(x);
  Vd  = this->memory.alloc();
  VectorType &d  = *Vd;
  d.reinit(x);

  deallog.push("Richardson");

  try
    {
      // Main loop
      while (conv==SolverControl::iterate)
        {
          // Do not use residual,
          // but do it in 2 steps
          A.vmult(r,x);
          r.sadd(-1.,1.,b);
          precondition.vmult(d,r);

          // The required norm of the
          // (preconditioned)
          // residual is computed in
          // criterion() and stored
          // in res.
          last_criterion = criterion();
          conv = this->iteration_status (iter, last_criterion, x);
          if (conv != SolverControl::iterate)
            break;

          x.add(additional_data.omega,d);
          print_vectors(iter,x,r,d);

          ++iter;
        }
    }
  catch (...)
    {
      this->memory.free(Vr);
      this->memory.free(Vd);
      deallog.pop();
      throw;
    }
  // Deallocate Memory
  this->memory.free(Vr);
  this->memory.free(Vd);
  deallog.pop();

  // in case of failure: throw exception
  if (conv != SolverControl::success)
    AssertThrow(false, SolverControl::NoConvergence (iter,
                                                     last_criterion));
  // otherwise exit as normal
}


template <class VectorType>
template <typename MatrixType, typename PreconditionerType>
void
SolverRichardson<VectorType>::Tsolve (const MatrixType         &A,
                                      VectorType               &x,
                                      const VectorType         &b,
                                      const PreconditionerType &precondition)
{
  SolverControl::State conv=SolverControl::iterate;
  double last_criterion = -std::numeric_limits<double>::max();

  unsigned int iter = 0;

  // Memory allocation
  Vr  = this->memory.alloc();
  VectorType &r  = *Vr;
  r.reinit(x);
  Vd  =this-> memory.alloc();
  VectorType &d  = *Vd;
  d.reinit(x);

  deallog.push("RichardsonT");

  try
    {
      // Main loop
      while (conv==SolverControl::iterate)
        {
          // Do not use Tresidual,
          // but do it in 2 steps
          A.Tvmult(r,x);
          r.sadd(-1.,1.,b);
          precondition.Tvmult(d,r);

          last_criterion = criterion();
          conv = this->iteration_status (iter, last_criterion, x);
          if (conv != SolverControl::iterate)
            break;

          x.add(additional_data.omega,d);
          print_vectors(iter,x,r,d);

          ++iter;
        }
    }
  catch (...)
    {
      this->memory.free(Vr);
      this->memory.free(Vd);
      deallog.pop();
      throw;
    }

  // Deallocate Memory
  this->memory.free(Vr);
  this->memory.free(Vd);
  deallog.pop();
  // in case of failure: throw exception
  if (conv != SolverControl::success)
    AssertThrow(false, SolverControl::NoConvergence (iter, last_criterion));

  // otherwise exit as normal
}


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



template <class VectorType>
inline typename VectorType::value_type
SolverRichardson<VectorType>::criterion()
{
  if (!additional_data.use_preconditioned_residual)
    res = Vr->l2_norm();
  else
    res = Vd->l2_norm();
  return res;
}


template <class VectorType>
inline void
SolverRichardson<VectorType>::set_omega (const double om)
{
  additional_data.omega=om;
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif