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 <dune/common/exceptions.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/simd.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;
    }
 
    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 SimdScalar<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) << norm << " ";
      s << std::setw(normSpacing) << 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) << 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!");
    }
 
    // #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;
    }
 
  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