Go to the documentation of this file.
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>
670 class DUNE_DEPRECATED_MSG("Use
SeqILU instead of SeqILU0!") SeqILU0:
public Preconditioner<X,Y> {
673 typedef typename std::remove_const<M>::type matrix_type;
675 typedef X domain_type;
677 typedef Y range_type;
679 typedef typename X::field_type field_type;
681 typedef Simd::Scalar<field_type> scalar_field_type;
689 SeqILU0 (
const M& A, scalar_field_type w)
709 SeqILU0 (
const M& A,
const ParameterTree& configuration)
712 _w = configuration.get<scalar_field_type>(
"relaxation");
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)
743 virtual void post (X& x)
745 DUNE_UNUSED_PARAMETER(x);
756 scalar_field_type _w;
777 template<
class M,
class X,
class Y,
int l=1>
778 class DUNE_DEPRECATED_MSG("Use SeqILU instead of SeqILUn!") SeqILUn :
public Preconditioner<X,Y> {
781 typedef typename std::remove_const<M>::type matrix_type;
783 typedef X domain_type;
785 typedef Y range_type;
787 typedef typename X::field_type field_type;
789 typedef Simd::Scalar<field_type> scalar_field_type;
798 SeqILUn (
const M& A,
int n, scalar_field_type w)
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");
813 _w = configuration.get<scalar_field_type>(
"relaxation");
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)
844 virtual void post (X& x)
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 );
Y range_type
range type of the preconditioner
Definition: preconditioners.hh:983
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:401
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:343
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:349
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:635
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
X domain_type
domain type of the preconditioner
Definition: preconditioners.hh:981
void bilu_backsolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:95
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:987
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:141
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
Sequential SSOR preconditioner.
Definition: preconditioners.hh:138
Col col
Definition: matrixmatrix.hh:349
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:509
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:944
virtual void apply(X &v, const Y &d)
Apply the preconditoner.
Definition: preconditioners.hh:607
Sequential ILU0 preconditioner.
Definition: preconditioners.hh:875
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:190
SeqJac(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:429
std::vector< block_type, typename matrix_type::allocator_type > inv_
Definition: preconditioners.hh:647
std::unique_ptr< matrix_type > ILU_
The ILU(n) decomposition of the matrix. As storage a BCRSMatrix is used.
Definition: preconditioners.hh:642
std::remove_const_t< M > matrix_type
type of matrix the preconditioner is for
Definition: preconditioners.hh:979
const bool wNotIdentity_
true if w != 1.0
Definition: preconditioners.hh:652
range_type::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:81
O InverseOperator
type of the wrapped inverse operator
Definition: preconditioners.hh:85
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:640
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:512
Some handy generic functions for ISTL matrices.
The sequential jacobian preconditioner.
Definition: preconditioners.hh:390
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:301
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:215
Turns an InverseOperator into a Preconditioner.
Definition: preconditioners.hh:72
Abstract base class for all solvers.
Definition: solver.hh:97
X::field_type field_type
field type of the preconditioner
Definition: preconditioners.hh:985
Richardson(const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:906
SeqJac(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:410
matrix_type ::block_type block_type
block type of matrix
Definition: preconditioners.hh:505
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:878
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:442
CRS upper_
Definition: preconditioners.hh:646
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:1073
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:258
const scalar_field_type w_
The relaxation factor to use.
Definition: preconditioners.hh:650
SeqSSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:177
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:312
Richardson(scalar_field_type w=1.0)
Constructor.
Definition: preconditioners.hh:891
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:256
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: preconditioners.hh:927
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:938
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:515
SeqSSOR(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:158
O::domain_type domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:77
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:596
virtual void post(domain_type &)
Clean up.
Definition: preconditioners.hh:108
SeqILU(const M &A, int n, scalar_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:560
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
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:52
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:503
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:112
InverseOperator2Preconditioner(InverseOperator &inverse_operator)
Construct the preconditioner from the solver.
Definition: preconditioners.hh:91
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:145
Define general, extensible interface for inverse operators.
void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: preconditioners.hh:1053
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:399
O::range_type range_type
The range type of the preconditioner.
Definition: preconditioners.hh:79
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:471
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:453
Statistics about the application of an inverse operator.
Definition: solver.hh:45
SeqILDL(const matrix_type &A, scalar_field_type relax=scalar_field_type(1))
constructor
Definition: preconditioners.hh:1013
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:916
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:202
SeqILU(const M &A, scalar_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:527
void bilu_decomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:169
ILU::CRS< block_type, typename M::allocator_type > CRS
type of ILU storage
Definition: preconditioners.hh:518
derive error class from the base class in common
Definition: istlexception.hh:16
void apply(X &v, const Y &d)
Apply the preconditioner in a special direction.
Definition: preconditioners.hh:326
SeqSOR(const M &A, int n, scalar_field_type w)
Constructor.
Definition: preconditioners.hh:269
Definition: allocator.hh:7
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:143
void post(X &x) override
Clean up.
Definition: preconditioners.hh:1067
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
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:507
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:652
void bilu0_decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:34
Incomplete LDL decomposition.
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:83
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:395
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:393
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:882
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:260
Sequential ILU preconditioner.
Definition: preconditioners.hh:500
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
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:254
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition: ildl.hh:86
Matrix & mat
Definition: matrixmatrix.hh:345
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:221
SeqILU(const M &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:546
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
SeqSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:288
compile-time parameter for block recursion depth
Definition: gsetc.hh:43
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:880
Sequential SOR preconditioner.
Definition: preconditioners.hh:249
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:465
SeqILDL(const matrix_type &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:1001
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:629
Category
Definition: solvercategory.hh:21
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
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:884
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:628
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:252
virtual void pre(domain_type &, range_type &)
Prepare the preconditioner.
Definition: preconditioners.hh:98
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:147
void apply(X &v, const Y &d) override
Apply one step of the preconditioner to the system A(v)=d.
Definition: preconditioners.hh:1060
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:149
CRS lower_
The ILU(n) decomposition of the matrix. As storage a CRS structure is used.
Definition: preconditioners.hh:645
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:397
sequential ILDL preconditioner
Definition: preconditioners.hh:971
void bildl_backsolve(const Matrix &A, X &v, const Y &d, bool isLowerTriangular=false)
Definition: ildl.hh:147