Compare revisions

Christian Engwer · Christian Engwer · Markus Blatt · Nils-Arne Dreier · Nils-Arne Dreier · Nils-Arne Dreier
--- a/dune/istl/cholmod.hh
+++ b/dune/istl/cholmod.hh
@@ -11,7 +11,7 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include<dune/istl/solver.hh>
-#include <dune/istl/solverfactory.hh>
+#include <dune/istl/solverregistry.hh>
 #include <dune/istl/foreach.hh>

 #include <vector>
@@ -250,7 +250,7 @@ public:

  /** @brief simple forward to apply(X&, Y&, InverseOperatorResult&)
    */
-  void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+  void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
  {
    apply(x,b,res);
  }
@@ -260,7 +260,7 @@ public:
   * The method assumes that setMatrix() was called before
   * In the case of a given ignore field the corresponding entries of both in x and b will stay untouched in this method.
  */
-  void apply(Vector& x, Vector& b, InverseOperatorResult& res)
+  void apply(Vector& x, Vector& b, InverseOperatorResult& res) override
  {
    // do nothing if N=0
    if ( nIsZero_ )
@@ -468,7 +468,7 @@ public:
    CholmodMethod::factorize(M.get(), L_, &c_);
  }

-  virtual SolverCategory::Category category() const
+  SolverCategory::Category category() const override
  {
    return SolverCategory::Category::sequential;
  }
@@ -524,34 +524,31 @@ private:
  std::vector<std::size_t> subIndices_;
 };

-  struct CholmodCreator{
-    template<class F> struct isValidBlock : std::false_type{};
-    template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
-    template<int k> struct isValidBlock<FieldVector<float,k>> : std::true_type{};
-
-    template<class TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator()(TL /*tl*/, const M& mat, const Dune::ParameterTree& /*config*/,
-               std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      using D = typename Dune::TypeListElement<1, TL>::type;
-      auto solver = std::make_shared<Dune::Cholmod<D>>();
-      solver->setMatrix(mat);
-      return solver;
-    }
-
-    // second version with SFINAE to validate the template parameters of Cholmod
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
-                std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      DUNE_THROW(UnsupportedType, "Unsupported Type in Cholmod");
-    }
-  };
-  DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator());
+  DUNE_REGISTER_SOLVER("cholmod",
+                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
+                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
+                       {
+                         using OpTraits = decltype(opTraits);
+                         using M = typename OpTraits::matrix_type;
+                         using D = typename OpTraits::domain_type;
+                         // works only for sequential operators
+                         if constexpr (OpTraits::isParallel){
+                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
+                             DUNE_THROW(Dune::InvalidStateException, "CholMod works only for sequential operators.");
+                         }
+                         if constexpr (OpTraits::isAssembled &&
+                                       // check whether the Matrix field_type is double or float
+                                       (std::is_same_v<typename FieldTraits<D>::field_type, double> ||
+                                       std::is_same_v<typename FieldTraits<D>::field_type, float>)){
+                           const auto& A = opTraits.getAssembledOpOrThrow(op);
+                           const M& mat = A->getmat();
+                           auto solver = std::make_shared<Dune::Cholmod<D>>();
+                           solver->setMatrix(mat);
+                           return solver;
+                         }
+                         DUNE_THROW(UnsupportedType,
+                                    "Unsupported Type in Cholmod (only double and float supported)");
+                       });

 } /* namespace Dune */


--- a/dune/istl/common/registry.hh
+++ b/dune/istl/common/registry.hh
@@ -53,7 +53,7 @@ namespace Dune {
    }

    /*
-      Register all creators from the registry in the Parameterizedobjectfactory An
+      Register all creators from the registry in the Parameterizedobjectfactory. An
      object of V is passed in the creator and should be used to determine the
      template arguments.
    */

--- a/dune/istl/eigenvalue/poweriteration.hh
+++ b/dune/istl/eigenvalue/poweriteration.hh
@@ -56,13 +56,13 @@ namespace Dune
          mutable_scaling_(mutable_scaling)
      {}

-      virtual void apply (const X& x, Y& y) const
+      void apply (const X& x, Y& y) const override
      {
        y = x;
        y *= immutable_scaling_*mutable_scaling_;
      }

-      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+      void applyscaleadd (field_type alpha, const X& x, Y& y) const override
      {
        X temp(x);
        temp *= immutable_scaling_*mutable_scaling_;
@@ -70,7 +70,7 @@ namespace Dune
      }

      //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
      {
        return SolverCategory::sequential;
      }
@@ -108,14 +108,14 @@ namespace Dune
          "Range type of both operators doesn't match!");
      }

-      virtual void apply (const domain_type& x, range_type& y) const
+      void apply (const domain_type& x, range_type& y) const override
      {
        op1_.apply(x,y);
        op2_.applyscaleadd(1.0,x,y);
      }

-      virtual void applyscaleadd (field_type alpha,
-        const domain_type& x, range_type& y) const
+      void applyscaleadd (field_type alpha,
+        const domain_type& x, range_type& y) const override
      {
        range_type temp(y);
        op1_.apply(x,temp);
@@ -124,7 +124,7 @@ namespace Dune
      }

      //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
      {
        return SolverCategory::sequential;
      }

--- a/dune/istl/ldl.hh
+++ b/dune/istl/ldl.hh
@@ -24,7 +24,7 @@ extern "C"
 #include <dune/istl/bccsmatrixinitializer.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/solvertype.hh>
-#include <dune/istl/solverfactory.hh>
+#include <dune/istl/solverregistry.hh>

 namespace Dune {
  /**
@@ -87,7 +87,7 @@ namespace Dune {
    typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;

    //! Category of the solver (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::Category::sequential;
    }
@@ -149,7 +149,7 @@ namespace Dune {
    }

    /** \copydoc InverseOperator::apply(X&, Y&, InverseOperatorResult&) */
-    virtual void apply(domain_type& x, range_type& b, InverseOperatorResult& res)
+    void apply(domain_type& x, range_type& b, InverseOperatorResult& res) override
    {
      const int dimMat(ldlMatrix_.N());
      ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
@@ -163,7 +163,7 @@ namespace Dune {
    }

    /** \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) */
-    virtual void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+    void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
    {
      apply(x,b,res);
    }
@@ -370,34 +370,33 @@ namespace Dune {
    enum {value = true};
  };

-  struct LDLCreator {
-    template<class F> struct isValidBlock : std::false_type{};
-    template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
-
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
-      std::enable_if_t<
-                isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      int verbose = config.get("verbose", 0);
-      return std::make_shared<Dune::LDL<M>>(mat,verbose);
-    }
-
-    // second version with SFINAE to validate the template parameters of LDL
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
-      std::enable_if_t<
-                !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      DUNE_THROW(UnsupportedType,
-        "Unsupported Type in LDL (only double and std::complex<double> supported)");
-    }
-  };
-  DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator());
+  DUNE_REGISTER_SOLVER("ldl",
+                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
+                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
+                       {
+                         using OpTraits = decltype(opTraits);
+                         using M = typename OpTraits::matrix_type;
+                         // works only for sequential operators
+                         if constexpr (OpTraits::isParallel){
+                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
+                             DUNE_THROW(Dune::InvalidStateException, "LDL works only for sequential operators.");
+                         }
+                         // check if LDL<M>* is convertible to
+                         // InverseOperator*. This allows only the explicit
+                         // specialized variants of LDL
+                         if constexpr (std::is_convertible_v<LDL<M>*,
+                                       Dune::InverseOperator<typename OpTraits::domain_type,
+                                       typename OpTraits::range_type>*> &&
+                                       std::is_same_v<typename FieldTraits<M>::field_type, double>
+                                       ){
+                           const auto& A = opTraits.getAssembledOpOrThrow(op);
+                           const M& mat = A->getmat();
+                           int verbose = config.get("verbose", 0);
+                           return std::make_shared<LDL<M>>(mat,verbose);
+                         }
+                         DUNE_THROW(UnsupportedType,
+                                    "Unsupported Type in LDL (only FieldMatrix<double,...> supported)");
+                       });

 } // end namespace Dune


--- a/dune/istl/novlpschwarz.hh
+++ b/dune/istl/novlpschwarz.hh
@@ -98,7 +98,7 @@ namespace Dune {
    {}

    //! apply operator to x:  \f$ y = A(x) \f$
-    virtual void apply (const X& x, Y& y) const
+    void apply (const X& x, Y& y) const override
    {
      y = 0;
      novlp_op_apply(x,y,1);
@@ -106,7 +106,7 @@ namespace Dune {
    }

    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
-    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+    void applyscaleadd (field_type alpha, const X& x, Y& y) const override
    {
      // only apply communication to alpha*A*x to make it consistent,
      // y already has to be consistent.
@@ -117,8 +117,8 @@ namespace Dune {
      y += y1;
    }

-    //! get matrix via *
-    virtual const matrix_type& getmat () const
+    //! get reference to matrix
+    const M& getmat () const override
    {
      return *_A_;
    }
@@ -232,7 +232,7 @@ namespace Dune {
    }

    //! Category of the linear operator (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::nonoverlapping;
    }
@@ -312,7 +312,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(domain_type&,range_type&)
     */
-    virtual void pre (domain_type& x, range_type& b)
+    void pre (domain_type& x, range_type& b) override
    {
      _preconditioner->pre(x,b);
    }
@@ -322,7 +322,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(domain_type&,const range_type&)
     */
-    virtual void apply (domain_type& v, const range_type& d)
+    void apply (domain_type& v, const range_type& d) override
    {
      // block preconditioner equivalent to WrappedPreconditioner from
      // pdelab/backend/ovlpistsolverbackend.hh,
@@ -343,13 +343,13 @@ namespace Dune {

       \copydoc Preconditioner::post(domain_type&)
     */
-    virtual void post (domain_type& x)
+    void post (domain_type& x) override
    {
      _preconditioner->post(x);
    }

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::nonoverlapping;
    }

--- a/dune/istl/operators.hh
+++ b/dune/istl/operators.hh
 // SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 // vi: set et ts=4 sw=2 sts=2:
 #ifndef DUNE_ISTL_OPERATORS_HH
 #define DUNE_ISTL_OPERATORS_HH
@@ -116,7 +115,7 @@ namespace Dune {
    typedef Y range_type;
    typedef typename X::field_type field_type;

-    //! get matrix via *
+    //! get reference to matrix
    virtual const M& getmat () const = 0;
  };

@@ -159,7 +158,7 @@ namespace Dune {
      _A_->usmv(alpha,x,y);
    }

-    //! get matrix via *
+    //! get reference to matrix
    const M& getmat () const override
    {
      return *_A_;
@@ -180,3 +179,10 @@ namespace Dune {
 } // end namespace

 #endif
+
+// Emacs configuration
+// Local Variables:
+// tab-width: 4
+// indent-tabs-mode: nil
+// c-basic-offset: 2
+// End:
--- a/dune/istl/overlappingschwarz.hh
+++ b/dune/istl/overlappingschwarz.hh
@@ -843,7 +843,7 @@ namespace Dune

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] X& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] X& b) override
    {}

    /*!
@@ -851,21 +851,21 @@ namespace Dune

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const X& d);
+    void apply (X& v, const X& d) override;

    /*!
       \brief Postprocess the preconditioner.

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    template<bool forward>
    void apply(X& v, const X& d);

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -16,6 +16,7 @@
 #include <dune/istl/superlu.hh>
 #include <dune/istl/umfpack.hh>
 #include <dune/istl/solvertype.hh>
+#include <dune/istl/solverregistry.hh>
 #include <dune/common/typetraits.hh>
 #include <dune/common/exceptions.hh>
 #include <dune/common/scalarvectorview.hh>
@@ -1167,8 +1168,8 @@ namespace Dune
  } // end namespace Amg

  struct AMGCreator{
-    template<class> struct isValidBlockType : std::false_type{};
-    template<class T, int n, int m> struct isValidBlockType<FieldMatrix<T,n,m>> : std::true_type{};
+    template<class> struct isValidMatrix : std::false_type{};
+    template<class T, int n, int m, class A> struct isValidMatrix<BCRSMatrix<FieldMatrix<T,n,m>, A>> : std::true_type{};

    template<class OP>
    std::shared_ptr<Dune::Preconditioner<typename OP::element_type::domain_type, typename OP::element_type::range_type> >
@@ -1242,13 +1243,13 @@ namespace Dune
        DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
    }

-    template<typename TL, typename OP>
-    std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
-                                         typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL tl, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
-                std::enable_if_t<isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
+    template<typename OpTraits, typename OP>
+    std::shared_ptr<Dune::Preconditioner<typename OpTraits::domain_type,
+                                         typename OpTraits::range_type>>
+    operator() (OpTraits opTraits, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
+                std::enable_if_t<isValidMatrix<typename OpTraits::matrix_type>::value,int> = 0) const
    {
-      using field_type = typename OP::matrix_type::field_type;
+      using field_type = typename OpTraits::matrix_type::field_type;
      using real_type = typename FieldTraits<field_type>::real_type;
      if (!std::is_convertible<field_type, real_type>())
        DUNE_THROW(UnsupportedType, "AMG needs field_type(" <<
@@ -1257,14 +1258,19 @@ namespace Dune
                   className<real_type>() <<
                   ").");
      std::string smoother = config.get("smoother", "ssor");
+      // we can irgnore the OpTraits here. As the AMG can only work
+      // with actual matrices, the operator op must be of type
+      // MatrixAdapter or *SchwarzOperator. In any case these
+      // operators provide all necessary information about matrix,
+      // domain and range type
      return makeAMG(op, smoother, config);
    }

-    template<typename TL, typename OP>
-    std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
-                                         typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const std::shared_ptr<OP>& /*mat*/, const Dune::ParameterTree& /*config*/,
-                std::enable_if_t<!isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
+    template<typename OpTraits, typename OP>
+    std::shared_ptr<Dune::Preconditioner<typename OpTraits::domain_type,
+                                         typename OpTraits::range_type>>
+    operator() (OpTraits opTraits, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
+                std::enable_if_t<!isValidMatrix<typename OpTraits::matrix_type>::value,int> = 0) const
    {
      DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
    }

--- a/dune/istl/paamg/twolevelmethod.hh
+++ b/dune/istl/paamg/twolevelmethod.hh
@@ -273,7 +273,7 @@ private:
      : amg_(op, crit,args), first_(true)
    {}

-    void apply(X& x, X& b, [[maybe_unused]] double reduction, [[maybe_unused]] InverseOperatorResult& res)
+    void apply(X& x, X& b, [[maybe_unused]] double reduction, [[maybe_unused]] InverseOperatorResult& res) override
    {
      if(first_)
      {
@@ -284,13 +284,13 @@ private:
      amg_.apply(x,b);
    }

-    void apply(X& x, X& b, InverseOperatorResult& res)
+    void apply(X& x, X& b, InverseOperatorResult& res) override
    {
      return apply(x,b,1e-8,res);
    }

    //! Category of the solver (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return amg_.category();
    }
@@ -428,15 +428,15 @@ public:
    delete coarseSolver_;
  }

-  void pre(FineDomainType& x, FineRangeType& b)
+  void pre(FineDomainType& x, FineRangeType& b) override
  {
    smoother_->pre(x,b);
  }

-  void post([[maybe_unused]] FineDomainType& x)
+  void post([[maybe_unused]] FineDomainType& x) override
  {}

-  void apply(FineDomainType& v, const FineRangeType& d)
+  void apply(FineDomainType& v, const FineRangeType& d) override
  {
    FineDomainType u(v);
    FineRangeType rhs(d);
@@ -463,7 +463,7 @@ public:
  }

  //! Category of the preconditioner (see SolverCategory::Category)
-  virtual SolverCategory::Category category() const
+  SolverCategory::Category category() const override
  {
    return SolverCategory::sequential;
  }

--- a/dune/istl/preconditioners.hh
+++ b/dune/istl/preconditioners.hh
@@ -99,21 +99,21 @@ namespace Dune {
        DUNE_THROW(InvalidStateException, "User-supplied solver category does not match that of the given inverse operator");
    }

-    virtual void pre(domain_type&,range_type&)
+    void pre([[maybe_unused]] domain_type&,[[maybe_unused]] range_type&) override
    {}

-    virtual void apply(domain_type& v, const range_type& d)
+    void apply(domain_type& v, const range_type& d) override
    {
      InverseOperatorResult res;
      range_type copy(d);
      inverse_operator_.apply(v, copy, res);
    }

-    virtual void post(domain_type&)
+    void post([[maybe_unused]] domain_type&) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::category(inverse_operator_);
    }
@@ -206,7 +206,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    /*!
@@ -214,7 +214,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      for (int i=0; i<_n; i++) {
        bsorf(_A_,v,d,_w,BL<l>());
@@ -227,11 +227,11 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -326,7 +326,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    /*!
@@ -334,7 +334,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      this->template apply<true>(v,d);
    }
@@ -365,11 +365,11 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -477,7 +477,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    /*!
@@ -485,7 +485,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      for (int i=0; i<_n; i++) {
        dbjac(_A_,v,d,_w,BL<l>());
@@ -497,11 +497,11 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -636,7 +636,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre([[maybe_unused]] X &x, [[maybe_unused]] Y &b)
+    void pre([[maybe_unused]] X &x, [[maybe_unused]] Y &b) override
    {
    }

@@ -645,7 +645,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply(X &v, const Y &d)
+    void apply(X &v, const Y &d) override
    {

      DILU::blockDILUBacksolve(_A_, Dinv_, v, d);
@@ -661,12 +661,12 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post([[maybe_unused]] X &x)
+    void post([[maybe_unused]] X &x) override
    {
    }

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -810,7 +810,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    /*!
@@ -818,7 +818,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      if( ILU_ )
      {
@@ -840,11 +840,11 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -917,7 +917,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
+    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    /*!
@@ -925,7 +925,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      v = d;
      v *= _w;
@@ -936,11 +936,11 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x)
+    void post ([[maybe_unused]] X& x) override
    {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }
@@ -949,11 +949,12 @@ namespace Dune {
    //! \brief The relaxation factor to use.
    real_field_type _w;
  };
-  DUNE_REGISTER_PRECONDITIONER("richardson", [](auto tl, const auto& /* mat */, const ParameterTree& config){
-                                               using D = typename Dune::TypeListElement<1, decltype(tl)>::type;
-                                               using R = typename Dune::TypeListElement<2, decltype(tl)>::type;
-                                               return std::make_shared<Richardson<D,R>>(config);
-                                             });
+  DUNE_REGISTER_PRECONDITIONER("richardson", [](auto opTraits, const auto& op, const ParameterTree& config){
+    using OpTraits = std::decay_t<decltype(opTraits)>;
+    using D = typename OpTraits::domain_type;
+    using R = typename OpTraits::range_type;
+    return std::make_shared<Richardson<D,R>>(config);
+  });


  /**

--- a/dune/istl/scalarproducts.hh
+++ b/dune/istl/scalarproducts.hh
@@ -136,7 +136,7 @@ namespace Dune {
       It is assumed that the vectors are consistent on the interior+border
       partition.
     */
-    virtual field_type dot (const X& x, const X& y) const override
+    field_type dot (const X& x, const X& y) const override
    {
      field_type result(0);
      _communication->dot(x,y,result); // explicitly loop and apply masking
@@ -146,13 +146,13 @@ namespace Dune {
    /*! \brief Norm of a right-hand side vector.
       The vector must be consistent on the interior+border partition
     */
-    virtual real_type norm (const X& x) const override
+    real_type norm (const X& x) const override
    {
      return _communication->norm(x);
    }

    //! Category of the scalar product (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const override
+    SolverCategory::Category category() const override
    {
      return _category;
    }

--- a/dune/istl/schwarz.hh
+++ b/dune/istl/schwarz.hh
@@ -113,7 +113,7 @@ namespace Dune {
    {}

    //! apply operator to x:  \f$ y = A(x) \f$
-    virtual void apply (const X& x, Y& y) const
+    void apply (const X& x, Y& y) const override
    {
      y = 0;
      _A_->umv(x,y);     // result is consistent on interior+border
@@ -122,7 +122,7 @@ namespace Dune {
    }

    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
-    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+    void applyscaleadd (field_type alpha, const X& x, Y& y) const override
    {
      _A_->usmv(alpha,x,y);     // result is consistent on interior+border
      communication.project(y);     // we want this here to avoid it before the preconditioner
@@ -130,13 +130,13 @@ namespace Dune {
    }

    //! get the sequential assembled linear operator.
-    virtual const matrix_type& getmat () const
+    const matrix_type& getmat () const override
    {
      return *_A_;
    }

    //! Category of the linear operator (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::overlapping;
    }
@@ -203,7 +203,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre (X& x, [[maybe_unused]] Y& b)
+    void pre (X& x, [[maybe_unused]] Y& b) override
    {
      communication.copyOwnerToAll(x,x);     // make dirichlet values consistent
    }
@@ -213,7 +213,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      for (int i=0; i<_n; i++) {
        bsorf(_A_,v,d,_w);
@@ -227,10 +227,10 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post ([[maybe_unused]] X& x) {}
+    void post ([[maybe_unused]] X& x) override {}

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::overlapping;
    }
@@ -323,7 +323,7 @@ namespace Dune {

       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre (X& x, Y& b)
+    void pre (X& x, Y& b) override
    {
      _communication.copyOwnerToAll(x,x);     // make dirichlet values consistent
      _preconditioner->pre(x,b);
@@ -334,7 +334,7 @@ namespace Dune {

       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (X& v, const Y& d)
+    void apply (X& v, const Y& d) override
    {
      _preconditioner->apply(v,d);
      _communication.copyOwnerToAll(v,v);
@@ -352,13 +352,13 @@ namespace Dune {

       \copydoc Preconditioner::post(X&)
     */
-    virtual void post (X& x)
+    void post (X& x) override
    {
      _preconditioner->post(x);
    }

    //! Category of the preconditioner (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::overlapping;
    }

--- a/dune/istl/solver.hh
+++ b/dune/istl/solver.hh
@@ -373,7 +373,7 @@ namespace Dune

       \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
     */
-    virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
+    void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
    {
      scalar_real_type saved_reduction = _reduction;
      _reduction = reduction;
@@ -382,7 +382,7 @@ namespace Dune
    }

    //! Category of the solver (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return _category;
    }

--- a/dune/istl/solverfactory.hh
+++ b/dune/istl/solverfactory.hh
@@ -11,9 +11,12 @@
 #include <memory>

 #include <dune/common/parametertree.hh>
+#include <dune/common/std/type_traits.hh>
 #include <dune/common/singleton.hh>
+#include <dune/common/parameterizedobject.hh>

-#include "solverregistry.hh"
+#include <dune/istl/solverregistry.hh>
+#include <dune/istl/common/registry.hh>
 #include <dune/istl/solver.hh>
 #include <dune/istl/schwarz.hh>
 #include <dune/istl/novlpschwarz.hh>
@@ -23,92 +26,107 @@ namespace Dune{
      @{
  */

-  // Direct solver factory:
-  template<class M, class X, class Y>
-  using DirectSolverSignature = std::shared_ptr<InverseOperator<X,Y>>(const M&, const ParameterTree&);
-  template<class M, class X, class Y>
-  using DirectSolverFactory = Singleton<ParameterizedObjectFactory<DirectSolverSignature<M,X,Y>>>;
-
  // Preconditioner factory:
-  template<class M, class X, class Y>
-  using PreconditionerSignature = std::shared_ptr<Preconditioner<X,Y>>(const std::shared_ptr<M>&, const ParameterTree&);
-  template<class M, class X, class Y>
-  using PreconditionerFactory = Singleton<ParameterizedObjectFactory<PreconditionerSignature<M,X,Y>>>;
+  template<class OP>
+  using PreconditionerSignature = std::shared_ptr<Preconditioner<typename OP::domain_type,typename OP::range_type>>(const std::shared_ptr<OP>&, const ParameterTree&);
+  template<class OP>
+  using PreconditionerFactory = Singleton<ParameterizedObjectFactory<PreconditionerSignature<OP>>>;

  // Iterative solver factory
-  template<class X, class Y>
-  using IterativeSolverSignature = std::shared_ptr<InverseOperator<X,Y>>(const std::shared_ptr<LinearOperator<X,Y>>&, const std::shared_ptr<ScalarProduct<X>>&, const std::shared_ptr<Preconditioner<X,Y>>, const ParameterTree&);
-  template<class X, class Y>
-  using IterativeSolverFactory = Singleton<ParameterizedObjectFactory<IterativeSolverSignature<X,Y>>>;
+  template<class OP>
+  using SolverSignature = std::shared_ptr<InverseOperator<typename OP::domain_type,typename OP::range_type>>(const std::shared_ptr<OP>&, const ParameterTree&);
+  template<class OP>
+  using SolverFactory = Singleton<ParameterizedObjectFactory<SolverSignature<OP>>>;
+
+  template<class Operator>
+  struct OperatorTraits{
+  private:
+    template<class O>
+    using _matrix_type = typename O::matrix_type;
+    template<class O>
+    using _comm_type = typename O::communication_type;
+  public:
+    using domain_type = typename Operator::domain_type;
+    using range_type = typename Operator::range_type;
+    using operator_type = Operator;
+    using solver_type = InverseOperator<domain_type, range_type>;
+    using matrix_type = Std::detected_or_t<int, _matrix_type, Operator>;
+    static constexpr bool isAssembled = !std::is_same<matrix_type, int>::value;
+    using comm_type = Std::detected_or_t<int, _comm_type, Operator>;
+    static constexpr bool isParallel = !std::is_same<comm_type, int>::value;
+
+    static const std::shared_ptr<AssembledLinearOperator<matrix_type, domain_type, range_type>>
+    getAssembledOpOrThrow(std::shared_ptr<LinearOperator<domain_type, range_type>> op){
+      std::shared_ptr<AssembledLinearOperator<matrix_type, domain_type, range_type>> aop
+        = std::dynamic_pointer_cast<AssembledLinearOperator<matrix_type, domain_type, range_type>>(op);
+      if(aop)
+          return aop;
+      DUNE_THROW(NoAssembledOperator, "Failed to cast to AssembledLinearOperator. Please pass in an AssembledLinearOperator.");
+    }
+
+    static const comm_type& getCommOrThrow(std::shared_ptr<LinearOperator<domain_type, range_type>> op){
+      std::shared_ptr<Operator> _op
+        = std::dynamic_pointer_cast<Operator>(op);
+      if constexpr (isParallel){
+        return _op->getCommunication();
+      }else{
+        DUNE_THROW(NoAssembledOperator, "Could not obtain communication object from operator. Please pass in a parallel operator.");
+      }
+    }
+
+    static std::shared_ptr<ScalarProduct<domain_type>> getScalarProduct(std::shared_ptr<LinearOperator<domain_type, range_type>> op)
+    {
+      if constexpr (isParallel){
+        return createScalarProduct<domain_type>(getCommOrThrow(op), op->category());
+      }else{
+        return std::make_shared<SeqScalarProduct<domain_type>>();
+      }
+    }
+  };
+
+  template<class Operator>
+  struct [[deprecated("Please change to the new solverfactory interface.")]]
+  TypeListElement<0, OperatorTraits<Operator>>{
+    using type [[deprecated]] = typename OperatorTraits<Operator>::matrix_type;
+    using Type [[deprecated]] = type;
+  };
+  template<class Operator>
+  struct [[deprecated("Please change to the new solverfactory interface.")]]
+  TypeListElement<1, OperatorTraits<Operator>>{
+    using type [[deprecated]] = typename OperatorTraits<Operator>::domain_type;
+    using Type [[deprecated]] = type;
+  };
+  template<class Operator>
+  struct [[deprecated("Please change to the new solverfactory interface.")]]
+  TypeListElement<2, OperatorTraits<Operator>>{
+    using type [[deprecated]] = typename OperatorTraits<Operator>::range_type;
+    using Type [[deprecated]] = type;
+  };

  // initSolverFactories differs in different compilation units, so we have it
  // in an anonymous namespace
  namespace {

-    /** initializes the direct solvers, preconditioners and iterative solvers in
-        the factories with the corresponding Matrix and Vector types.
+    /** initializes the preconditioners and solvers in
+        the factories with the corresponding Operator and Vector types.

       @tparam O the assembled linear operator type
    */
    template<class O>
    int initSolverFactories(){
-      using M  = typename O::matrix_type;
-      using X  = typename O::range_type;
-      using Y  = typename O::domain_type;
-      using TL = Dune::TypeList<M,X,Y>;
-      auto& dsfac=Dune::DirectSolverFactory<M,X,Y>::instance();
-      addRegistryToFactory<TL>(dsfac, DirectSolverTag{});
-      auto& pfac=Dune::PreconditionerFactory<O,X,Y>::instance();
-      addRegistryToFactory<TL>(pfac, PreconditionerTag{});
-      using TLS = Dune::TypeList<X,Y>;
-      auto& isfac=Dune::IterativeSolverFactory<X,Y>::instance();
-      return addRegistryToFactory<TLS>(isfac, IterativeSolverTag{});
+      using OpInfo = OperatorTraits<O>;
+      auto& pfac=Dune::PreconditionerFactory<O>::instance();
+      addRegistryToFactory<OpInfo>(pfac, PreconditionerTag{});
+      auto& isfac=Dune::SolverFactory<O>::instance();
+      return addRegistryToFactory<OpInfo>(isfac, SolverTag{});
    }
  } // end anonymous namespace

-
-  template<class O, class Preconditioner>
-  std::shared_ptr<Preconditioner> wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec,
-                                                              const O&)
-  {
-    return prec;
-  }
-
-  template<class M, class X, class Y, class C, class Preconditioner>
-  std::shared_ptr<Preconditioner>
-  wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec,
-                              const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C> >& op)
-  {
-    return std::make_shared<BlockPreconditioner<X,Y,C,Preconditioner> >(prec, op->getCommunication());
-  }
-
-  template<class M, class X, class Y, class C, class Preconditioner>
-  std::shared_ptr<Preconditioner>
-  wrapPreconditioner4Parallel(const std::shared_ptr<Preconditioner>& prec,
-                              const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C> >& op)
-  {
-    return std::make_shared<NonoverlappingBlockPreconditioner<C,Preconditioner> >(prec, op->getCommunication());
-  }
-
-  template<class M, class X, class Y>
-  std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<MatrixAdapter<M,X,Y> >&)
-  {
-    return std::make_shared<SeqScalarProduct<X>>();
-  }
-  template<class M, class X, class Y, class C>
-  std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C> >& op)
-  {
-    return createScalarProduct<X>(op->getCommunication(), op->category());
-  }
-
-  template<class M, class X, class Y, class C>
-  std::shared_ptr<ScalarProduct<X>> createScalarProduct(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C> >& op)
-  {
-    return createScalarProduct<X>(op->getCommunication(), op->category());
-  }
-
  /**
-     @brief Factory to assembly solvers configured by a `ParameterTree`.
+     \brief Instantiates an `InverseOperator` from an Operator and a
+     configuration given as a ParameterTree.
+     \param op Operator
+     \param config `ParameterTree` with configuration

     Example ini File that can be passed in to construct a CGSolver with a SSOR
     preconditioner:
@@ -125,95 +143,51 @@ namespace Dune{
     \endverbatim

     \tparam Operator type of the operator, necessary to deduce the matrix type etc.
+
   */
  template<class Operator>
-  class SolverFactory {
-    using Domain = typename Operator::domain_type;
-    using Range = typename Operator::range_type;
-    using Solver = Dune::InverseOperator<Domain,Range>;
-    using Preconditioner = Dune::Preconditioner<Domain, Range>;
-
-    template<class O>
-    using _matrix_type = typename O::matrix_type;
-    using matrix_type = Std::detected_or_t<int, _matrix_type, Operator>;
-    static constexpr bool isAssembled = !std::is_same<matrix_type, int>::value;
-
-    static const matrix_type* getmat(std::shared_ptr<Operator> op){
-      std::shared_ptr<AssembledLinearOperator<matrix_type, Domain, Range>> aop
-        = std::dynamic_pointer_cast<AssembledLinearOperator<matrix_type, Domain, Range>>(op);
-      if(aop)
-        return &aop->getmat();
-      return nullptr;
-    }
-
-  public:
-
-    /** @brief get a solver from the factory
-     */
-    static std::shared_ptr<Solver> get(std::shared_ptr<Operator> op,
-                                       const ParameterTree& config,
-                                       std::shared_ptr<Preconditioner> prec = nullptr){
-      std::string type = config.get<std::string>("type");
-      std::shared_ptr<Solver> result;
-      const matrix_type* mat = getmat(op);
-      if(mat){
-        if (DirectSolverFactory<matrix_type, Domain, Range>::instance().contains(type)) {
-          if(op->category()!=SolverCategory::sequential){
-            DUNE_THROW(NotImplemented, "The solver factory does not support parallel direct solvers!");
-          }
-          result = DirectSolverFactory<matrix_type, Domain, Range>::instance().create(type, *mat, config);
-          return result;
-        }
-      }
-      // if no direct solver is found it might be an iterative solver
-      if (!IterativeSolverFactory<Domain, Range>::instance().contains(type)) {
-        DUNE_THROW(Dune::InvalidStateException, "Solver not found in the factory.");
-      }
-      if(!prec){
-        const ParameterTree& precConfig = config.sub("preconditioner");
-        std::string prec_type = precConfig.get<std::string>("type");
-        prec = PreconditionerFactory<Operator, Domain, Range>::instance().create(prec_type, op, precConfig);
-        if (prec->category() != op->category() && prec->category() == SolverCategory::sequential)
-          // try to wrap to a parallel preconditioner
-          prec = wrapPreconditioner4Parallel(prec, op);
-      }
-      std::shared_ptr<ScalarProduct<Domain>> sp = createScalarProduct(op);
-      result = IterativeSolverFactory<Domain, Range>::instance().create(type, op, sp, prec, config);
-      return result;
-    }
-
-    /**
-      @brief Construct a Preconditioner for a given Operator
-     */
-    static std::shared_ptr<Preconditioner> getPreconditioner(std::shared_ptr<Operator> op,
-                                                             const ParameterTree& config){
-      const matrix_type* mat = getmat(op);
-      if(mat){
-        std::string prec_type = config.get<std::string>("type");
-        return PreconditionerFactory<Operator, Domain, Range>::instance().create(prec_type, op, config);
-      }else{
-        DUNE_THROW(InvalidStateException, "Could not obtain matrix from operator. Please pass in an AssembledLinearOperator.");
-      }
+  std::shared_ptr<InverseOperator<typename Operator::domain_type,
+                                  typename Operator::range_type>>
+  getSolverFromFactory(std::shared_ptr<Operator> op, const ParameterTree& config,
+                       std::shared_ptr<Preconditioner<typename Operator::domain_type, typename Operator::range_type>> prec = nullptr)
+  {
+    if(prec){
+      PreconditionerFactory<Operator>::instance().define("__passed at runtime__",
+                                                         [=](auto...){
+                                                           return prec;
+                                                         });
+      ParameterTree config_tmp = config;
+      config_tmp.sub("preconditioner")["type"] = std::string("__passed at runtime__");
+      return SolverFactory<Operator>::instance().
+        create(config.get<std::string>("type"),op, config_tmp);
    }
-  };
+    return SolverFactory<Operator>::instance().
+      create(config.get<std::string>("type"),op, config);
+  }

+  class UnknownSolverCategory : public InvalidStateException{};
  /**
-     \brief Instantiates an `InverseOperator` from an Operator and a
-     configuration given as a ParameterTree.
-     \param op Operator
-     \param config `ParameterTree` with configuration
-     \param prec Custom `Preconditioner` (optional). If not given it will be
-     created with the `PreconditionerFactory` and the configuration given in
-     subKey "preconditioner".
-
-   */
+     @brief Construct a Preconditioner for a given Operator
+  */
  template<class Operator>
-  std::shared_ptr<InverseOperator<typename Operator::domain_type,
-                                  typename Operator::range_type>> getSolverFromFactory(std::shared_ptr<Operator> op,
-                               const ParameterTree& config,
-                               std::shared_ptr<Preconditioner<typename Operator::domain_type,
-                               typename Operator::range_type>> prec = nullptr){
-    return SolverFactory<Operator>::get(op, config, prec);
+  std::shared_ptr<Preconditioner<typename Operator::domain_type,
+                                 typename Operator::range_type>>
+  getPreconditionerFromFactory(std::shared_ptr<Operator> op,
+                               const ParameterTree& config){
+    using Domain = typename Operator::domain_type;
+    using Range = typename Operator::range_type;
+    std::string prec_type = config.get<std::string>("type");
+    std::shared_ptr<Preconditioner<typename Operator::domain_type,
+                                 typename Operator::range_type>> prec = PreconditionerFactory<Operator>::instance().create(prec_type, op, config);
+    if constexpr (OperatorTraits<Operator>::isParallel){
+      using Comm = typename OperatorTraits<Operator>::comm_type;
+      const Comm& comm = OperatorTraits<Operator>::getCommOrThrow(op);
+      if(op->category() == SolverCategory::overlapping && prec->category() == SolverCategory::sequential)
+        return std::make_shared<BlockPreconditioner<Domain,Range,Comm> >(prec, comm);
+      else if(op->category() == SolverCategory::nonoverlapping && prec->category() == SolverCategory::sequential)
+        return std::make_shared<NonoverlappingBlockPreconditioner<Comm, Preconditioner<Domain, Range>> >(prec, comm);
+    }
+    return prec;
  }

  /**

--- a/dune/istl/solverregistry.hh
+++ b/dune/istl/solverregistry.hh
@@ -10,65 +10,105 @@
 #include <dune/istl/preconditioner.hh>
 #include <dune/istl/solver.hh>

-#define DUNE_REGISTER_DIRECT_SOLVER(name, ...)                \
-  DUNE_REGISTRY_PUT(DirectSolverTag, name, __VA_ARGS__)
+namespace Dune::Impl {
+  template<class C>
+  [[deprecated("DUNE_REGISTER_ITERATIVE_SOLVER is deprecated. Please use DUNE_REGISTER_SOLVER instead.")]]
+  auto translateToOldIterativeSolverInterface(C oldCreator){
+    return [=](auto optraits, auto&&... args){
+      using OpTraits = decltype(optraits);
+      using TL = TypeList<typename OpTraits::domain_type, typename OpTraits::range_type>;
+      return oldCreator(TL{}, args...);
+    };
+  }
+}
+
+#define DUNE_REGISTER_DIRECT_SOLVER(name, ...)  \
+  _Pragma ("GCC warning \"'DUNE_REGISTER_DIRECT_SOLVER' macro is deprecated and will be removed after the release of DUNE 2.9. Please use 'DUNE_REGISTER_SOLVER'\"") \
+  DUNE_REGISTER_SOLVER(name, __VA_ARGS__);
+
+#define DUNE_REGISTER_ITERATIVE_SOLVER(name, ...)  \
+  _Pragma ("GCC warning \"'DUNE_REGISTER_ITERATIVE_SOLVER' macro is deprecated and will be removed after the release of DUNE 2.9. Please use 'DUNE_REGISTER_SOLVER'\"") \
+  DUNE_REGISTER_SOLVER(name, Impl::translateToOldIterativeSolverInterface(__VA_ARGS__));

 #define DUNE_REGISTER_PRECONDITIONER(name, ...)                \
  DUNE_REGISTRY_PUT(PreconditionerTag, name, __VA_ARGS__)

-#define DUNE_REGISTER_ITERATIVE_SOLVER(name, ...)                \
-  DUNE_REGISTRY_PUT(IterativeSolverTag, name, __VA_ARGS__)
+#define DUNE_REGISTER_SOLVER(name, ...)                \
+  DUNE_REGISTRY_PUT(SolverTag, name, __VA_ARGS__)

 namespace Dune{
+
  /** @addtogroup ISTL_Factory
      @{
  */

  namespace {
-    struct DirectSolverTag {};
    struct PreconditionerTag {};
-    struct IterativeSolverTag {};
+    struct SolverTag {};
  }
+
+  //! This exception is thrown if the requested solver or preconditioner needs an assembled matrix
+  class NoAssembledOperator : public InvalidStateException{};
+
  template<template<class,class,class,int>class Preconditioner, int blockLevel=1>
  auto defaultPreconditionerBlockLevelCreator(){
-    return [](auto typeList, const auto& matrix, const Dune::ParameterTree& config)
-           {
-             using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type;
-             using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type;
-             using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type;
-             std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
-               = std::make_shared<Preconditioner<Matrix, Domain, Range, blockLevel>>(matrix, config);
-             return preconditioner;
-           };
+    return [](auto opInfo, const auto& linearOperator, const Dune::ParameterTree& config)
+    {
+      using OpInfo = std::decay_t<decltype(opInfo)>;
+      using Matrix = typename OpInfo::matrix_type;
+      using Domain = typename OpInfo::domain_type;
+      using Range = typename OpInfo::range_type;
+      std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner;
+      if constexpr (OpInfo::isAssembled){
+        const auto& A = opInfo.getAssembledOpOrThrow(linearOperator);
+        // const Matrix& matrix = A->getmat();
+        preconditioner
+          = std::make_shared<Preconditioner<Matrix, Domain, Range, blockLevel>>(A, config);
+      }else{
+        DUNE_THROW(NoAssembledOperator, "Could not obtain matrix from operator. Please pass in an AssembledLinearOperator.");
+      }
+      return preconditioner;
+    };
  }

  template<template<class,class,class>class Preconditioner>
  auto defaultPreconditionerCreator(){
-    return [](auto typeList, const auto& matrix, const Dune::ParameterTree& config)
-           {
-             using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type;
-             using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type;
-             using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type;
-             std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
-               = std::make_shared<Preconditioner<Matrix, Domain, Range>>(matrix, config);
-             return preconditioner;
-           };
+    return [](auto opInfo, const auto& linearOperator, const Dune::ParameterTree& config)
+    {
+      using OpInfo = std::decay_t<decltype(opInfo)>;
+      using Matrix = typename OpInfo::matrix_type;
+      using Domain = typename OpInfo::domain_type;
+      using Range = typename OpInfo::range_type;
+      std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner;
+      if constexpr (OpInfo::isAssembled){
+        const auto& A = opInfo.getAssembledOpOrThrow(linearOperator);
+        // const Matrix& matrix = A->getmat();
+        preconditioner
+          = std::make_shared<Preconditioner<Matrix, Domain, Range>>(A, config);
+      }else{
+        DUNE_THROW(NoAssembledOperator, "Could not obtain matrix from operator. Please pass in an AssembledLinearOperator.");
+      }
+      return preconditioner;
+    };
  }

  template<template<class...>class Solver>
  auto defaultIterativeSolverCreator(){
-    return [](auto typeList,
+    return [](auto opInfo,
              const auto& linearOperator,
-              const auto& scalarProduct,
-              const auto& preconditioner,
              const Dune::ParameterTree& config)
-           {
-             using Domain = typename Dune::TypeListElement<0, decltype(typeList)>::type;
-             using Range = typename Dune::TypeListElement<1, decltype(typeList)>::type;
-             std::shared_ptr<Dune::InverseOperator<Domain, Range>> solver
-               = std::make_shared<Solver<Domain>>(linearOperator, scalarProduct, preconditioner, config);
-             return solver;
-           };
+    {
+      using OpInfo = std::decay_t<decltype(opInfo)>;
+      using Operator = typename OpInfo::operator_type;
+      using Domain = typename OpInfo::domain_type;
+      using Range = typename OpInfo::range_type;
+      std::shared_ptr<Operator> _op = std::dynamic_pointer_cast<Operator>(linearOperator);
+      std::shared_ptr<Preconditioner<Domain,Range>> preconditioner = getPreconditionerFromFactory(_op, config.sub("preconditioner"));
+      std::shared_ptr<ScalarProduct<Range>> scalarProduct = opInfo.getScalarProduct(linearOperator);
+      std::shared_ptr<Dune::InverseOperator<Domain, Range>> solver
+        = std::make_shared<Solver<Domain>>(linearOperator, scalarProduct, preconditioner, config);
+      return solver;
+    };
  }

  /* This exception is thrown, when the requested solver is in the factory but

--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -70,7 +70,7 @@ namespace Dune {
    using IterativeSolver<X,X>::apply;

    //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      Iteration iteration(*this, res);
      _prec->pre(x,b);
@@ -115,7 +115,7 @@ namespace Dune {
    using IterativeSolver<X,X>::_verbose;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("loopsolver", defaultIterativeSolverCreator<Dune::LoopSolver>());
+  DUNE_REGISTER_SOLVER("loopsolver", defaultIterativeSolverCreator<Dune::LoopSolver>());


  // all these solvers are taken from the SUMO library
@@ -139,7 +139,7 @@ namespace Dune {

       \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
     */
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      Iteration iteration(*this, res);
      _prec->pre(x,b);             // prepare preconditioner
@@ -186,7 +186,7 @@ namespace Dune {
    using IterativeSolver<X,X>::_verbose;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("gradientsolver", defaultIterativeSolverCreator<Dune::GradientSolver>());
+  DUNE_REGISTER_SOLVER("gradientsolver", defaultIterativeSolverCreator<Dune::GradientSolver>());

  //! \brief conjugate gradient method
  template<class X>
@@ -276,7 +276,7 @@ namespace Dune {
             E.g. numeric_limits<double>::quiet_NaN()*0.0==0.0 with gcc-5.3
             -ffast-math.
     */
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      Iteration iteration(*this,res);
      _prec->pre(x,b);             // prepare preconditioner
@@ -411,7 +411,7 @@ namespace Dune {
    using IterativeSolver<X,X>::_verbose;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("cgsolver", defaultIterativeSolverCreator<Dune::CGSolver>());
+  DUNE_REGISTER_SOLVER("cgsolver", defaultIterativeSolverCreator<Dune::CGSolver>());

  // Ronald Kriemanns BiCG-STAB implementation from Sumo
  //! \brief Bi-conjugate Gradient Stabilized (BiCG-STAB)
@@ -436,7 +436,7 @@ namespace Dune {

       \note Currently, the BiCGSTABSolver aborts when it detects a breakdown.
     */
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      using std::abs;
      const Simd::Scalar<real_type> EPSILON=1e-80;
@@ -597,7 +597,7 @@ namespace Dune {
    template<class CountType>
    using Iteration = typename IterativeSolver<X,X>::template Iteration<CountType>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("bicgstabsolver", defaultIterativeSolverCreator<Dune::BiCGSTABSolver>());
+  DUNE_REGISTER_SOLVER("bicgstabsolver", defaultIterativeSolverCreator<Dune::BiCGSTABSolver>());

  /*! \brief Minimal Residual Method (MINRES)

@@ -624,7 +624,7 @@ namespace Dune {

       \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
     */
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      using std::sqrt;
      using std::abs;
@@ -807,7 +807,7 @@ namespace Dune {
    using IterativeSolver<X,X>::_verbose;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("minressolver", defaultIterativeSolverCreator<Dune::MINRESSolver>());
+  DUNE_REGISTER_SOLVER("minressolver", defaultIterativeSolverCreator<Dune::MINRESSolver>());

  /**
     \brief implements the Generalized Minimal Residual (GMRes) method
@@ -907,7 +907,7 @@ namespace Dune {
       \note Currently, the RestartedGMResSolver aborts when it detects a
             breakdown.
     */
-    virtual void apply (X& x, Y& b, InverseOperatorResult& res)
+    void apply (X& x, Y& b, InverseOperatorResult& res) override
    {
      apply(x,b,Simd::max(_reduction),res);
    }
@@ -920,7 +920,7 @@ namespace Dune {
       \note Currently, the RestartedGMResSolver aborts when it detects a
             breakdown.
     */
-    virtual void apply (X& x, Y& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+    void apply (X& x, Y& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
    {
      using std::abs;
      const Simd::Scalar<real_type> EPSILON = 1e-80;
@@ -1119,7 +1119,7 @@ namespace Dune {
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
    int _restart;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("restartedgmressolver", defaultIterativeSolverCreator<Dune::RestartedGMResSolver>());
+  DUNE_REGISTER_SOLVER("restartedgmressolver", defaultIterativeSolverCreator<Dune::RestartedGMResSolver>());

  /**
     \brief implements the Flexible Generalized Minimal Residual (FGMRes) method (right preconditioned)
@@ -1287,7 +1287,7 @@ private:
    using RestartedGMResSolver<X,Y>::_restart;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("restartedflexiblegmressolver", defaultIterativeSolverCreator<Dune::RestartedFlexibleGMResSolver>());
+  DUNE_REGISTER_SOLVER("restartedflexiblegmressolver", defaultIterativeSolverCreator<Dune::RestartedFlexibleGMResSolver>());

  /**
   * @brief Generalized preconditioned conjugate gradient solver.
@@ -1388,7 +1388,7 @@ private:

       \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
     */
-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      Iteration iteration(*this, res);
      _prec->pre(x,b);                 // prepare preconditioner
@@ -1483,7 +1483,7 @@ private:
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
    int _restart;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("generalizedpcgsolver", defaultIterativeSolverCreator<Dune::GeneralizedPCGSolver>());
+  DUNE_REGISTER_SOLVER("generalizedpcgsolver", defaultIterativeSolverCreator<Dune::GeneralizedPCGSolver>());

  /*! \brief Accelerated flexible conjugate gradient method

@@ -1580,7 +1580,7 @@ private:
             -ffast-math.
     */

-    virtual void apply (X& x, X& b, InverseOperatorResult& res)
+    void apply (X& x, X& b, InverseOperatorResult& res) override
    {
      using rAlloc = ReboundAllocatorType<X,field_type>;
      res.clear();
@@ -1667,7 +1667,7 @@ private:
    using IterativeSolver<X,X>::_verbose;
    using Iteration = typename IterativeSolver<X,X>::template Iteration<unsigned int>;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("restartedfcgsolver", defaultIterativeSolverCreator<Dune::RestartedFCGSolver>());
+  DUNE_REGISTER_SOLVER("restartedfcgsolver", defaultIterativeSolverCreator<Dune::RestartedFCGSolver>());

  /*! \brief Complete flexible conjugate gradient method

@@ -1690,7 +1690,8 @@ private:
    using RestartedFCGSolver<X>::apply;

    // just a minor part of the RestartedFCGSolver apply method will be modified
-    virtual void apply (X& x, X& b, InverseOperatorResult& res) override {
+    void apply (X& x, X& b, InverseOperatorResult& res) override
+    {
      // reset limiter of orthogonalization loop
      _k_limit = 0;
      this->RestartedFCGSolver<X>::apply(x,b,res);
@@ -1698,7 +1699,7 @@ private:

  private:
    // This function is called every iteration to orthogonalize against the last search directions.
-    virtual void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<field_type,ReboundAllocatorType<X,field_type>>& ddotAd,std::vector<X>& d) override {
+    void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<field_type,ReboundAllocatorType<X,field_type>>& ddotAd,std::vector<X>& d) override {
      // This FCGSolver uses values with higher array indexes too, if existent.
      for (int k = 0; k < _k_limit; k++) {
        if(i_bounded!=k)
@@ -1711,7 +1712,7 @@ private:
    };

    // This function is called every mmax iterations to handle limited array sizes.
-    virtual void cycle(std::vector<X>& Ad, [[maybe_unused]] std::vector<X>& d, [[maybe_unused]] std::vector<field_type,ReboundAllocatorType<X,field_type> >& ddotAd,int& i_bounded) override {
+    void cycle(std::vector<X>& Ad, [[maybe_unused]] std::vector<X>& d, [[maybe_unused]] std::vector<field_type,ReboundAllocatorType<X,field_type> >& ddotAd,int& i_bounded) override {
      // Only the loop index i_bounded return to 0, if it reached mmax.
      i_bounded = 0;
      // Now all arrays are filled and the loop in void orthogonalizations can use the whole arrays.
@@ -1729,7 +1730,7 @@ private:
    using RestartedFCGSolver<X>::_maxit;
    using RestartedFCGSolver<X>::_verbose;
  };
-  DUNE_REGISTER_ITERATIVE_SOLVER("completefcgsolver", defaultIterativeSolverCreator<Dune::CompleteFCGSolver>());
+  DUNE_REGISTER_SOLVER("completefcgsolver", defaultIterativeSolverCreator<Dune::CompleteFCGSolver>());
  /** @} end documentation */
 } // end namespace


--- a/dune/istl/spqr.hh
+++ b/dune/istl/spqr.hh
@@ -79,7 +79,7 @@ namespace Dune {
    typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;

    //! Category of the solver (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::Category::sequential;
    }
@@ -152,7 +152,7 @@ namespace Dune {
    }

    /** \copydoc InverseOperator::apply(X&, Y&, InverseOperatorResult&) */
-    virtual void apply(domain_type& x, range_type& b, InverseOperatorResult& res)
+    void apply(domain_type& x, range_type& b, InverseOperatorResult& res) override
    {
      const std::size_t numRows(spqrMatrix_.N());
      // fill B
@@ -188,7 +188,7 @@ namespace Dune {
    }

    /** \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) */
-    virtual void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+    void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
    {
      apply(x, b, res);
    }
@@ -330,35 +330,34 @@ namespace Dune {
    enum {value = true};
  };

-  struct SPQRCreator {
-    template<class> struct isValidBlock : std::false_type{};
-
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
-      std::enable_if_t<
-                isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      int verbose = config.get("verbose", 0);
-      return std::make_shared<Dune::SPQR<M>>(mat,verbose);
-    }
-
-    // second version with SFINAE to validate the template parameters of SPQR
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
-      std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      DUNE_THROW(UnsupportedType,
-        "Unsupported Type in SPQR (only double and std::complex<double> supported)");
-    }
-  };
-  template<> struct SPQRCreator::isValidBlock<FieldVector<double,1>> : std::true_type{};
-  // std::complex is temporary disabled, because it fails if libc++ is used
-  //template<> struct SPQRCreator::isValidMatrixBlock<FieldMatrix<std::complex<double>,1,1>> : std::true_type{};
-  DUNE_REGISTER_DIRECT_SOLVER("spqr", Dune::SPQRCreator());
+  DUNE_REGISTER_SOLVER("spqr",
+                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
+                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
+                       {
+                         using OpTraits = decltype(opTraits);
+                         using M = typename OpTraits::matrix_type;
+                         // works only for sequential operators
+                         if constexpr (OpTraits::isParallel){
+                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
+                             DUNE_THROW(Dune::InvalidStateException, "SPQR works only for sequential operators.");
+                         }
+                         // check if SPQR<M>* is convertible to
+                         // InverseOperator*. This allows only the explicit
+                         // specialized variants of SPQR
+                         // std::complex is temporary disabled, because it fails if libc++ is used
+                         if constexpr (std::is_convertible_v<SPQR<M>*,
+                                       Dune::InverseOperator<typename OpTraits::domain_type,
+                                       typename OpTraits::range_type>*> &&
+                                       std::is_same_v<typename FieldTraits<M>::field_type, double>
+                                       ){
+                           const auto& A = opTraits.getAssembledOpOrThrow(op);
+                           const M& mat = A->getmat();
+                           int verbose = config.get("verbose", 0);
+                           return std::make_shared<Dune::SPQR<M>>(mat,verbose);
+                         }
+                         DUNE_THROW(UnsupportedType,
+                                    "Unsupported Type in SPQR (only FieldMatrix<double,...> supported)");
+                       });

 } // end namespace Dune


--- a/dune/istl/superlu.hh
+++ b/dune/istl/superlu.hh
@@ -284,7 +284,7 @@ namespace Dune
    using range_type = typename Impl::SuperLUVectorChooser<M>::range_type;

    //! Category of the solver (see SolverCategory::Category)
-    virtual SolverCategory::Category category() const
+    SolverCategory::Category category() const override
    {
      return SolverCategory::Category::sequential;
    }
@@ -334,12 +334,12 @@ namespace Dune
    /**
     *  \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
     */
-    void apply(domain_type& x, range_type& b, InverseOperatorResult& res);
+    void apply(domain_type& x, range_type& b, InverseOperatorResult& res) override;

    /**
     *  \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
     */
-    void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+    void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
    {
      apply(x,b,res);
    }
@@ -724,35 +724,37 @@ namespace Dune
    enum { value = true };
  };

-  struct SuperLUCreator {
-    template<class> struct isValidBlock : std::false_type{};
-    template<int k> struct isValidBlock<Dune::FieldVector<double,k>> : std::true_type{};
-    template<int k> struct isValidBlock<Dune::FieldVector<std::complex<double>,k>> : std::true_type{};
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
-                std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      int verbose = config.get("verbose", 0);
-      return std::make_shared<Dune::SuperLU<M>>(mat,verbose);
-    }
-
-    // second version with SFINAE to validate the template parameters of SuperLU
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
-      std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      DUNE_THROW(UnsupportedType,
-        "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
-    }
-  };
-  template<> struct SuperLUCreator::isValidBlock<double> : std::true_type{};
-  template<> struct SuperLUCreator::isValidBlock<std::complex<double>> : std::true_type{};

-  DUNE_REGISTER_DIRECT_SOLVER("superlu", SuperLUCreator());
+  DUNE_REGISTER_SOLVER("superlu",
+                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
+                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
+                       {
+                         using OpTraits = decltype(opTraits);
+                         using M = typename OpTraits::matrix_type;
+                         // works only for sequential operators
+                         if constexpr (OpTraits::isParallel){
+                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
+                             DUNE_THROW(Dune::InvalidStateException, "SuperLU works only for sequential operators.");
+                         }
+                         if constexpr (OpTraits::isAssembled &&
+                                       // check whether the Matrix field_type is double or complex<double>
+                                       std::is_same_v<typename FieldTraits<M>::real_type, double>){
+                           // check if SuperLU<M>* is convertible to
+                           // InverseOperator*. This checks compatibility of the
+                           // domain and range types
+                           if constexpr (std::is_convertible_v<SuperLU<M>*,
+                                         Dune::InverseOperator<typename OpTraits::domain_type,
+                                         typename OpTraits::range_type>*>
+                                         ){
+                             const auto& A = opTraits.getAssembledOpOrThrow(op);
+                             const M& mat = A->getmat();
+                             int verbose = config.get("verbose", 0);
+                             return std::make_shared<Dune::SuperLU<M>>(mat,verbose);
+                           }
+                         }
+                         DUNE_THROW(UnsupportedType,
+                                    "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
+                       });
 } // end namespace DUNE

 // undefine macros from SuperLU's slu_util.h

--- a/dune/istl/supermatrix.hh
+++ b/dune/istl/supermatrix.hh
@@ -300,7 +300,7 @@ namespace Dune
    }

    /** @brief free allocated space. */
-    virtual void free()
+    void free() override
    {
      ISTL::Impl::BCCSMatrix<typename BCRSMatrix<B,TA>::field_type, int>::free();
      SUPERLU_FREE(A.Store);
@@ -326,7 +326,7 @@ namespace Dune
    SuperMatrixInitializer() : ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<B,A>, int>()
    {}

-    virtual void createMatrix() const
+    void createMatrix() const override
    {
      ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<B,A>, int>::createMatrix();
      SuperMatrixCreateSparseChooser<typename Matrix::field_type>

--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -158,6 +158,11 @@ dune_add_test(SOURCES matrixmarkettest.cc)

 dune_add_test(SOURCES iluildltest.cc)

+dune_add_test(SOURCES matrixfreesolverfactorytest.cc)
+
 dune_add_test(SOURCES blocklevel.cc COMPILE_ONLY)

+dune_add_test(SOURCES registermacrosbackwardcompatibilitytest.cc
+              CMAKE_GUARD "SuiteSparse_UMFPACK_FOUND")
+
 exclude_from_headercheck(complexdata.hh)
No results found