seqistlsolverbackend.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_PDELAB_BACKEND_ISTL_SEQISTLSOLVERBACKEND_HH
#define DUNE_PDELAB_BACKEND_ISTL_SEQISTLSOLVERBACKEND_HH
 
#include <dune/common/deprecated.hh>
#include <dune/common/parallel/mpihelper.hh>
 
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/io.hh>
#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
 
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
#include <dune/pdelab/backend/solver.hh>
#include <dune/pdelab/backend/istl/vector.hh>
#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
 
namespace Dune {
  namespace PDELab {
 
 
 
    template<typename X, typename Y, typename GO>
    class OnTheFlyOperator : public Dune::LinearOperator<X,Y>
    {
    public:
      typedef X domain_type;
      typedef Y range_type;
      typedef typename X::field_type field_type;
      static constexpr bool isLinear = GO::LocalAssembler::isLinear();
 
 
      OnTheFlyOperator (const GO& go_)
        : go(go_)
        , u_(nullptr)
      {}
 
      void setLinearizationPoint(const X& u)
      {
        u_ = &u;
      }
 
      virtual void apply (const X& x, Y& y) const override
      {
        y = 0.0;
        if (isLinear)
          go.jacobian_apply(x,y);
        else {
          if (u_ == nullptr)
            DUNE_THROW(Dune::InvalidStateException, "You seem to apply a nonlinear operator without setting the linearization point first!");
          go.jacobian_apply(*u_, x, y);
        }
      }
 
      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
      {
        Y temp(y);
        temp = 0.0;
        if (isLinear)
          go.jacobian_apply(x,temp);
        else {
          if (u_ == nullptr)
            DUNE_THROW(Dune::InvalidStateException, "You seem to apply a nonlinear operator without setting the linearization point first!");
          go.jacobian_apply(*u_, x, temp);
        }
        y.axpy(alpha,temp);
      }
 
      SolverCategory::Category category() const override
      {
        return SolverCategory::sequential;
      }
 
    private:
      const GO& go;
      const X* u_;
    };
 
    //==============================================================================
    // Here we add some standard linear solvers conforming to the linear solver
    // interface required to solve linear and nonlinear problems.
    //==============================================================================
 
    //=====================================================================
    // First we add some base classes where the preconditioner is specified
    //=====================================================================
 
    template<template<class> class Solver>
    class ISTLBackend_SEQ_Richardson
      : public SequentialNorm, public LinearResultStorage
    {
    public:
      explicit ISTLBackend_SEQ_Richardson(unsigned maxiter_=5000, int verbose_=1)
        : maxiter(maxiter_), verbose(verbose_)
      {}
 
      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)
      {
        using Backend::Native;
        using Backend::native;
 
        Dune::MatrixAdapter<Native<M>,
                            Native<V>,
                            Native<W>> opa(native(A));
        Dune::Richardson<Native<V>,Native<W> > prec(0.7);
        Solver<Native<V> > solver(opa, prec, reduction, maxiter, verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }
 
    private:
      unsigned maxiter;
      int verbose;
    };
 
    template<class GO, template<class> class Solver>
    class ISTLBackend_SEQ_MatrixFree_Richardson
      : public SequentialNorm, public LinearResultStorage
    {
      using V = typename GO::Traits::Domain;
      using W = typename GO::Traits::Range;
    public:
      explicit ISTLBackend_SEQ_MatrixFree_Richardson(const GO& go, unsigned maxiter=5000, int verbose=1)
        : opa_(go)
        , maxiter_(maxiter)
        , verbose_(verbose)
      {}
 
      void apply(V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
      {
        Dune::Richardson<V,W> prec(1.0);
        Solver<V> solver(opa_, prec, reduction, maxiter_, verbose_);
        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;
      }
 
      void setLinearizationPoint(const V& u)
      {
        opa_.setLinearizationPoint(u);
      }
 
    private:
      Dune::PDELab::OnTheFlyOperator<V,W,GO> opa_;
      unsigned maxiter_;
      int verbose_;
    };
 
    template<template<class,class,class,int> class Preconditioner,
             template<class> class Solver>
    class ISTLBackend_SEQ_Base
      : public SequentialNorm, public LinearResultStorage
    {
    public:
      explicit ISTLBackend_SEQ_Base(unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : maxiter(maxiter_), verbose(verbose_), preconditioner_steps(preconditioner_steps_)
      {}
 
 
 
      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)
      {
        using Backend::Native;
        using Backend::native;
 
        Dune::MatrixAdapter<Native<M>,
                            Native<V>,
                            Native<W>> opa(native(A));
        Preconditioner<Native<M>,
                       Native<V>,
                       Native<W>,
                       1> prec(native(A), preconditioner_steps, 1.0);
        Solver<Native<V>> solver(opa, prec, reduction, maxiter, verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }
 
    private:
      unsigned maxiter;
      int verbose;
      unsigned preconditioner_steps;
    };
 
    template<template<typename> class Solver>
    class ISTLBackend_SEQ_ILU0
      :  public SequentialNorm, public LinearResultStorage
    {
    public:
      explicit ISTLBackend_SEQ_ILU0 (unsigned maxiter_=5000, int verbose_=1)
        : maxiter(maxiter_), verbose(verbose_)
       {}
      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)
      {
        using Backend::Native;
        using Backend::native;
        Dune::MatrixAdapter<Native<M>,
                            Native<V>,
                            Native<W>> opa(native(A));
        Dune::SeqILU<Native<M>,
                      Native<V>,
                      Native<W>
                      > ilu0(native(A), 1.0);
        Solver<Native<V>> solver(opa, ilu0, reduction, maxiter, verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
       }
    private:
      unsigned maxiter;
      int verbose;
    };
 
    template<template<typename> class Solver>
    class ISTLBackend_SEQ_ILUn
      :  public SequentialNorm, public LinearResultStorage
    {
    public:
      ISTLBackend_SEQ_ILUn (int n, double w, unsigned maxiter_=5000, int verbose_=1)
        : n_(n), w_(w), maxiter(maxiter_), verbose(verbose_)
       {}
      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)
      {
        using Backend::Native;
        using Backend::native;
        Dune::MatrixAdapter<Native<M>,
                            Native<V>,
                            Native<W>
                            > opa(native(A));
        Dune::SeqILU<Native<M>,
                      Native<V>,
                      Native<W>
                      > ilun(native(A), n_, w_, false);
        Solver<Native<V>> solver(opa, ilun, reduction, maxiter, verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
       }
    private:
      int n_;
      double w_;
 
      unsigned maxiter;
      int verbose;
    };
 
    //========================================
    // Start with matrix based solver backends
    //========================================
 
 
    class ISTLBackend_SEQ_LOOP_Jac
      : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::LoopSolver>
    {
    public:
      explicit ISTLBackend_SEQ_LOOP_Jac (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::LoopSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
    class ISTLBackend_SEQ_BCGS_Richardson
      : public ISTLBackend_SEQ_Richardson<Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_BCGS_Richardson (unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_Richardson<Dune::BiCGSTABSolver>(maxiter_, verbose_)
      {}
    };
 
    class ISTLBackend_SEQ_BCGS_Jac
      : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_BCGS_Jac (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::BiCGSTABSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
    class ISTLBackend_SEQ_BCGS_SSOR
      : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_BCGS_SSOR (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::BiCGSTABSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
    class ISTLBackend_SEQ_BCGS_ILU0
      : public ISTLBackend_SEQ_ILU0<Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_BCGS_ILU0 (unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_ILU0<Dune::BiCGSTABSolver>(maxiter_, verbose_)
      {}
    };
 
    class ISTLBackend_SEQ_CG_ILU0
      : public ISTLBackend_SEQ_ILU0<Dune::CGSolver>
    {
    public:
      explicit ISTLBackend_SEQ_CG_ILU0 (unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_ILU0<Dune::CGSolver>(maxiter_, verbose_)
      {}
    };
 
    class ISTLBackend_SEQ_BCGS_ILUn
      : public ISTLBackend_SEQ_ILUn<Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_BCGS_ILUn (int n_, double w_=1.0, unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_ILUn<Dune::BiCGSTABSolver>(n_, w_, maxiter_, verbose_)
      {}
    };
 
    class ISTLBackend_SEQ_CG_ILUn
      : public ISTLBackend_SEQ_ILUn<Dune::CGSolver>
    {
    public:
      explicit ISTLBackend_SEQ_CG_ILUn (int n_, double w_=1.0, unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_ILUn<Dune::CGSolver>(n_, w_, maxiter_, verbose_)
      {}
    };
 
    class ISTLBackend_SEQ_CG_SSOR
      : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::CGSolver>
    {
    public:
      explicit ISTLBackend_SEQ_CG_SSOR (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::CGSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
    class ISTLBackend_SEQ_MINRES_SSOR
      : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::MINRESSolver>
    {
    public:
      explicit ISTLBackend_SEQ_MINRES_SSOR (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::MINRESSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
    class ISTLBackend_SEQ_CG_Jac
      : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::CGSolver>
    {
    public:
      explicit ISTLBackend_SEQ_CG_Jac (unsigned maxiter_=5000, int verbose_=1, unsigned preconditioner_steps_=1)
        : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::CGSolver>(maxiter_, verbose_, preconditioner_steps_)
      {}
    };
 
#if HAVE_SUPERLU || DOXYGEN
    class ISTLBackend_SEQ_SuperLU
      : public SequentialNorm, public LinearResultStorage
    {
    public:
      explicit ISTLBackend_SEQ_SuperLU (int verbose_=1)
        : verbose(verbose_)
      {}
 
 
      ISTLBackend_SEQ_SuperLU (int maxiter, int verbose_)
        : verbose(verbose_)
      {}
 
      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)
      {
        using Backend::Native;
        using Backend::native;
        using ISTLM = Native<M>;
        Dune::SuperLU<ISTLM> solver(native(A), verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }
 
    private:
      int verbose;
    };
#endif // HAVE_SUPERLU || DOXYGEN
 
#if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
    class ISTLBackend_SEQ_UMFPack
      : public SequentialNorm, public LinearResultStorage
    {
    public:
      explicit ISTLBackend_SEQ_UMFPack (int verbose_=1)
        : verbose(verbose_)
      {}
 
 
      ISTLBackend_SEQ_UMFPack (int maxiter, int verbose_)
        : verbose(verbose_)
      {}
 
      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)
      {
        using Backend::native;
        using ISTLM = Backend::Native<M>;
        Dune::UMFPack<ISTLM> solver(native(A), verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }
 
    private:
      int verbose;
    };
#endif // HAVE_SUITESPARSE_UMFPACK || DOXYGEN
 
    class ISTLBackend_SEQ_ExplicitDiagonal
      : public SequentialNorm, public LinearResultStorage
    {
    public:
      ISTLBackend_SEQ_ExplicitDiagonal ()
      {}
 
      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)
      {
        using Backend::Native;
        Dune::SeqJac<Native<M>,
                     Native<V>,
                     Native<W>
                     > jac(Backend::native(A),1,1.0);
        jac.pre(z,r);
        jac.apply(z,r);
        jac.post(z);
        res.converged  = true;
        res.iterations = 1;
        res.elapsed    = 0.0;
        res.reduction  = reduction;
        res.conv_rate  = reduction; // pow(reduction,1.0/1)
      }
    };
 
 
    struct ISTLAMGStatistics
    {
      double tprepare;
      int levels;
      double tsolve;
      double tsetup;
      int iterations;
      bool directCoarseLevelSolver;
    };
 
    template<class GO, template<class,class,class,int> class Preconditioner, template<class> class Solver,
              bool skipBlocksizeCheck = false>
    class ISTLBackend_SEQ_AMG : public LinearResultStorage
    {
      typedef typename GO::Traits::TrialGridFunctionSpace GFS;
      typedef typename GO::Traits::Jacobian M;
      typedef Backend::Native<M> MatrixType;
      typedef typename GO::Traits::Domain V;
      typedef Backend::Native<V> VectorType;
      typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
      typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
      typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
      typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG;
      typedef Dune::Amg::Parameters Parameters;
 
    public:
      ISTLBackend_SEQ_AMG(unsigned maxiter_=5000, int verbose_=1,
                          bool reuse_=false, bool usesuperlu_=true)
        : maxiter(maxiter_), params(15,2000), verbose(verbose_),
          reuse(reuse_), firstapply(true), usesuperlu(usesuperlu_)
      {
        params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
        params.setDebugLevel(verbose_);
#if !HAVE_SUPERLU
        if (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 setparams(Parameters params_)
      {
        params = params_;
      }
 
      void setReuse(bool reuse_)
      {
        reuse = reuse_;
      }
 
      bool getReuse() const
      {
        return reuse;
      }
 
      typename V::ElementType norm (const V& v) const
      {
        return Backend::native(v).two_norm();
      }
 
      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;
        SmootherArgs smootherArgs;
        smootherArgs.iterations = 1;
        smootherArgs.relaxationFactor = 1;
 
        Criterion criterion(params);
        //only construct a new AMG if the matrix changes
        if (reuse==false || firstapply==true){
          oop.reset(new Operator(mat));
          amg.reset(new AMG(*oop, criterion, smootherArgs));
          firstapply = false;
          stats.tsetup = watch.elapsed();
          stats.levels = amg->maxlevels();
          stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
        }
        watch.reset();
        Dune::InverseOperatorResult stat;
 
        Solver<VectorType> solver(*oop,*amg,reduction,maxiter,verbose);
        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:
      unsigned maxiter;
      Parameters params;
      int verbose;
      bool reuse;
      bool firstapply;
      bool usesuperlu;
      std::shared_ptr<Operator> oop;
      std::shared_ptr<AMG> amg;
      ISTLAMGStatistics stats;
    };
 
 
    template<class GO>
    class ISTLBackend_SEQ_CG_AMG_SSOR
      : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::CGSolver>
    {
 
    public:
      ISTLBackend_SEQ_CG_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
                                  bool reuse_=false, bool usesuperlu_=true)
        : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::CGSolver>
          (maxiter_, verbose_, reuse_, usesuperlu_)
      {}
    };
 
    template<class GO>
    class ISTLBackend_SEQ_BCGS_AMG_SSOR
      : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::BiCGSTABSolver>
    {
 
    public:
      ISTLBackend_SEQ_BCGS_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
                                    bool reuse_=false, bool usesuperlu_=true)
        : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::BiCGSTABSolver>
          (maxiter_, verbose_, reuse_, usesuperlu_)
      {}
    };
 
    template<class GO>
    class ISTLBackend_SEQ_BCGS_AMG_SOR
      : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::BiCGSTABSolver>
    {
 
    public:
      ISTLBackend_SEQ_BCGS_AMG_SOR(unsigned maxiter_=5000, int verbose_=1,
                                   bool reuse_=false, bool usesuperlu_=true)
        : ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::BiCGSTABSolver>
          (maxiter_, verbose_, reuse_, usesuperlu_)
      {}
    };
 
    template<class GO>
    class ISTLBackend_SEQ_LS_AMG_SSOR
      : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::LoopSolver>
    {
 
    public:
      ISTLBackend_SEQ_LS_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
                                  bool reuse_=false, bool usesuperlu_=true)
        : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::LoopSolver>
          (maxiter_, verbose_, reuse_, usesuperlu_)
      {}
    };
 
    template<class GO>
    class ISTLBackend_SEQ_LS_AMG_SOR
      : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::LoopSolver>
    {
 
    public:
      ISTLBackend_SEQ_LS_AMG_SOR(unsigned maxiter_=5000, int verbose_=1,
                                 bool reuse_=false, bool usesuperlu_=true)
        : ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::LoopSolver>
          (maxiter_, verbose_, reuse_, usesuperlu_)
      {}
    };
 
    class ISTLBackend_SEQ_GMRES_ILU0
      : public SequentialNorm, public LinearResultStorage
    {
    public :
 
      explicit ISTLBackend_SEQ_GMRES_ILU0(int restart_ = 200, int maxiter_ = 5000, int verbose_ = 1)
        : restart(restart_), maxiter(maxiter_), verbose(verbose_)
      {}
 
      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)
      {
        using Backend::Native;
        using Backend::native;
        Dune::MatrixAdapter<
          Native<M>,
          Native<V>,
          Native<W>
          > opa(native(A));
        Dune::SeqILU<
          Native<M>,
          Native<V>,
          Native<W>
          > ilu0(native(A), 1.0);
        Dune::RestartedGMResSolver<Native<V>> solver(opa,ilu0,reduction,restart,maxiter,verbose);
        Dune::InverseOperatorResult stat;
        solver.apply(native(z), native(r), stat);
        res.converged  = stat.converged;
        res.iterations = stat.iterations;
        res.elapsed    = stat.elapsed;
        res.reduction  = stat.reduction;
        res.conv_rate  = stat.conv_rate;
      }
 
    private :
      int restart, maxiter, verbose;
    };
 
    //============================
    // Matrix free solver backends
    //============================
 
    template <class GO>
    class ISTLBackend_SEQ_MatrixFree_BCGS_Richardson
      : public ISTLBackend_SEQ_MatrixFree_Richardson<GO, Dune::BiCGSTABSolver>
    {
    public:
      explicit ISTLBackend_SEQ_MatrixFree_BCGS_Richardson (const GO& go_, unsigned maxiter_=5000, int verbose_=1)
        : ISTLBackend_SEQ_MatrixFree_Richardson<GO, Dune::BiCGSTABSolver>(go_, maxiter_, verbose_)
      {}
    };
 
 
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_BACKEND_ISTL_SEQISTLSOLVERBACKEND_HH