CHANGELOG.md

@@ -37,6 +37,9 @@
 
 - Deprecated support for SuperLU 4.x. It will be removed after Dune 2.7.
 
+- The interface methods `dot()` and `norm()` of ScalarProduct are now `const`. You will
+  have to adjust the method signatures in your own scalar product implementations.
+
 # Release 2.6
 
 - `BDMatrix` objects can now be constructed and assigned from `std::initializer_list`.

dune/istl/colcompmatrix.hh

@@ -450,7 +450,7 @@ namespace Dune
   {
     for(size_type i=0; i<n; i++) {
       for(size_type j=0; j<m; j++) {
-        assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
+        assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<(size_type)mat->colstart[colindex*m+j+1]);
         assert((size_type)marker[colindex*m+j]<mat->Nnz_);
         mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
         mat->values[marker[colindex*m+j]]=Impl::asMatrix(*col)[i][j];

dune/istl/paamg/amg.hh

@@ -123,6 +123,32 @@ namespace Dune
           const SmootherArgs& smootherArgs=SmootherArgs(),
           const ParallelInformation& pinfo=ParallelInformation());
 
+      /*!
+         \brief Constructor an AMG via ParameterTree.
+
+         \param fineOperator The operator on the fine level.
+         \param configuration ParameterTree containing AMG parameters.
+         \param pinfo Optionally, specify ParallelInformation
+
+         ParameterTree Key         | Meaning
+         --------------------------|------------
+         smootherIterations        | The number of iterations to perform.
+         smootherRelaxation        | Common parameter defined [here](@ref ISTL_Factory_Common_Params).
+         maxLevel                  | Maximum number of levels allowed in the hierarchy.
+         coarsenTarget             | Maximum number of unknowns on the coarsest level.
+         prolongationDampingFactor | Damping factor for the prolongation.
+         alpha                     | Scaling avlue for marking connections as strong.
+         beta                      | Treshold for marking nodes as isolated.
+         additive                  | Whether to use additive multigrid.
+         gamma                     | 1 for V-cycle, 2 for W-cycle.
+         preSteps                  | Number of presmoothing steps.
+         postSteps                 | Number of postsmoothing steps.
+         verbosity                 | Output verbosity.
+
+         See \ref ISTL_Factory for the ParameterTree layout and examples.
+       */
+      AMG(std::shared_ptr<const Operator> fineOperator, const ParameterTree& configuration, const ParallelInformation& pinfo=ParallelInformation());
+
       /**
        * @brief Copy constructor.
        */
@@ -354,6 +380,60 @@ namespace Dune
       createHierarchies(criterion, matrixptr, pinfo);
     }
 
+    template<class M, class X, class S, class PI, class A>
+    AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
+                         const ParameterTree& configuration,
+                         const ParallelInformation& pinfo) :
+      smoothers_(new Hierarchy<Smoother,A>),
+      solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
+      coarsesolverconverged(true), coarseSmoother_(),
+      category_(SolverCategory::category(pinfo))
+    {
+
+      if (configuration.hasKey ("smootherIterations"))
+        smootherArgs_.iterations = configuration.get<int>("smootherIterations");
+
+      if (configuration.hasKey ("smootherRelaxation"))
+        smootherArgs_.relaxationFactor = configuration.get<typename SmootherArgs::RelaxationFactor>("smootherRelaxation");
+
+      typedef Amg::RowSum Norm;
+      typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<typename M::matrix_type,Norm> > Criterion;
+
+      Criterion criterion (configuration.get<int>("maxLevel"));
+
+      if (configuration.hasKey ("coarsenTarget"))
+        criterion.setCoarsenTarget (configuration.get<int>("coarsenTarget"));
+      if (configuration.hasKey ("prolongationDampingFactor"))
+        criterion.setProlongationDampingFactor (configuration.get<double>("prolongationDampingFactor"));
+
+      if (configuration.hasKey ("alpha"))
+        criterion.setAlpha (configuration.get<double> ("alpha"));
+
+      if (configuration.hasKey ("beta"))
+        criterion.setBeta (configuration.get<double> ("beta"));
+
+      if (configuration.hasKey ("gamma"))
+        criterion.setGamma (configuration.get<std::size_t> ("gamma"));
+      gamma_ = criterion.getGamma();
+
+      if (configuration.hasKey ("additive"))
+        criterion.setAdditive (configuration.get<bool>("additive"));
+      additive = criterion.getAdditive();
+
+      if (configuration.hasKey ("preSteps"))
+        criterion.setNoPreSmoothSteps (configuration.get<std::size_t> ("preSteps"));
+      preSteps_ = criterion.getNoPreSmoothSteps ();
+
+      if (configuration.hasKey ("postSteps"))
+        criterion.setNoPostSmoothSteps (configuration.get<std::size_t> ("preSteps"));
+      postSteps_ = criterion.getNoPostSmoothSteps ();
+
+      verbosity_ = configuration.get("verbosity",2);
+      criterion.setDebugLevel (verbosity_);
+
+      createHierarchies(criterion, matrixptr, pinfo);
+    }
+
     template <class Matrix,
               class Vector>
     struct DirectSolverSelector
@@ -423,7 +503,7 @@ namespace Dune
     template<class C>
     void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
       const std::shared_ptr<const Operator>& matrixptr,
-      const PI& pinfo)
+                                            const PI& pinfo)
     {
       Timer watch;
       matrices_ = std::make_shared<OperatorHierarchy>(
@@ -550,7 +630,7 @@ namespace Dune
       for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
         bool isDirichlet = true;
         bool hasDiagonal = false;
-        Block diagonal;
+        Block diagonal{};
         for(ColIter col=row->begin(); col!=row->end(); ++col) {
           if(row.index()==col.index()) {
             diagonal = *col;

dune/istl/paamg/smoother.hh

@@ -67,19 +67,22 @@ namespace Dune
 
     };
 
-    template<class X, class Y, class C, class T>
-    struct SmootherTraits<BlockPreconditioner<X,Y,C,T> >
+    template<class X, class Y>
+    struct SmootherTraits<Richardson<X,Y>>
     {
-      typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
+      typedef DefaultSmootherArgs<typename X::field_type> Arguments;
 
     };
 
+    template<class X, class Y, class C, class T>
+    struct SmootherTraits<BlockPreconditioner<X,Y,C,T> >
+        : public SmootherTraits<T>
+    {};
+
     template<class C, class T>
     struct SmootherTraits<NonoverlappingBlockPreconditioner<C,T> >
-    {
-      typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
-
-    };
+        : public SmootherTraits<T>
+    {};
 
     /**
      * @brief Construction Arguments for the default smoothers
@@ -168,6 +171,49 @@ namespace Dune
     };
 
 
+    template<class X, class Y>
+    class DefaultConstructionArgs<Richardson<X,Y>>
+    {
+      typedef Richardson<X,Y> T;
+
+      typedef typename SmootherTraits<T>::Arguments SmootherArgs;
+
+    public:
+      virtual ~DefaultConstructionArgs()
+      {}
+
+      template <class... Args>
+      void setMatrix(const Args&...)
+      {}
+
+      void setArgs(const SmootherArgs& args)
+      {
+        args_=&args;
+      }
+
+      template<class T1>
+      void setComm(T1& comm)
+      {
+        DUNE_UNUSED_PARAMETER(comm);
+      }
+
+      const SequentialInformation& getComm()
+      {
+        return comm_;
+      }
+
+      const SmootherArgs getArgs() const
+      {
+        return *args_;
+      }
+
+    private:
+      const SmootherArgs* args_;
+      SequentialInformation comm_;
+    };
+
+
+
     template<class T>
     class ConstructionTraits;
 
@@ -233,6 +279,21 @@ namespace Dune
       }
     };
 
+    /**
+     * @brief Policy for the construction of the Richardson smoother
+     */
+    template<class X, class Y>
+    struct ConstructionTraits<Richardson<X,Y> >
+    {
+      typedef DefaultConstructionArgs<Richardson<X,Y> > Arguments;
+
+      static inline std::shared_ptr<Richardson<X,Y>> construct(Arguments& args)
+      {
+        return std::make_shared<Richardson<X,Y>>
+          (args.getArgs().relaxationFactor);
+      }
+    };
+
 
 DUNE_NO_DEPRECATED_BEGIN // for deprecated SeqILU0
     /**

dune/istl/preconditioners.hh

@@ -12,6 +12,7 @@
 
 #include <dune/common/simd/simd.hh>
 #include <dune/common/unused.hh>
+#include <dune/common/parametertree.hh>
 
 #include "preconditioner.hh"
 #include "solver.hh"
@@ -159,6 +160,27 @@ namespace Dune {
       CheckIfDiagonalPresent<M,l>::check(_A_);
     }
 
+    /*!
+       \brief Constructor.
+
+       \param A The matrix to operate on.
+       \param configuration ParameterTree containing preconditioner parameters.
+
+       ParameterTree Key | Meaning
+       ------------------|------------
+       iterations        | The number of iterations to perform.
+       relaxation        | Common parameter defined [here](@ref ISTL_Factory_Common_Params).
+
+       See \ref ISTL_Factory for the ParameterTree layout and examples.
+     */
+    SeqSSOR (const M& A, const ParameterTree& configuration)
+      : _A_(A)
+    {
+      _n = configuration.get<field_type>("iterations");
+      _w = configuration.get<field_type>("relaxation");
+      CheckIfDiagonalPresent<M,l>::check(_A_);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -249,6 +271,17 @@ namespace Dune {
       CheckIfDiagonalPresent<M,l>::check(_A_);
     }
 
+    /*!
+       \copydoc SeqSSOR::SeqSSOR(const M&,const ParameterTree&)
+     */
+    SeqSOR (const M& A, const ParameterTree& configuration)
+      : _A_(A)
+    {
+      _n = configuration.get<field_type>("iterations");
+      _w = configuration.get<field_type>("relaxation");
+      CheckIfDiagonalPresent<M,l>::check(_A_);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -354,6 +387,17 @@ namespace Dune {
       CheckIfDiagonalPresent<M,l>::check(_A_);
     }
 
+    /*!
+       \copydoc SeqSSOR::SeqSSOR(const M&,const ParameterTree&)
+     */
+    SeqGS (const M& A, const ParameterTree& configuration)
+      : _A_(A)
+    {
+      _n = configuration.get<field_type>("iterations");
+      _w = configuration.get<field_type>("relaxation");
+      CheckIfDiagonalPresent<M,l>::check(_A_);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -440,6 +484,17 @@ namespace Dune {
       CheckIfDiagonalPresent<M,l>::check(_A_);
     }
 
+    /*!
+       \copydoc SeqSSOR::SeqSSOR(const M&,const ParameterTree&)
+     */
+    SeqJac (const M& A, const ParameterTree& configuration)
+      : _A_(A)
+    {
+      _n = configuration.get<field_type>("iterations");
+      _w = configuration.get<field_type>("relaxation");
+      CheckIfDiagonalPresent<M,l>::check(_A_);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -538,7 +593,7 @@ namespace Dune {
 
        Constructor invoking ILU(n).
        \param A The matrix to operate on.
-       \param n The number of iterations to perform.
+       \param n The order of the ILU decomposition.
        \param w The relaxation factor.
        \param resort true if a resort of the computed ILU for improved performance should be done.
      */
@@ -679,6 +734,25 @@ namespace Dune {
       bilu0_decomposition(ILU);
     }
 
+    /*!
+       \brief Constructor.
+
+       \param A The matrix to operate on.
+       \param configuration ParameterTree containing preconditioner parameters.
+
+       ParameterTree Key | Meaning
+       ------------------|------------
+       relaxation        | Common parameter defined [here](@ref ISTL_Factory_Common_Params).
+
+       See \ref ISTL_Factory for the ParameterTree layout and examples.
+     */
+    SeqILU0 (const M& A, const ParameterTree& configuration)
+      : ILU(A) // copy A
+    {
+      _w = configuration.get<field_type>("relaxation");
+      bilu0_decomposition(ILU);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -758,7 +832,7 @@ namespace Dune {
 
        Constructor gets all parameters to operate the prec.
        \param A The matrix to operate on.
-       \param n The number of iterations to perform.
+       \param n The order of the ILU decomposition.
        \param w The relaxation factor.
      */
     SeqILUn (const M& A, int n, scalar_field_type w)
@@ -769,6 +843,17 @@ namespace Dune {
       bilu_decomposition(A,n,ILU);
     }
 
+    /*!
+       \copydoc SeqSSOR::SeqSSOR(const M&,const ParameterTree&)
+     */
+    SeqILUn (const M& A, const ParameterTree& configuration)
+      : ILU(A.N(),A.M(),M::row_wise)
+    {
+      _n = configuration.get<field_type>("iterations");
+      _w = configuration.get<field_type>("relaxation");
+      bilu_decomposition(A,_n,ILU);
+    }
+
     /*!
        \brief Prepare the preconditioner.
 
@@ -847,6 +932,14 @@ namespace Dune {
       _w(w)
     {}
 
+    /*!
+       \copydoc SeqILU0::SeqILU0(const M&,const ParameterTree&)
+     */
+    Richardson (const ParameterTree& configuration)
+    {
+      _w = configuration.get<field_type>("relaxation");
+    }
+
     /*!
        \brief Prepare the preconditioner.
 

dune/istl/solver.hh

@@ -6,11 +6,15 @@
 
 #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"
@@ -193,7 +197,7 @@ namespace Dune
      constructors, which are then only imported in the actual solver
      implementation.
    */
-  template<class X, class Y>
+  template<class X, class Y, class CountType = unsigned int>
   class IterativeSolver : public InverseOperator<X,Y>{
   public:
     using typename InverseOperator<X,Y>::domain_type;
@@ -266,6 +270,56 @@ namespace Dune
         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,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<real_type>("reduction"),
+        configuration.get<int>("maxit"),
+        configuration.get<int>("verbose"))
+    {}
+
+    /**
+        \brief General constructor to initialize an iterative solver
+
+        \param op The operator we solve.
+        \param sp The scalar product to use, e. g. SeqScalarproduct.
+        \param prec The preconditioner to apply in each iteration of the loop.
+        Has to inherit from Preconditioner.
+        \param reduction The relative defect reduction to achieve when applying
+        the operator.
+        \param maxit The maximum number of iteration steps allowed when applying
+        the operator.
+        \param verbose The verbosity level.
+
+        Verbose levels are:
+        <ul>
+        <li> 0 : print nothing </li>
+        <li> 1 : print initial and final defect and statistics </li>
+        <li> 2 : print line for each iteration </li>
+        </ul>
+     */
+    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)
@@ -298,6 +352,120 @@ namespace Dune
       return _category;
     }
 
+    std::string name() const{
+      std::string name = className(*this);
+      return name.substr(0, name.find("<"));
+    }
+
+      /*!
+     \brief Class for controlling iterative methods
+
+     This class provides building blocks for a iterative method. It does all
+     things that have to do with output, residual checking (NaN, infinite,
+     convergence) and sets also the fields of InverseOperatorResult.
+
+     Instances of this class are meant to create with
+     IterativeSolver::startIteration and stored as a local variable in the apply
+     method. If the scope of the apply method is left the destructor of this
+     class sets all the solver statistics in the InverseOperatorResult and
+     prints the final output.
+
+     During the iteration in every step Iteration::step should be called with
+     the current iteration count and norm of the residual. It returns true if
+     convergence is achieved.
+   */
+    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();
+      }
+
+      /*! \brief registers the iteration step, checks for invalid defect norm
+          and convergence.
+
+        \param i The current iteration count
+        \param def The current norm of the defect
+
+        \return true is convergence is achieved
+
+        \throw SolverAbort when `def` contains inf or NaN
+       */
+      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;