3#ifndef DUNE_ISTL_PRECONDITIONERS_HH
4#define DUNE_ISTL_PRECONDITIONERS_HH
13#include <dune/common/simd/simd.hh>
14#include <dune/common/parametertree.hh>
70 template<
class O,
int c = -1>
72 public Preconditioner<typename O::domain_type, typename O::range_type>
93 : inverse_operator_(inverse_operator)
96 DUNE_THROW(InvalidStateException,
"User-supplied solver category does not match that of the given inverse operator");
106 inverse_operator_.apply(v, copy, res);
138 template<
class M,
class X,
class Y,
int l=1>
162 : _A_(A), _n(n), _w(w)
181 :
SeqSSOR(A->getmat(), configuration)
197 SeqSSOR (
const M& A,
const ParameterTree& configuration)
206 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
214 virtual void apply (X& v,
const Y& d)
216 for (
int i=0; i<_n; i++) {
227 virtual void post ([[maybe_unused]] X& x)
258 template<
class M,
class X,
class Y,
int l=1>
282 : _A_(A), _n(n), _w(w)
301 :
SeqSOR(A->getmat(), configuration)
317 SeqSOR (
const M& A,
const ParameterTree& configuration)
326 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
334 virtual void apply (X& v,
const Y& d)
336 this->
template apply<true>(v,d);
347 template<
bool forward>
351 for (
int i=0; i<_n; i++) {
355 for (
int i=0; i<_n; i++) {
365 virtual void post ([[maybe_unused]] X& x)
395 template<
class M,
class X,
class Y,
int l=1>
409 template<
class M,
class X,
class Y,
int l=1>
433 : _A_(A), _n(n), _w(w)
452 :
SeqJac(A->getmat(), configuration)
468 SeqJac (
const M& A,
const ParameterTree& configuration)
477 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
485 virtual void apply (X& v,
const Y& d)
487 for (
int i=0; i<_n; i++) {
497 virtual void post ([[maybe_unused]] X& x)
529 template<
class M,
class X,
class Y,
int l=1>
560 :
SeqILU( A, 0, w, resort )
579 :
SeqILU(A->getmat(), configuration)
596 SeqILU(
const M& A,
const ParameterTree& config)
599 config.
get(
"resort", false))
616 wNotIdentity_([w]{
using std::abs;
return abs(w -
real_field_type(1)) > 1e-15;}() )
621 ILU_.reset(
new matrix_type( A ) );
628 ILU_.reset(
new matrix_type( A.
N(), A.
M(), matrix_type::row_wise) );
646 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
654 virtual void apply (X& v,
const Y& d)
676 virtual void post ([[maybe_unused]] X& x)
687 std::unique_ptr< matrix_type >
ILU_;
692 std::vector< block_type, typename matrix_type::allocator_type >
inv_;
710 template<
class X,
class Y>
753 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
761 virtual void apply (X& v,
const Y& d)
772 virtual void post ([[maybe_unused]] X& x)
786 using D =
typename Dune::TypeListElement<1,
decltype(tl)>::type;
787 using R =
typename Dune::TypeListElement<2,
decltype(tl)>::type;
788 return std::make_shared<Richardson<D,R>>(config);
802 template<
class M,
class X,
class Y >
836 :
SeqILDL(A->getmat(), configuration)
864 : decomposition_( A.N(), A.M(),
matrix_type::random ),
868 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
870 const auto &A_i = *i;
871 const auto ij = A_i.find( i.index() );
872 if( ij != A_i.end() )
873 decomposition_.setrowsize( i.index(), ij.offset()+1 );
875 DUNE_THROW(
ISTLError,
"diagonal entry missing" );
877 decomposition_.endrowsizes();
880 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
882 const auto &A_i = *i;
883 for(
auto ij = A_i.begin(); ij.index() < i.index() ; ++ij )
884 decomposition_.addindex( i.index(), ij.index() );
885 decomposition_.addindex( i.index(), i.index() );
887 decomposition_.endindices();
891 for(
auto row = decomposition_.begin(), rowend = decomposition_.end(); row != rowend; ++row, ++i )
893 auto ij = i->begin();
894 for(
auto col = row->begin(), colend = row->end();
col != colend; ++
col, ++ij )
903 void pre ([[maybe_unused]] X &x, [[maybe_unused]] Y &b)
override
907 void apply ( X &v,
const Y &d )
override
914 void post ([[maybe_unused]] X &x)
override
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Incomplete LDL decomposition.
The incomplete LU factorization kernels.
Some handy generic functions for ISTL matrices.
Define general, extensible interface for inverse operators.
#define DUNE_REGISTER_PRECONDITIONER(name,...)
Definition: solverregistry.hh:14
Col col
Definition: matrixmatrix.hh:349
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:644
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:656
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:632
Definition: allocator.hh:9
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:291
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:306
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:93
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:32
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:166
compile-time parameter for block recursion depth
Definition: gsetc.hh:43
a simple compressed row storage matrix class
Definition: ilu.hh:258
derive error class from the base class in common
Definition: istlexception.hh:17
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition: matrixutils.hh:51
A linear operator exporting itself in matrix form.
Definition: operators.hh:107
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Turns an InverseOperator into a Preconditioner.
Definition: preconditioners.hh:73
O::range_type range_type
The range type of the preconditioner.
Definition: preconditioners.hh:78
O::domain_type domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:76
virtual void post(domain_type &)
Clean up.
Definition: preconditioners.hh:109
range_type::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:80
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:84
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:82
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:113
virtual void pre(domain_type &, range_type &)
Prepare the preconditioner.
Definition: preconditioners.hh:99
InverseOperator2Preconditioner(InverseOperator &inverse_operator)
Construct the preconditioner from the solver.
Definition: preconditioners.hh:92
O InverseOperator
type of the wrapped inverse operator
Definition: preconditioners.hh:86
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:102
Sequential SSOR preconditioner.
Definition: preconditioners.hh:139
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:227
SeqSSOR(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:180
SeqSSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:197
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:231
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:148
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:150
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:144
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:142
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:214
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:206
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:146
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:152
SeqSSOR(const M &A, int n, real_field_type w)
Constructor.
Definition: preconditioners.hh:161
Sequential SOR preconditioner.
Definition: preconditioners.hh:259
SeqSOR(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:300
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:262
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:272
void apply(X &v, const Y &d)
Apply the preconditioner in a special direction.
Definition: preconditioners.hh:348
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:264
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:365
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:326
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:270
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:369
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:334
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:266
SeqSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:317
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:268
SeqSOR(const M &A, int n, real_field_type w)
Constructor.
Definition: preconditioners.hh:281
The sequential jacobian preconditioner.
Definition: preconditioners.hh:410
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:497
SeqJac(const M &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:468
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:485
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:413
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:421
SeqJac(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:451
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:419
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:477
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:415
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:423
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:501
SeqJac(const M &A, int n, real_field_type w)
Constructor.
Definition: preconditioners.hh:432
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:417
Sequential ILU preconditioner.
Definition: preconditioners.hh:530
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:676
SeqILU(const M &A, int n, real_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:610
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:646
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition: preconditioners.hh:654
ILU::CRS< block_type, typename M::allocator_type > CRS
type of ILU storage
Definition: preconditioners.hh:550
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:539
CRS lower_
The ILU(n) decomposition of the matrix. As storage a CRS structure is used.
Definition: preconditioners.hh:690
const bool wNotIdentity_
true if w != 1.0
Definition: preconditioners.hh:697
SeqILU(const M &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:596
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:533
matrix_type::block_type block_type
block type of matrix
Definition: preconditioners.hh:535
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:547
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:542
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:680
SeqILU(const M &A, real_field_type w, const bool resort=false)
Constructor.
Definition: preconditioners.hh:559
const real_field_type w_
The relaxation factor to use.
Definition: preconditioners.hh:695
SeqILU(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:578
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:537
std::vector< block_type, typename matrix_type::allocator_type > inv_
Definition: preconditioners.hh:692
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:545
std::unique_ptr< matrix_type > ILU_
The ILU(n) decomposition of the matrix. As storage a BCRSMatrix is used.
Definition: preconditioners.hh:687
CRS upper_
Definition: preconditioners.hh:691
Richardson preconditioner.
Definition: preconditioners.hh:711
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:718
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:776
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:716
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:753
Richardson(real_field_type w=1.0)
Constructor.
Definition: preconditioners.hh:729
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:772
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:722
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:720
Richardson(const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:744
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:714
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: preconditioners.hh:761
sequential ILDL preconditioner
Definition: preconditioners.hh:805
SeqILDL(const matrix_type &A, const ParameterTree &config)
Constructor.
Definition: preconditioners.hh:851
SeqILDL(const matrix_type &A, real_field_type relax=real_field_type(1))
constructor
Definition: preconditioners.hh:863
X domain_type
domain type of the preconditioner
Definition: preconditioners.hh:813
void post(X &x) override
Clean up.
Definition: preconditioners.hh:914
Y range_type
range type of the preconditioner
Definition: preconditioners.hh:815
std::remove_const_t< M > matrix_type
type of matrix the preconditioner is for
Definition: preconditioners.hh:811
void apply(X &v, const Y &d) override
Apply one step of the preconditioner to the system A(v)=d.
Definition: preconditioners.hh:907
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition: preconditioners.hh:821
SeqILDL(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition: preconditioners.hh:835
void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: preconditioners.hh:903
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition: preconditioners.hh:819
X::field_type field_type
field type of the preconditioner
Definition: preconditioners.hh:817
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:918
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