"Fossies" - the Fresh Open Source Software Archive

Member "dune-istl-2.7.1/dune/istl/solver.hh" (26 Nov 2020, 19624 Bytes) of package /linux/misc/dune/dune-istl-2.7.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "solver.hh" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2.6.0_vs_2.7.0.

    1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
    2 // vi: set et ts=4 sw=2 sts=2:
    3 
    4 #ifndef DUNE_ISTL_SOLVER_HH
    5 #define DUNE_ISTL_SOLVER_HH
    6 
    7 #include <iomanip>
    8 #include <ostream>
    9 #include <string>
   10 #include <functional>
   11 
   12 #include <dune/common/exceptions.hh>
   13 #include <dune/common/shared_ptr.hh>
   14 #include <dune/common/simd/io.hh>
   15 #include <dune/common/simd/simd.hh>
   16 #include <dune/common/parametertree.hh>
   17 #include <dune/common/timer.hh>
   18 
   19 #include "solvertype.hh"
   20 #include "preconditioner.hh"
   21 #include "operators.hh"
   22 #include "scalarproducts.hh"
   23 
   24 namespace Dune
   25 {
   26 /**
   27  * @addtogroup ISTL_Solvers
   28  * @{
   29  */
   30 /** \file
   31 
   32       \brief   Define general, extensible interface for inverse operators.
   33 
   34       Implementation here covers only inversion of linear operators,
   35       but the implementation might be used for nonlinear operators
   36       as well.
   37    */
   38   /**
   39       \brief Statistics about the application of an inverse operator
   40 
   41       The return value of an application of the inverse
   42       operator delivers some important information about
   43       the iteration.
   44    */
   45   struct InverseOperatorResult
   46   {
   47     /** \brief Default constructor */
   48     InverseOperatorResult ()
   49     {
   50       clear();
   51     }
   52 
   53     /** \brief Resets all data */
   54     void clear ()
   55     {
   56       iterations = 0;
   57       reduction = 0;
   58       converged = false;
   59       conv_rate = 1;
   60       elapsed = 0;
   61       condition_estimate = -1;
   62     }
   63 
   64     /** \brief Number of iterations */
   65     int iterations;
   66 
   67     /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */
   68     double reduction;
   69 
   70     /** \brief True if convergence criterion has been met */
   71     bool converged;
   72 
   73     /** \brief Convergence rate (average reduction per step) */
   74     double conv_rate;
   75 
   76     /** \brief Estimate of condition number */
   77     double condition_estimate = -1;
   78 
   79     /** \brief Elapsed time in seconds */
   80     double elapsed;
   81   };
   82 
   83 
   84   //=====================================================================
   85   /*!
   86      \brief Abstract base class for all solvers.
   87 
   88      An InverseOperator computes the solution of \f$ A(x)=b\f$ where
   89      \f$ A : X \to Y \f$ is an operator.
   90      Note that the solver "knows" which operator
   91      to invert and which preconditioner to apply (if any). The
   92      user is only interested in inverting the operator.
   93      InverseOperator might be a Newton scheme, a Krylov subspace method,
   94      or a direct solver or just anything.
   95    */
   96   template<class X, class Y>
   97   class InverseOperator {
   98   public:
   99     //! \brief Type of the domain of the operator to be inverted.
  100     typedef X domain_type;
  101 
  102     //! \brief Type of the range of the operator to be inverted.
  103     typedef Y range_type;
  104 
  105     /** \brief The field type of the operator. */
  106     typedef typename X::field_type field_type;
  107 
  108     //! \brief The real type of the field type (is the same if using real numbers, but differs for std::complex)
  109     typedef typename FieldTraits<field_type>::real_type real_type;
  110 
  111     //! \brief scalar type underlying the field_type
  112     typedef Simd::Scalar<real_type> scalar_real_type;
  113 
  114     /**
  115         \brief Apply inverse operator,
  116 
  117         \warning Note: right hand side b may be overwritten!
  118 
  119         \param x The left hand side to store the result in.
  120         \param b The right hand side
  121         \param res Object to store the statistics about applying the operator.
  122 
  123         \throw SolverAbort When the solver detects a problem and cannot
  124                            continue
  125      */
  126     virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
  127 
  128     /*!
  129        \brief apply inverse operator, with given convergence criteria.
  130 
  131        \warning Right hand side b may be overwritten!
  132 
  133        \param x The left hand side to store the result in.
  134        \param b The right hand side
  135        \param reduction The minimum defect reduction to achieve.
  136        \param res Object to store the statistics about applying the operator.
  137 
  138        \throw SolverAbort When the solver detects a problem and cannot
  139                           continue
  140      */
  141     virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
  142 
  143     //! Category of the solver (see SolverCategory::Category)
  144     virtual SolverCategory::Category category() const
  145 #ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
  146     {
  147       DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
  148     }
  149 #else
  150     = 0;
  151 #endif
  152 
  153     //! \brief Destructor
  154     virtual ~InverseOperator () {}
  155 
  156   protected:
  157     // spacing values
  158     enum { iterationSpacing = 5 , normSpacing = 16 };
  159 
  160     //! helper function for printing header of solver output
  161     void printHeader(std::ostream& s) const
  162     {
  163       s << std::setw(iterationSpacing)  << " Iter";
  164       s << std::setw(normSpacing) << "Defect";
  165       s << std::setw(normSpacing) << "Rate" << std::endl;
  166     }
  167 
  168     //! helper function for printing solver output
  169     template <typename CountType, typename DataType>
  170     void printOutput(std::ostream& s,
  171                      const CountType& iter,
  172                      const DataType& norm,
  173                      const DataType& norm_old) const
  174     {
  175       const DataType rate = norm/norm_old;
  176       s << std::setw(iterationSpacing)  << iter << " ";
  177       s << std::setw(normSpacing) << Simd::io(norm) << " ";
  178       s << std::setw(normSpacing) << Simd::io(rate) << std::endl;
  179     }
  180 
  181     //! helper function for printing solver output
  182     template <typename CountType, typename DataType>
  183     void printOutput(std::ostream& s,
  184                      const CountType& iter,
  185                      const DataType& norm) const
  186     {
  187       s << std::setw(iterationSpacing)  << iter << " ";
  188       s << std::setw(normSpacing) << Simd::io(norm) << std::endl;
  189     }
  190   };
  191 
  192   /*!
  193      \brief Base class for all implementations of iterative solvers
  194 
  195      This class provides all storage, which is needed by the usual
  196      iterative solvers. In additional it provides all the necessary
  197      constructors, which are then only imported in the actual solver
  198      implementation.
  199    */
  200   template<class X, class Y>
  201   class IterativeSolver : public InverseOperator<X,Y>{
  202   public:
  203     using typename InverseOperator<X,Y>::domain_type;
  204     using typename InverseOperator<X,Y>::range_type;
  205     using typename InverseOperator<X,Y>::field_type;
  206     using typename InverseOperator<X,Y>::real_type;
  207     using typename InverseOperator<X,Y>::scalar_real_type;
  208 
  209     /*!
  210        \brief General constructor to initialize an iterative solver.
  211 
  212        \param op The operator we solve.
  213        \param prec The preconditioner to apply in each iteration of the loop.
  214        Has to inherit from Preconditioner.
  215        \param reduction The relative defect reduction to achieve when applying
  216        the operator.
  217        \param maxit The maximum number of iteration steps allowed when applying
  218        the operator.
  219        \param verbose The verbosity level.
  220 
  221        Verbose levels are:
  222        <ul>
  223        <li> 0 : print nothing </li>
  224        <li> 1 : print initial and final defect and statistics </li>
  225        <li> 2 : print line for each iteration </li>
  226        </ul>
  227      */
  228     IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
  229       _op(stackobject_to_shared_ptr(op)),
  230       _prec(stackobject_to_shared_ptr(prec)),
  231       _sp(new SeqScalarProduct<X>),
  232       _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
  233     {
  234       if(SolverCategory::category(op) != SolverCategory::sequential)
  235         DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
  236       if(SolverCategory::category(prec) != SolverCategory::sequential)
  237         DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
  238     }
  239 
  240     /**
  241         \brief General constructor to initialize an iterative solver
  242 
  243         \param op The operator we solve.
  244         \param sp The scalar product to use, e. g. SeqScalarproduct.
  245         \param prec The preconditioner to apply in each iteration of the loop.
  246         Has to inherit from Preconditioner.
  247         \param reduction The relative defect reduction to achieve when applying
  248         the operator.
  249         \param maxit The maximum number of iteration steps allowed when applying
  250         the operator.
  251         \param verbose The verbosity level.
  252 
  253         Verbose levels are:
  254         <ul>
  255         <li> 0 : print nothing </li>
  256         <li> 1 : print initial and final defect and statistics </li>
  257         <li> 2 : print line for each iteration </li>
  258         </ul>
  259      */
  260     IterativeSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
  261       scalar_real_type reduction, int maxit, int verbose) :
  262       _op(stackobject_to_shared_ptr(op)),
  263       _prec(stackobject_to_shared_ptr(prec)),
  264       _sp(stackobject_to_shared_ptr(sp)),
  265       _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
  266     {
  267       if(SolverCategory::category(op) != SolverCategory::category(prec))
  268         DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
  269       if(SolverCategory::category(op) != SolverCategory::category(sp))
  270         DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
  271     }
  272 
  273         /*!
  274        \brief Constructor.
  275 
  276        \param op The operator we solve
  277        \param prec The preconditioner to apply in each iteration of the loop.
  278        \param configuration ParameterTree containing iterative solver parameters.
  279 
  280        ParameterTree Key | Meaning
  281        ------------------|------------
  282        reduction         | The relative defect reduction to achieve when applying the operator
  283        maxit             | The maximum number of iteration steps allowed when applying the operator
  284        verbose           | The verbosity level
  285 
  286        See \ref ISTL_Factory for the ParameterTree layout and examples.
  287      */
  288     IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
  289       IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
  290         configuration.get<real_type>("reduction"),
  291         configuration.get<int>("maxit"),
  292         configuration.get<int>("verbose"))
  293     {}
  294 
  295     /*!
  296        \brief Constructor.
  297 
  298        \param op The operator we solve
  299        \param sp The scalar product to use, e. g. SeqScalarproduct.
  300        \param prec The preconditioner to apply in each iteration of the loop.
  301        \param configuration ParameterTree containing iterative solver parameters.
  302 
  303        ParameterTree Key | Meaning
  304        ------------------|------------
  305        reduction         | The relative defect reduction to achieve when applying the operator
  306        maxit             | The maximum number of iteration steps allowed when applying the operator
  307        verbose           | The verbosity level
  308 
  309        See \ref ISTL_Factory for the ParameterTree layout and examples.
  310      */
  311     IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
  312       IterativeSolver(op,sp,prec,
  313         configuration.get<scalar_real_type>("reduction"),
  314         configuration.get<int>("maxit"),
  315         configuration.get<int>("verbose"))
  316     {}
  317 
  318     /**
  319         \brief General constructor to initialize an iterative solver
  320 
  321         \param op The operator we solve.
  322         \param sp The scalar product to use, e. g. SeqScalarproduct.
  323         \param prec The preconditioner to apply in each iteration of the loop.
  324         Has to inherit from Preconditioner.
  325         \param reduction The relative defect reduction to achieve when applying
  326         the operator.
  327         \param maxit The maximum number of iteration steps allowed when applying
  328         the operator.
  329         \param verbose The verbosity level.
  330 
  331         Verbose levels are:
  332         <ul>
  333         <li> 0 : print nothing </li>
  334         <li> 1 : print initial and final defect and statistics </li>
  335         <li> 2 : print line for each iteration </li>
  336         </ul>
  337      */
  338     IterativeSolver (std::shared_ptr<LinearOperator<X,Y>> op,
  339                      std::shared_ptr<ScalarProduct<X>> sp,
  340                      std::shared_ptr<Preconditioner<X,Y>> prec,
  341                      scalar_real_type reduction, int maxit, int verbose) :
  342       _op(op),
  343       _prec(prec),
  344       _sp(sp),
  345       _reduction(reduction), _maxit(maxit), _verbose(verbose),
  346       _category(SolverCategory::category(*op))
  347     {
  348       if(SolverCategory::category(*op) != SolverCategory::category(*prec))
  349         DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
  350       if(SolverCategory::category(*op) != SolverCategory::category(*sp))
  351         DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
  352     }
  353 
  354     // #warning actually we want to have this as the default and just implement the second one
  355     // //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
  356     // virtual void apply (X& x, Y& b, InverseOperatorResult& res)
  357     // {
  358     //   apply(x,b,_reduction,res);
  359     // }
  360 
  361 #ifndef DOXYGEN
  362     // make sure the three-argument apply from the base class does not get shadowed
  363     // by the redefined four-argument version below
  364     using InverseOperator<X,Y>::apply;
  365 #endif
  366 
  367     /*!
  368        \brief Apply inverse operator with given reduction factor.
  369 
  370        \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
  371      */
  372     virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
  373     {
  374       scalar_real_type saved_reduction = _reduction;
  375       _reduction = reduction;
  376       this->apply(x,b,res);
  377       _reduction = saved_reduction;
  378     }
  379 
  380     //! Category of the solver (see SolverCategory::Category)
  381     virtual SolverCategory::Category category() const
  382     {
  383       return _category;
  384     }
  385 
  386     std::string name() const{
  387       std::string name = className(*this);
  388       return name.substr(0, name.find("<"));
  389     }
  390 
  391       /*!
  392      \brief Class for controlling iterative methods
  393 
  394      This class provides building blocks for a iterative method. It does all
  395      things that have to do with output, residual checking (NaN, infinite,
  396      convergence) and sets also the fields of InverseOperatorResult.
  397 
  398      Instances of this class are meant to create with
  399      IterativeSolver::startIteration and stored as a local variable in the apply
  400      method. If the scope of the apply method is left the destructor of this
  401      class sets all the solver statistics in the InverseOperatorResult and
  402      prints the final output.
  403 
  404      During the iteration in every step Iteration::step should be called with
  405      the current iteration count and norm of the residual. It returns true if
  406      convergence is achieved.
  407    */
  408     template<class CountType = unsigned int>
  409     class Iteration {
  410     public:
  411       Iteration(const IterativeSolver& parent, InverseOperatorResult& res)
  412         : _i(0)
  413         , _res(res)
  414         , _parent(parent)
  415         , _valid(true)
  416       {
  417         res.clear();
  418         if(_parent._verbose>0){
  419           std::cout << "=== " << parent.name() << std::endl;
  420           if(_parent._verbose > 1)
  421             _parent.printHeader(std::cout);
  422         }
  423       }
  424 
  425       Iteration(const Iteration&) = delete;
  426       Iteration(Iteration&& other)
  427         : _def0(other._def0)
  428         , _def(other._def)
  429         , _i(other._i)
  430         , _watch(other._watch)
  431         , _res(other._res)
  432         , _parent(other._parent)
  433         , _valid(other._valid)
  434       {
  435         other._valid = false;
  436       }
  437 
  438       ~Iteration(){
  439         if(_valid)
  440           finalize();
  441       }
  442 
  443       /*! \brief registers the iteration step, checks for invalid defect norm
  444           and convergence.
  445 
  446         \param i The current iteration count
  447         \param def The current norm of the defect
  448 
  449         \return true is convergence is achieved
  450 
  451         \throw SolverAbort when `def` contains inf or NaN
  452        */
  453       bool step(CountType i, real_type def){
  454         if (!Simd::allTrue(isFinite(def))) // check for inf or NaN
  455         {
  456           if (_parent._verbose>0)
  457             std::cout << "=== " << _parent.name() << ": abort due to infinite or NaN defect"
  458                       << std::endl;
  459           DUNE_THROW(SolverAbort,
  460                      _parent.name() << ": defect=" << Simd::io(def)
  461                      << " is infinite or NaN");
  462         }
  463         if(i == 0)
  464           _def0 = def;
  465         if(_parent._verbose > 1){
  466           if(i!=0)
  467             _parent.printOutput(std::cout,i,def,_def);
  468           else
  469             _parent.printOutput(std::cout,i,def);
  470         }
  471         _def = def;
  472         _i = i;
  473         _res.converged = (Simd::allTrue(def<_def0*_parent._reduction) || Simd::max(def)<1E-30);    // convergence check
  474         return _res.converged;
  475       }
  476 
  477     protected:
  478       void finalize(){
  479         _res.converged = (Simd::allTrue(_def<_def0*_parent._reduction) || Simd::max(_def)<1E-30);
  480         _res.iterations = _i;
  481         _res.reduction = static_cast<double>(Simd::max(_def/_def0));
  482         _res.conv_rate  = pow(_res.reduction,1.0/_i);
  483         _res.elapsed = _watch.elapsed();
  484         if (_parent._verbose>0)                 // final print
  485           {
  486             std::cout << "=== rate=" << _res.conv_rate
  487                       << ", T=" << _res.elapsed
  488                       << ", TIT=" << _res.elapsed/_res.iterations
  489                       << ", IT=" << _res.iterations << std::endl;
  490           }
  491       }
  492 
  493       real_type _def0 = 0.0, _def = 0.0;
  494       CountType _i;
  495       Timer _watch;
  496       InverseOperatorResult& _res;
  497       const IterativeSolver& _parent;
  498       bool _valid;
  499     };
  500 
  501   protected:
  502     std::shared_ptr<LinearOperator<X,Y>> _op;
  503     std::shared_ptr<Preconditioner<X,Y>> _prec;
  504     std::shared_ptr<ScalarProduct<X>> _sp;
  505     scalar_real_type _reduction;
  506     int _maxit;
  507     int _verbose;
  508     SolverCategory::Category _category;
  509   };
  510 
  511   /**
  512    * \brief Helper class for notifying a DUNE-ISTL linear solver about
  513    *        a change of the iteration matrix object in a unified way,
  514    *        i.e. independent from the solver's type (direct/iterative).
  515    *
  516    * \author Sebastian Westerheide.
  517    */
  518   template <typename ISTLLinearSolver, typename BCRSMatrix>
  519   class SolverHelper
  520   {
  521   public:
  522     static void setMatrix (ISTLLinearSolver& solver,
  523                            const BCRSMatrix& matrix)
  524     {
  525       static const bool is_direct_solver
  526         = Dune::IsDirectSolver<ISTLLinearSolver>::value;
  527       SolverHelper<ISTLLinearSolver,BCRSMatrix>::
  528         Implementation<is_direct_solver>::setMatrix(solver,matrix);
  529     }
  530 
  531   protected:
  532     /**
  533      * \brief Implementation that works together with iterative ISTL
  534      *        solvers, e.g. Dune::CGSolver or Dune::BiCGSTABSolver.
  535      */
  536     template <bool is_direct_solver, typename Dummy = void>
  537     struct Implementation
  538     {
  539       static void setMatrix (ISTLLinearSolver&,
  540                              const BCRSMatrix&)
  541       {}
  542     };
  543 
  544     /**
  545      * \brief Implementation that works together with direct ISTL
  546      *        solvers, e.g. Dune::SuperLU or Dune::UMFPack.
  547      */
  548     template <typename Dummy>
  549     struct Implementation<true,Dummy>
  550     {
  551       static void setMatrix (ISTLLinearSolver& solver,
  552                              const BCRSMatrix& matrix)
  553       {
  554         solver.setMatrix(matrix);
  555       }
  556     };
  557   };
  558 
  559 /**
  560  * @}
  561  */
  562 }
  563 
  564 #endif