From b60572ac3ddfcbd83b4e793d5d8fc782f8548793 Mon Sep 17 00:00:00 2001
From: Nils-Arne Dreier <n.dreier@uni-muenster.de>
Date: Mon, 13 Sep 2021 14:49:46 +0200
Subject: [PATCH] simplify and cleanup solverfactory
---
dune/istl/cholmod.hh | 28 ++--
dune/istl/common/registry.hh | 2 +-
dune/istl/ldl.hh | 26 ++--
dune/istl/paamg/amg.hh | 28 ++--
dune/istl/preconditioners.hh | 11 +-
dune/istl/solverfactory.hh | 251 ++++++++++++++---------------------
dune/istl/solverregistry.hh | 75 ++++++-----
dune/istl/solvers.hh | 20 +--
dune/istl/spqr.hh | 25 ++--
dune/istl/superlu.hh | 25 ++--
dune/istl/umfpack.hh | 36 +++--
11 files changed, 250 insertions(+), 277 deletions(-)
diff --git a/dune/istl/cholmod.hh b/dune/istl/cholmod.hh
index 87c29eb69..836584f8d 100644
--- 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>
@@ -529,29 +529,31 @@ private:
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
+ template<class OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
- using D = typename Dune::TypeListElement<1, TL>::type;
+ using D = typename OpTraits::domain_type;
+ using M = typename OpTraits::matrix_type;
+ const M& mat = opTraits.getMatOrThrow(op);
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
+ template<class OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<!isValidBlock<typename OpTraits::matrix_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", Dune::CholmodCreator());
} /* namespace Dune */
diff --git a/dune/istl/common/registry.hh b/dune/istl/common/registry.hh
index 33ddee83d..54dc9a296 100644
--- 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.
*/
diff --git a/dune/istl/ldl.hh b/dune/istl/ldl.hh
index 0c2ed6bab..5e6de06b0 100644
--- 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 {
/**
@@ -374,30 +374,32 @@ namespace Dune {
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,
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
- isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
+ isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
+ using M = typename OpTraits::matrix_type;
+ const M& mat = opTraits.getMatOrThrow(op);
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*/,
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
- !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
+ !isValidBlock<typename OpTraits::matrix_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", Dune::LDLCreator());
} // end namespace Dune
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 8cb6aa374..dcac8e745 100644
--- 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>
@@ -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<isValidBlockType<typename OpTraits::matrix_type::block_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<!isValidBlockType<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
}
diff --git a/dune/istl/preconditioners.hh b/dune/istl/preconditioners.hh
index edacb143d..f62daa81d 100644
--- a/dune/istl/preconditioners.hh
+++ b/dune/istl/preconditioners.hh
@@ -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);
+ });
/**
diff --git a/dune/istl/solverfactory.hh b/dune/istl/solverfactory.hh
index ba05aed95..3f6a5364b 100644
--- 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,88 @@ 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>>>;
+
+ class NoAssembledOperator : public InvalidStateException{};
+
+ 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 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 matrix_type& getMatOrThrow(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->getmat();
+ DUNE_THROW(NoAssembledOperator, "Could not obtain matrix from operator. 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>>();
+ }
+ }
+ };
// 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 +124,39 @@ 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){
+ 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;
}
/**
diff --git a/dune/istl/solverregistry.hh b/dune/istl/solverregistry.hh
index 390df0de2..4d9d57bdd 100644
--- a/dune/istl/solverregistry.hh
+++ b/dune/istl/solverregistry.hh
@@ -10,65 +10,70 @@
#include <dune/istl/preconditioner.hh>
#include <dune/istl/solver.hh>
-#define DUNE_REGISTER_DIRECT_SOLVER(name, ...) \
- DUNE_REGISTRY_PUT(DirectSolverTag, name, __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 {};
}
+
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;
+ const Matrix& matrix = opInfo.getMatOrThrow(linearOperator);
+ std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
+ = std::make_shared<Preconditioner<Matrix, Domain, Range, blockLevel>>(matrix, config);
+ 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;
+ const Matrix& matrix = opInfo.getMatOrThrow(linearOperator);
+ std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
+ = std::make_shared<Preconditioner<Matrix, Domain, Range>>(matrix, config);
+ 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
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index 98dbc9fdf..cb09353ad 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -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
@@ -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>
@@ -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)
@@ -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)
@@ -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
@@ -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.
@@ -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
@@ -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
@@ -1730,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
diff --git a/dune/istl/spqr.hh b/dune/istl/spqr.hh
index 61c05510c..4e2d6b665 100644
--- a/dune/istl/spqr.hh
+++ b/dune/istl/spqr.hh
@@ -333,23 +333,26 @@ namespace Dune {
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,
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
- isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
+ isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
+ using M = typename OpTraits::matrix_type;
+ const M& mat = opTraits.getMatOrThrow(op);
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
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
+ !isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
DUNE_THROW(UnsupportedType,
"Unsupported Type in SPQR (only double and std::complex<double> supported)");
@@ -358,7 +361,7 @@ namespace Dune {
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", Dune::SPQRCreator());
} // end namespace Dune
diff --git a/dune/istl/superlu.hh b/dune/istl/superlu.hh
index 783e9d194..9ae8e6328 100644
--- a/dune/istl/superlu.hh
+++ b/dune/istl/superlu.hh
@@ -728,22 +728,25 @@ namespace Dune
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
+
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<typename OpTraits::domain_type,
+ typename OpTraits::range_type> >
+ operator() (OpTraits opTraits, const OP& op, const Dune::ParameterTree& config,
+ std::enable_if_t<isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
+ using M = typename OpTraits::matrix_type;
+ const M& mat = opTraits.getMatOrThrow(op);
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
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<typename OpTraits::domain_type,
+ typename OpTraits::range_type>>
+ operator() (OpTraits opTraits, const OP& op, const Dune::ParameterTree& config,
+ std::enable_if_t<!isValidBlock<typename OpTraits::matrix_type::block_type>::value,int> = 0) const
{
DUNE_THROW(UnsupportedType,
"Unsupported Type in SuperLU (only double and std::complex<double> supported)");
@@ -752,7 +755,7 @@ namespace Dune
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", SuperLUCreator());
} // end namespace DUNE
// undefine macros from SuperLU's slu_util.h
diff --git a/dune/istl/umfpack.hh b/dune/istl/umfpack.hh
index 54cb264b0..17bda13d8 100644
--- a/dune/istl/umfpack.hh
+++ b/dune/istl/umfpack.hh
@@ -796,28 +796,36 @@ namespace Dune {
struct UMFPackCreator {
- template<class TL, class M,class=void> struct isValidBlock : std::false_type{};
- template<class TL, class M> struct isValidBlock<TL,M,
+ template<class OpTraits, class=void> struct isValidBlock : std::false_type{};
+ template<class OpTraits> struct isValidBlock<OpTraits,
std::enable_if_t<
- std::is_same_v<Impl::UMFPackDomainType<M>, typename Dune::TypeListElement<1,TL>::type>
- && std::is_same_v<Impl::UMFPackRangeType<M>, typename Dune::TypeListElement<2,TL>::type>
+ std::is_same_v<Impl::UMFPackDomainType<typename OpTraits::matrix_type>, typename OpTraits::domain_type>
+ && std::is_same_v<Impl::UMFPackRangeType<typename OpTraits::matrix_type>, typename OpTraits::range_type>
>> : std::true_type {};
- template<typename TL, typename M>
- std::shared_ptr<Dune::InverseOperator<Impl::UMFPackDomainType<M>,Impl::UMFPackRangeType<M>>>
- operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
- std::enable_if_t<isValidBlock<TL, M>::value,int> = 0) const
+ template<typename OpTraits, typename OP>
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
+ isValidBlock<OpTraits>::value
+ && Simd::lanes<typename OpTraits::domain_type::field_type>() == 1,int> = 0) const
{
+ using M = typename OpTraits::matrix_type;
+ const M& mat = opTraits.getMatOrThrow(op);
int verbose = config.get("verbose", 0);
return std::make_shared<Dune::UMFPack<M>>(mat,verbose);
}
// second version with SFINAE to validate the template parameters of UMFPack
- 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<TL, M>::value,int> = 0) const
+ template<typename OpTraits, typename OP>
+ std::shared_ptr<Dune::InverseOperator<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<
+ !isValidBlock<OpTraits>::value
+ || Simd::lanes<typename OpTraits::domain_type::field_type>() != 1,int> = 0) const
{
using D = typename Dune::TypeListElement<1,TL>::type;
using R = typename Dune::TypeListElement<2,TL>::type;
@@ -833,7 +841,7 @@ namespace Dune {
);
}
};
- DUNE_REGISTER_DIRECT_SOLVER("umfpack",Dune::UMFPackCreator());
+ DUNE_REGISTER_SOLVER("umfpack",Dune::UMFPackCreator());
} // end namespace Dune
#endif // HAVE_SUITESPARSE_UMFPACK
--
GitLab