3 #ifndef DUNE_ISTL_PRECONDITIONERS_HH
4 #define DUNE_ISTL_PRECONDITIONERS_HH
13 #include <dune/common/simd/simd.hh>
14 #include <dune/common/unused.hh>
15 #include <dune/common/parametertree.hh>
71 template<
class O,
int c = -1>
73 public Preconditioner<typename O::domain_type, typename O::range_type>
92 : inverse_operator_(inverse_operator)
95 DUNE_THROW(InvalidStateException,
"User supplied solver category does not match that of the supplied iverser operator");
105 inverse_operator_.apply(v, copy, res);
137 template<
class M,
class X,
class Y,
int l=1>
159 : _A_(A), _n(n), _w(w)
177 SeqSSOR (
const M& A,
const ParameterTree& configuration)
180 _n = configuration.get<
int>(
"iterations",1);
190 virtual void pre (X& x, Y& b)
192 DUNE_UNUSED_PARAMETER(x);
193 DUNE_UNUSED_PARAMETER(b);
202 virtual void apply (X& v,
const Y& d)
204 for (
int i=0; i<_n; i++) {
217 DUNE_UNUSED_PARAMETER(x);
248 template<
class M,
class X,
class Y,
int l=1>
270 : _A_(A), _n(n), _w(w)
288 SeqSOR (
const M& A,
const ParameterTree& configuration)
291 _n = configuration.get<
int>(
"iterations",1);
301 virtual void pre (X& x, Y& b)
303 DUNE_UNUSED_PARAMETER(x);
304 DUNE_UNUSED_PARAMETER(b);
312 virtual void apply (X& v,
const Y& d)
314 this->
template apply<true>(v,d);
325 template<
bool forward>
329 for (
int i=0; i<_n; i++) {
333 for (
int i=0; i<_n; i++) {
345 DUNE_UNUSED_PARAMETER(x);
375 template<
class M,
class X,
class Y,
int l=1>
389 template<
class M,
class X,
class Y,
int l=1>
411 : _A_(A), _n(n), _w(w)
429 SeqJac (
const M& A,
const ParameterTree& configuration)
432 _n = configuration.get<
int>(
"iterations",1);
442 virtual void pre (X& x, Y& b)
444 DUNE_UNUSED_PARAMETER(x);
445 DUNE_UNUSED_PARAMETER(b);
453 virtual void apply (X& v,
const Y& d)
455 for (
int i=0; i<_n; i++) {
467 DUNE_UNUSED_PARAMETER(x);
482 scalar_field_type _w;
499 template<
class M,
class X,
class Y,
int l=1>
528 :
SeqILU( A, 0, w, resort )
546 SeqILU(
const M& A,
const ParameterTree& config)
549 config.
get(
"resort", false))
596 virtual void pre (X& x, Y& b)
598 DUNE_UNUSED_PARAMETER(x);
599 DUNE_UNUSED_PARAMETER(b);
607 virtual void apply (X& v,
const Y& d)
631 DUNE_UNUSED_PARAMETER(x);
642 std::unique_ptr< matrix_type >
ILU_;
647 std::vector< block_type, typename matrix_type::allocator_type >
inv_;
669 template<
class M,
class X,
class Y,
int l=1>
709 SeqILU0 (
const M& A,
const ParameterTree& configuration)
721 virtual void pre (X& x, Y& b)
723 DUNE_UNUSED_PARAMETER(x);
724 DUNE_UNUSED_PARAMETER(b);
732 virtual void apply (X& v,
const Y& d)
745 DUNE_UNUSED_PARAMETER(x);
756 scalar_field_type _w;
777 template<
class M,
class X,
class Y,
int l=1>
799 : ILU(A.N(),A.M(),M::row_wise),
809 SeqILUn (
const M& A,
const ParameterTree& configuration)
810 : ILU(A.N(),A.M(),M::row_wise)
812 _n = configuration.get<
int>(
"iterations");
822 virtual void pre (X& x, Y& b)
824 DUNE_UNUSED_PARAMETER(x);
825 DUNE_UNUSED_PARAMETER(b);
833 virtual void apply (X& v,
const Y& d)
846 DUNE_UNUSED_PARAMETER(x);
861 scalar_field_type _w;
874 template<
class X,
class Y>
916 virtual void pre (X& x, Y& b)
918 DUNE_UNUSED_PARAMETER(x);
919 DUNE_UNUSED_PARAMETER(b);
927 virtual void apply (X& v,
const Y& d)
940 DUNE_UNUSED_PARAMETER(x);
954 using D =
typename Dune::TypeListElement<1, decltype(tl)>::type;
955 using R =
typename Dune::TypeListElement<2, decltype(tl)>::type;
956 return std::make_shared<Richardson<D,R>>(config);
970 template<
class M,
class X,
class Y >
1014 : decomposition_( A.N(), A.M(),
matrix_type::random ),
1018 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
1020 const auto &A_i = *i;
1021 const auto ij = A_i.find( i.index() );
1022 if( ij != A_i.end() )
1023 decomposition_.setrowsize( i.index(), ij.offset()+1 );
1025 DUNE_THROW(
ISTLError,
"diagonal entry missing" );
1027 decomposition_.endrowsizes();
1030 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
1032 const auto &A_i = *i;
1033 for(
auto ij = A_i.begin(); ij.index() < i.index() ; ++ij )
1034 decomposition_.addindex( i.index(), ij.index() );
1035 decomposition_.addindex( i.index(), i.index() );
1037 decomposition_.endindices();
1041 for(
auto row = decomposition_.begin(), rowend = decomposition_.end(); row != rowend; ++row, ++i )
1043 auto ij = i->begin();
1044 for(
auto col = row->begin(), colend = row->end();
col != colend; ++
col, ++ij )
1053 void pre ( X &x, Y &b )
override
1055 DUNE_UNUSED_PARAMETER( x );
1056 DUNE_UNUSED_PARAMETER( b );
1060 void apply ( X &v,
const Y &d )
override
1069 DUNE_UNUSED_PARAMETER( x );
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Incomplete LDL decomposition.
Some handy generic functions for ISTL matrices.
Define general, extensible interface for inverse operators.
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:640
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:652
void bilu_backsolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:95
void bilu_decomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:169
void bilu0_decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:34
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:628
Definition: allocator.hh:7
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition: ildl.hh:86
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:292
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
void bildl_backsolve(const Matrix &A, X &v, const Y &d, bool isLowerTriangular=false)
Definition: ildl.hh:147
void convertToCRS(const M &A, CRS &lower, CRS &upper, InvVector &inv)
convert ILU decomposition into CRS format for lower and upper triangular and inverse.
Definition: ilu.hh:310
void bilu_backsolve(const CRS &lower, const CRS &upper, const InvVector &inv, X &v, const Y &d)
LU backsolve with stored inverse in CRS format for lower and upper triangular.
Definition: ilu.hh:377
compile-time parameter for block recursion depth
Definition: gsetc.hh:43
derive error class from the base class in common
Definition: istlexception.hh:16
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:52
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Turns an InverseOperator into a Preconditioner.
Definition: preconditioners.hh:74
O::range_type range_type
The range type of the preconditioner.
Definition: preconditioners.hh:79
O::domain_type domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:77
virtual void post(domain_type &)
Clean up.
Definition: preconditioners.hh:108
range_type::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:81
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:83
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:112
virtual void pre(domain_type &, range_type &)
Prepare the preconditioner.
Definition: preconditioners.hh:98
InverseOperator2Preconditioner(InverseOperator &inverse_operator)
Construct the preconditioner from the solver.
Definition: preconditioners.hh:91
O InverseOperator
type of the wrapped inverse operator
Definition: preconditioners.hh:85
virtual void apply(domain_type &v, const range_type &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: preconditioners.hh:101
Sequential SSOR preconditioner.
Definition: preconditioners.hh:138
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:215
SeqSSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:177
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:221
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:147
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:149
SeqSSOR(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:158
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:143
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:141
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:202
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:190
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:145
Sequential SOR preconditioner.
Definition: preconditioners.hh:249
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:252
void apply(X &v, const Y &d)
Apply the preconditioner in a special direction.
Definition: preconditioners.hh:326
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:254
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:343
SeqSOR(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:269
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:301
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:260
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:349
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:312
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:256
SeqSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:288
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:258
The sequential jacobian preconditioner.
Definition: preconditioners.hh:390
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:465
SeqJac(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:429
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:453
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:393
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:401
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:399
SeqJac(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:410
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:442
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:395
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:471
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:397
Sequential ILU preconditioner.
Definition: preconditioners.hh:500
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:629
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:596
matrix_type ::block_type block_type
block type of matrix
Definition: preconditioners.hh:505
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:607
ILU::CRS< block_type, typename M::allocator_type > CRS
type of ILU storage
Definition: preconditioners.hh:518
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:509
const scalar_field_type w_
The relaxation factor to use.
Definition: preconditioners.hh:650
CRS lower_
The ILU(n) decomposition of the matrix. As storage a CRS structure is used.
Definition: preconditioners.hh:645
const bool wNotIdentity_
true if w != 1.0
Definition: preconditioners.hh:652
SeqILU(const M &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:546
SeqILU(const M &A, int n, scalar_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:560
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:503
SeqILU(const M &A, scalar_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:527
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:512
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:635
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:507
std::vector< block_type, typename matrix_type::allocator_type > inv_
Definition: preconditioners.hh:647
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:515
std::unique_ptr< matrix_type > ILU_
The ILU(n) decomposition of the matrix. As storage a BCRSMatrix is used.
Definition: preconditioners.hh:642
CRS upper_
Definition: preconditioners.hh:646
Sequential ILU0 preconditioner.
Definition: preconditioners.hh:670
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:743
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:681
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:721
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:749
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:732
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:677
SeqILU0(const M &A, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:689
SeqILU0(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:709
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:675
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:679
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:673
Sequential ILU(n) preconditioner.
Definition: preconditioners.hh:778
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:789
SeqILUn(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:798
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:822
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:781
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:844
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:783
SeqILUn(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:809
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:785
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: preconditioners.hh:833
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:787
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:850
Richardson preconditioner.
Definition: preconditioners.hh:875
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:882
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:944
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:880
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:916
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:938
Richardson(scalar_field_type w=1.0)
Constructor.
Definition: preconditioners.hh:891
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:884
Richardson(const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:906
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:878
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: preconditioners.hh:927
sequential ILDL preconditioner
Definition: preconditioners.hh:973
SeqILDL(const matrix_type &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:1001
SeqILDL(const matrix_type &A, scalar_field_type relax=scalar_field_type(1))
constructor
Definition: preconditioners.hh:1013
X domain_type
domain type of the preconditioner
Definition: preconditioners.hh:981
void post(X &x) override
Clean up.
Definition: preconditioners.hh:1067
Y range_type
range type of the preconditioner
Definition: preconditioners.hh:983
std::remove_const_t< M > matrix_type
type of matrix the preconditioner is for
Definition: preconditioners.hh:979
void apply(X &v, const Y &d) override
Apply one step of the preconditioner to the system A(v)=d.
Definition: preconditioners.hh:1060
void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: preconditioners.hh:1053
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:987
X::field_type field_type
field type of the preconditioner
Definition: preconditioners.hh:985
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:1073
Statistics about the application of an inverse operator.
Definition: solver.hh:46
Abstract base class for all solvers.
Definition: solver.hh:97
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:32