solver.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_ISTL_SOLVER_HH
#define DUNE_ISTL_SOLVER_HH
 
#include <iomanip>
#include <ostream>
#include <string>
#include <functional>
 
#include <dune/common/exceptions.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/simd/io.hh>
#include <dune/common/simd/simd.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/timer.hh>
 
#include "solvertype.hh"
#include "preconditioner.hh"
#include "operators.hh"
#include "scalarproducts.hh"
 
namespace Dune
{
  struct InverseOperatorResult
  {
    InverseOperatorResult ()
    {
      clear();
    }
 
    void clear ()
    {
      iterations = 0;
      reduction = 0;
      converged = false;
      conv_rate = 1;
      elapsed = 0;
      condition_estimate = -1;
    }
 
    int iterations;
 
    double reduction;
 
    bool converged;
 
    double conv_rate;
 
    double condition_estimate = -1;
 
    double elapsed;
  };
 
 
  //=====================================================================
  template<class X, class Y>
  class InverseOperator {
  public:
    typedef X domain_type;
 
    typedef Y range_type;
 
    typedef typename X::field_type field_type;
 
    typedef typename FieldTraits<field_type>::real_type real_type;
 
    typedef Simd::Scalar<real_type> scalar_real_type;
 
    virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
 
    virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
 
    virtual SolverCategory::Category category() const
#ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
    {
      DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
    }
#else
    = 0;
#endif
 
    virtual ~InverseOperator () {}
 
  protected:
    // spacing values
    enum { iterationSpacing = 5 , normSpacing = 16 };
 
    void printHeader(std::ostream& s) const
    {
      s << std::setw(iterationSpacing)  << " Iter";
      s << std::setw(normSpacing) << "Defect";
      s << std::setw(normSpacing) << "Rate" << std::endl;
    }
 
    template <typename CountType, typename DataType>
    void printOutput(std::ostream& s,
                     const CountType& iter,
                     const DataType& norm,
                     const DataType& norm_old) const
    {
      const DataType rate = norm/norm_old;
      s << std::setw(iterationSpacing)  << iter << " ";
      s << std::setw(normSpacing) << Simd::io(norm) << " ";
      s << std::setw(normSpacing) << Simd::io(rate) << std::endl;
    }
 
    template <typename CountType, typename DataType>
    void printOutput(std::ostream& s,
                     const CountType& iter,
                     const DataType& norm) const
    {
      s << std::setw(iterationSpacing)  << iter << " ";
      s << std::setw(normSpacing) << Simd::io(norm) << std::endl;
    }
  };
 
  template<class X, class Y>
  class IterativeSolver : public InverseOperator<X,Y>{
  public:
    using typename InverseOperator<X,Y>::domain_type;
    using typename InverseOperator<X,Y>::range_type;
    using typename InverseOperator<X,Y>::field_type;
    using typename InverseOperator<X,Y>::real_type;
    using typename InverseOperator<X,Y>::scalar_real_type;
 
    IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
      _op(stackobject_to_shared_ptr(op)),
      _prec(stackobject_to_shared_ptr(prec)),
      _sp(new SeqScalarProduct<X>),
      _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
    {
      if(SolverCategory::category(op) != SolverCategory::sequential)
        DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
      if(SolverCategory::category(prec) != SolverCategory::sequential)
        DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
    }
 
    IterativeSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
      scalar_real_type reduction, int maxit, int verbose) :
      _op(stackobject_to_shared_ptr(op)),
      _prec(stackobject_to_shared_ptr(prec)),
      _sp(stackobject_to_shared_ptr(sp)),
      _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
    {
      if(SolverCategory::category(op) != SolverCategory::category(prec))
        DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
      if(SolverCategory::category(op) != SolverCategory::category(sp))
        DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
    }
 
    IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
      IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
        configuration.get<real_type>("reduction"),
        configuration.get<int>("maxit"),
        configuration.get<int>("verbose"))
    {}
 
    IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
      IterativeSolver(op,sp,prec,
        configuration.get<scalar_real_type>("reduction"),
        configuration.get<int>("maxit"),
        configuration.get<int>("verbose"))
    {}
 
    IterativeSolver (std::shared_ptr<LinearOperator<X,Y>> op,
                     std::shared_ptr<ScalarProduct<X>> sp,
                     std::shared_ptr<Preconditioner<X,Y>> prec,
                     scalar_real_type reduction, int maxit, int verbose) :
      _op(op),
      _prec(prec),
      _sp(sp),
      _reduction(reduction), _maxit(maxit), _verbose(verbose),
      _category(SolverCategory::category(*op))
    {
      if(SolverCategory::category(*op) != SolverCategory::category(*prec))
        DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
      if(SolverCategory::category(*op) != SolverCategory::category(*sp))
        DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
    }
 
    // #warning actually we want to have this as the default and just implement the second one
    // //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
    // virtual void apply (X& x, Y& b, InverseOperatorResult& res)
    // {
    //   apply(x,b,_reduction,res);
    // }
 
#ifndef DOXYGEN
    // make sure the three-argument apply from the base class does not get shadowed
    // by the redefined four-argument version below
    using InverseOperator<X,Y>::apply;
#endif
 
    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
    {
      scalar_real_type saved_reduction = _reduction;
      _reduction = reduction;
      this->apply(x,b,res);
      _reduction = saved_reduction;
    }
 
    virtual SolverCategory::Category category() const
    {
      return _category;
    }
 
    std::string name() const{
      std::string name = className(*this);
      return name.substr(0, name.find("<"));
    }
 
    template<class CountType = unsigned int>
    class Iteration {
    public:
      Iteration(const IterativeSolver& parent, InverseOperatorResult& res)
        : _i(0)
        , _res(res)
        , _parent(parent)
        , _valid(true)
      {
        res.clear();
        if(_parent._verbose>0){
          std::cout << "=== " << parent.name() << std::endl;
          if(_parent._verbose > 1)
            _parent.printHeader(std::cout);
        }
      }
 
      Iteration(const Iteration&) = delete;
      Iteration(Iteration&& other)
        : _def0(other._def0)
        , _def(other._def)
        , _i(other._i)
        , _watch(other._watch)
        , _res(other._res)
        , _parent(other._parent)
        , _valid(other._valid)
      {
        other._valid = false;
      }
 
      ~Iteration(){
        if(_valid)
          finalize();
      }
 
      bool step(CountType i, real_type def){
        if (!Simd::allTrue(isFinite(def))) // check for inf or NaN
        {
          if (_parent._verbose>0)
            std::cout << "=== " << _parent.name() << ": abort due to infinite or NaN defect"
                      << std::endl;
          DUNE_THROW(SolverAbort,
                     _parent.name() << ": defect=" << Simd::io(def)
                     << " is infinite or NaN");
        }
        if(i == 0)
          _def0 = def;
        if(_parent._verbose > 1){
          if(i!=0)
            _parent.printOutput(std::cout,i,def,_def);
          else
            _parent.printOutput(std::cout,i,def);
        }
        _def = def;
        _i = i;
        _res.converged = (Simd::allTrue(def<_def0*_parent._reduction) || Simd::max(def)<1E-30);    // convergence check
        return _res.converged;
      }
 
    protected:
      void finalize(){
        _res.converged = (Simd::allTrue(_def<_def0*_parent._reduction) || Simd::max(_def)<1E-30);
        _res.iterations = _i;
        _res.reduction = static_cast<double>(Simd::max(_def/_def0));
        _res.conv_rate  = pow(_res.reduction,1.0/_i);
        _res.elapsed = _watch.elapsed();
        if (_parent._verbose>0)                 // final print
          {
            std::cout << "=== rate=" << _res.conv_rate
                      << ", T=" << _res.elapsed
                      << ", TIT=" << _res.elapsed/_res.iterations
                      << ", IT=" << _res.iterations << std::endl;
          }
      }
 
      real_type _def0 = 0.0, _def = 0.0;
      CountType _i;
      Timer _watch;
      InverseOperatorResult& _res;
      const IterativeSolver& _parent;
      bool _valid;
    };
 
  protected:
    std::shared_ptr<LinearOperator<X,Y>> _op;
    std::shared_ptr<Preconditioner<X,Y>> _prec;
    std::shared_ptr<ScalarProduct<X>> _sp;
    scalar_real_type _reduction;
    int _maxit;
    int _verbose;
    SolverCategory::Category _category;
  };
 
  template <typename ISTLLinearSolver, typename BCRSMatrix>
  class SolverHelper
  {
  public:
    static void setMatrix (ISTLLinearSolver& solver,
                           const BCRSMatrix& matrix)
    {
      static const bool is_direct_solver
        = Dune::IsDirectSolver<ISTLLinearSolver>::value;
      SolverHelper<ISTLLinearSolver,BCRSMatrix>::
        Implementation<is_direct_solver>::setMatrix(solver,matrix);
    }
 
  protected:
    template <bool is_direct_solver, typename Dummy = void>
    struct Implementation
    {
      static void setMatrix (ISTLLinearSolver&,
                             const BCRSMatrix&)
      {}
    };
 
    template <typename Dummy>
    struct Implementation<true,Dummy>
    {
      static void setMatrix (ISTLLinearSolver& solver,
                             const BCRSMatrix& matrix)
      {
        solver.setMatrix(matrix);
      }
    };
  };
 
}
 
#endif