istl/novlpistlsolverbackend.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_NOVLPISTLSOLVERBACKEND_HH
#define DUNE_NOVLPISTLSOLVERBACKEND_HH

#include <cstddef>

#include <dune/common/deprecated.hh>
#include <dune/common/parallel/mpihelper.hh>

#include <dune/grid/common/gridenums.hh>

#include <dune/istl/io.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/superlu.hh>

#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
#include <dune/pdelab/backend/istl/vector.hh>
#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
#include <dune/pdelab/backend/istl/blockmatrixdiagonal.hh>
#include <dune/pdelab/backend/istl/parallelhelper.hh>
#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>

namespace Dune {
  namespace PDELab {


    //========================================================
    // Generic support for nonoverlapping grids
    //========================================================


    template<typename GFS, typename M, typename X, typename Y>
    class NonoverlappingOperator
      : public Dune::AssembledLinearOperator<M,X,Y>
    {
    public:
      using matrix_type = Backend::Native<M>;
      using domain_type = Backend::Native<X>;
      using range_type = Backend::Native<Y>;
      typedef typename X::field_type field_type;

      //redefine the category, that is the only difference
      enum {category=Dune::SolverCategory::nonoverlapping};


      NonoverlappingOperator (const GFS& gfs_, const M& A)
        : gfs(gfs_), _A_(A)
      { }


      virtual void apply (const X& x, Y& y) const
      {
        using Backend::native;
        // apply local operator; now we have sum y_p = sequential y
        native(_A_).mv(native(x),native(y));

        // accumulate y on border
        Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
        if (gfs.gridView().comm().size()>1)
          gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
      }


      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
      {
        using Backend::native;
        // apply local operator; now we have sum y_p = sequential y
        native(_A_).usmv(alpha,native(x),native(y));

        // accumulate y on border
        Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
        if (gfs.gridView().comm().size()>1)
          gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
      }

      virtual const M& getmat () const
      {
        return _A_;
      }

    private:
      const GFS& gfs;
      const M& _A_;
    };

    // parallel scalar product assuming no overlap
    template<class GFS, class X>
    class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
    {
    public:
      typedef X domain_type;
      typedef typename X::ElementType field_type;

      enum {category=Dune::SolverCategory::nonoverlapping};

      NonoverlappingScalarProduct (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
        : gfs(gfs_), helper(helper_)
      {}

      virtual field_type dot (const X& x, const X& y)
      {
        // do local scalar product on unique partition
        field_type sum = helper.disjointDot(x,y);

        // do global communication
        return gfs.gridView().comm().sum(sum);
      }

      virtual double norm (const X& x)
      {
        return sqrt(static_cast<double>(this->dot(x,x)));
      }

      void make_consistent (X& x) const
      {
        Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,x);
        if (gfs.gridView().comm().size()>1)
          gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
      }

    private:
      const GFS& gfs;
      const istl::ParallelHelper<GFS>& helper;
    };

    // parallel Richardson preconditioner
    template<class GFS, class X, class Y>
    class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
    {
    public:
      typedef X domain_type;
      typedef Y range_type;
      typedef typename X::ElementType field_type;

      // define the category
      enum {
        category=Dune::SolverCategory::nonoverlapping
      };

      NonoverlappingRichardson (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
        : gfs(gfs_), helper(helper_)
      {
      }

      virtual void pre (X& x, Y& b) {}

      virtual void apply (X& v, const Y& d)
      {
        v = d;
      }

      virtual void post (X& x) {}

    private:
      const GFS& gfs;
      const istl::ParallelHelper<GFS>& helper;
    };


    template<typename A, typename X, typename Y>
    class NonoverlappingJacobi
      : public Dune::Preconditioner<X,Y>
    {

      typedef typename istl::BlockMatrixDiagonal<A>::MatrixElementVector Diagonal;

      Diagonal _inverse_diagonal;

    public:

      typedef X domain_type;

      typedef Y range_type;
      typedef typename X::ElementType field_type;

      enum {
        category=Dune::SolverCategory::nonoverlapping
      };


      template<typename GFS>
      NonoverlappingJacobi(const GFS& gfs, const A &m)
        : _inverse_diagonal(m)
      {
        // make the diagonal consistent...
        typename istl::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal);
        gfs.gridView().communicate(addDH,
                                   InteriorBorder_InteriorBorder_Interface,
                                   ForwardCommunication);

        // ... and then invert it
        _inverse_diagonal.invert();

      }

      virtual void pre (X& x, Y& b) {}

      /*
       * For this preconditioner, this method works with both consistent and
       * inconsistent vectors: if d is consistent, v will be consistent, if d
       * is inconsistent, v will be inconsistent.
       */
      virtual void apply (X& v, const Y& d)
      {
        _inverse_diagonal.mv(d,v);
      }

      virtual void post (X& x) {}
    };


    template<class GFS>
    class ISTLBackend_NOVLP_CG_NOPREC
    {
      typedef istl::ParallelHelper<GFS> PHELPER;

    public:
      explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
                                            unsigned maxiter_=5000,
                                            int verbose_=1)
        : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
      {}

      template<class V>
      typename V::ElementType norm (const V& v) const
      {
        V x(v); // make a copy because it has to be made consistent
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
      {
        typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP;
        POP pop(gfs,A);
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
        PRICH prich(gfs,phelper);
        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
        Dune::InverseOperatorResult stat;
        solver.apply(z,r,stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const Dune::PDELab::LinearSolverResult<double>& result() const
      {
        return res;
      }

    private:
      const GFS& gfs;
      PHELPER phelper;
      Dune::PDELab::LinearSolverResult<double> res;
      unsigned maxiter;
      int verbose;
    };

    template<class GFS>
    class ISTLBackend_NOVLP_CG_Jacobi
    {
      typedef istl::ParallelHelper<GFS> PHELPER;

      const GFS& gfs;
      PHELPER phelper;
      LinearSolverResult<double> res;
      unsigned maxiter;
      int verbose;

    public:

      explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
                                           unsigned maxiter_ = 5000,
                                           int verbose_ = 1) :
        gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
      {}


      template<class V>
      typename V::ElementType norm (const V& v) const
      {
        V x(v); // make a copy because it has to be made consistent
        typedef NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }


      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
      {
        typedef NonoverlappingOperator<GFS,M,V,W> POP;
        POP pop(gfs,A);
        typedef NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);

        typedef NonoverlappingJacobi<M,V,W> PPre;
        PPre ppre(gfs,Backend::native(A));

        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
        InverseOperatorResult stat;
        solver.apply(z,r,stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const LinearSolverResult<double>& result() const
      { return res; }
    };

    template<class GFS>
    class ISTLBackend_NOVLP_BCGS_NOPREC
    {
      typedef istl::ParallelHelper<GFS> PHELPER;

    public:
      explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
        : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
      {}

      template<class V>
      typename V::ElementType norm (const V& v) const
      {
        V x(v); // make a copy because it has to be made consistent
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
      {
        typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP;
        POP pop(gfs,A);
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
        PRICH prich(gfs,phelper);
        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
        Dune::InverseOperatorResult stat;
        solver.apply(z,r,stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const Dune::PDELab::LinearSolverResult<double>& result() const
      {
        return res;
      }

    private:
      const GFS& gfs;
      PHELPER phelper;
      Dune::PDELab::LinearSolverResult<double> res;
      unsigned maxiter;
      int verbose;
    };


    template<class GFS>
    class ISTLBackend_NOVLP_BCGS_Jacobi
    {
      typedef istl::ParallelHelper<GFS> PHELPER;

    public:
      explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
        : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
      {}

      template<class V>
      typename V::ElementType norm (const V& v) const
      {
        V x(v); // make a copy because it has to be made consistent
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
      {
        typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP;
        POP pop(gfs,A);
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);

        typedef NonoverlappingJacobi<M,V,W> PPre;
        PPre ppre(gfs,A);

        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
        Dune::InverseOperatorResult stat;
        solver.apply(z,r,stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const Dune::PDELab::LinearSolverResult<double>& result() const
      {
        return res;
      }

    private:
      const GFS& gfs;
      PHELPER phelper;
      Dune::PDELab::LinearSolverResult<double> res;
      unsigned maxiter;
      int verbose;
    };

    template<typename GFS>
    class ISTLBackend_NOVLP_ExplicitDiagonal
    {
      typedef istl::ParallelHelper<GFS> PHELPER;

      const GFS& gfs;
      PHELPER phelper;
      Dune::PDELab::LinearSolverResult<double> res;

    public:
      explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
        : gfs(gfs_), phelper(gfs)
      {}

      template<class V>
      typename V::ElementType norm (const V& v) const
      {
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        V x(v); // make a copy because it has to be made consistent
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
      {
        Dune::SeqJac<M,V,W> jac(A,1,1.0);
        jac.pre(z,r);
        jac.apply(z,r);
        jac.post(z);
        if (gfs.gridView().comm().size()>1)
        {
          Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,z);
          gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
        }
        res.converged  = true;
        res.iterations = 1;
        res.elapsed    = 0.0;
        res.reduction  = reduction;
        res.conv_rate  = reduction; // pow(reduction,1.0/1)
      }

      const Dune::PDELab::LinearSolverResult<double>& result() const
      {
        return res;
      }
    };


    template<class GO,
             template<class,class,class,int> class Preconditioner,
             template<class> class Solver>
    class ISTLBackend_NOVLP_BASE_PREC
    {
      typedef typename GO::Traits::TrialGridFunctionSpace GFS;
      typedef istl::ParallelHelper<GFS> PHELPER;

    public:
      explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
        : _grid_operator(grid_operator)
        , gfs(grid_operator.trialGridFunctionSpace())
        , phelper(gfs,verbose_)
        , maxiter(maxiter_)
        , steps(steps_)
        , verbose(verbose_)
      {}

      template<class Vector>
      typename Vector::ElementType norm (const Vector& v) const
      {
        Vector x(v); // make a copy because it has to be made consistent
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,Vector> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      template<class M, class V, class W>
      void apply(M& A, V& z, W& r, typename V::ElementType reduction)
      {
        using MatrixType = Backend::Native<M>;
        MatrixType& mat = Backend::native(A);
        using VectorType = Backend::Native<W>;
#if HAVE_MPI
        typedef typename istl::CommSelector<96,Dune::MPIHelper::isFake>::type Comm;
        _grid_operator.make_consistent(A);
        Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
        phelper.createIndexSetAndProjectForAMG(mat, oocc);
        typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
        Smoother smoother(mat, steps, 1.0);
        typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP;
        PSP psp(oocc);
        typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
        Operator oop(mat,oocc);
        typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother;
        ParSmoother parsmoother(smoother, oocc);
#else
        typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
        ParSmoother parsmoother(mat, steps, 1.0);
        typedef Dune::SeqScalarProduct<VectorType> PSP;
        PSP psp;
        typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
        Operator oop(mat);
#endif
        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
        Dune::InverseOperatorResult stat;
        //make r consistent
        if (gfs.gridView().comm().size()>1){
          Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
          gfs.gridView().communicate(adddh,
                                     Dune::InteriorBorder_InteriorBorder_Interface,
                                     Dune::ForwardCommunication);
        }

        solver.apply(z,r,stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const Dune::PDELab::LinearSolverResult<double>& result() const
      {
        return res;
      }

    private:
      const GO& _grid_operator;
      const GFS& gfs;
      PHELPER phelper;
      Dune::PDELab::LinearSolverResult<double> res;
      unsigned maxiter;
      unsigned steps;
      int verbose;
    };


  template<class GO>
  class ISTLBackend_NOVLP_BCGS_SSORk
    : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
  {

  public:
    explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
                                           int steps_=5, int verbose_=1)
      : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_)
    {}
  };

  template<class GO>
  class ISTLBackend_NOVLP_CG_SSORk
    : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
  {

  public:
    explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
                                         int steps_=5, int verbose_=1)
      : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_)
    {}
  };

    template<class GO,int s, template<class,class,class,int> class Preconditioner,
             template<class> class Solver>
    class ISTLBackend_AMG_NOVLP : public LinearResultStorage
    {
      typedef typename GO::Traits::TrialGridFunctionSpace GFS;
      typedef typename istl::ParallelHelper<GFS> PHELPER;
      typedef typename GO::Traits::Jacobian M;
      typedef Backend::Native<M> MatrixType;
      typedef typename GO::Traits::Domain V;
      typedef Backend::Native<V> VectorType;
      typedef typename istl::CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
#if HAVE_MPI
      typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
      typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother;
      typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
#else
      typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
      typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
#endif
      typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
      typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
      typedef Dune::Amg::Parameters Parameters;

    public:
      ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000,
                            int verbose_=1, bool reuse_=false,
                            bool usesuperlu_=true)
        : _grid_operator(grid_operator)
        , gfs(grid_operator.trialGridFunctionSpace())
        , phelper(gfs,verbose_)
        , maxiter(maxiter_)
        , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu)
        , verbose(verbose_)
        , reuse(reuse_)
        , firstapply(true)
        , usesuperlu(usesuperlu_)
      {
        params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
        params.setDebugLevel(verbose_);
#if !HAVE_SUPERLU
        if (phelper.rank() == 0 && usesuperlu == true)
          {
            std::cout << "WARNING: You are using AMG without SuperLU!"
                      << " Please consider installing SuperLU,"
                      << " or set the usesuperlu flag to false"
                      << " to suppress this warning." << std::endl;
          }
#endif
      }

      void setParameters(const Parameters& params_)
      {
        params = params_;
      }

      void setparams(Parameters params_) DUNE_DEPRECATED_MSG("setparams() is deprecated, use setParameters() instead")
      {
        params = params_;
      }

      const Parameters& parameters() const
      {
        return params;
      }

      typename V::ElementType norm (const V& v) const
      {
        V x(v); // make a copy because it has to be made consistent
        typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
        PSP psp(gfs,phelper);
        psp.make_consistent(x);
        return psp.norm(x);
      }

      void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
      {
        Timer watch;
        MatrixType& mat = Backend::native(A);
        typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
          Dune::Amg::FirstDiagonal> > Criterion;
#if HAVE_MPI
        Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
        _grid_operator.make_consistent(A);
        phelper.createIndexSetAndProjectForAMG(A, oocc);
        Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
        Operator oop(mat, oocc);
#else
        Comm oocc(gfs.gridView().comm());
        Operator oop(mat);
        Dune::SeqScalarProduct<VectorType> sp;
#endif
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;
        //use noAccu or atOnceAccu
        Criterion criterion(params);
        stats.tprepare=watch.elapsed();
        watch.reset();

        int verb=0;
        if (gfs.gridView().comm().rank()==0) verb=verbose;
        //only construct a new AMG if the matrix changes
        if (reuse==false || firstapply==true){
          amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
          firstapply = false;
          stats.tsetup = watch.elapsed();
          stats.levels = amg->maxlevels();
          stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
        }

        Dune::InverseOperatorResult stat;
        // make r consistent
        if (gfs.gridView().comm().size()>1) {
          Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
          gfs.gridView().communicate(adddh,
                                     Dune::InteriorBorder_InteriorBorder_Interface,
                                     Dune::ForwardCommunication);
        }
        watch.reset();
        Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
        solver.apply(Backend::native(z),Backend::native(r),stat);
        stats.tsolve= watch.elapsed();
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }

      const ISTLAMGStatistics& statistics() const
      {
        return stats;
      }

    private:
      const GO& _grid_operator;
      const GFS& gfs;
      PHELPER phelper;
      unsigned maxiter;
      Parameters params;
      int verbose;
      bool reuse;
      bool firstapply;
      bool usesuperlu;
      std::shared_ptr<AMG> amg;
      ISTLAMGStatistics stats;
    };


    template<class GO, int s=96>
    class ISTLBackend_NOVLP_CG_AMG_SSOR
      : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
    {

    public:
      ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
                                    int verbose_=1, bool reuse_=false,
                                    bool usesuperlu_=true)
        : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
      {}
    };

    template<class GO, int s=96>
    class ISTLBackend_NOVLP_BCGS_AMG_SSOR
      : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
    {

    public:
      ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
                                      int verbose_=1, bool reuse_=false,
                                      bool usesuperlu_=true)
        : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
      {}
    };

    template<class GO, int s=96>
    class ISTLBackend_NOVLP_LS_AMG_SSOR
      : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
    {

    public:
      ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
                                    int verbose_=1, bool reuse_=false,
                                    bool usesuperlu_=true)
        : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
      {}
    };

  } // namespace PDELab
} // namespace Dune

#endif