5#ifndef DUNE_ISTL_PRECONDITIONERS_HH
6#define DUNE_ISTL_PRECONDITIONERS_HH
15#include <dune/common/simd/simd.hh>
16#include <dune/common/parametertree.hh>
72 template<
class O,
int c = -1>
74 public Preconditioner<typename O::domain_type, typename O::range_type>
95 : inverse_operator_(inverse_operator)
98 DUNE_THROW(InvalidStateException,
"User-supplied solver category does not match that of the given inverse operator");
108 inverse_operator_.apply(v, copy, res);
140 template<
class M,
class X,
class Y,
int l=1>
164 : _A_(A), _n(n), _w(w)
183 :
SeqSSOR(A->getmat(), configuration)
199 SeqSSOR (
const M& A,
const ParameterTree& configuration)
208 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
216 virtual void apply (X& v,
const Y& d)
218 for (
int i=0; i<_n; i++) {
229 virtual void post ([[maybe_unused]] X& x)
260 template<
class M,
class X,
class Y,
int l=1>
284 : _A_(A), _n(n), _w(w)
303 :
SeqSOR(A->getmat(), configuration)
319 SeqSOR (
const M& A,
const ParameterTree& configuration)
328 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
336 virtual void apply (X& v,
const Y& d)
338 this->
template apply<true>(v,d);
349 template<
bool forward>
353 for (
int i=0; i<_n; i++) {
357 for (
int i=0; i<_n; i++) {
367 virtual void post ([[maybe_unused]] X& x)
397 template<
class M,
class X,
class Y,
int l=1>
411 template<
class M,
class X,
class Y,
int l=1>
435 : _A_(A), _n(n), _w(w)
454 :
SeqJac(A->getmat(), configuration)
470 SeqJac (
const M& A,
const ParameterTree& configuration)
479 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
487 virtual void apply (X& v,
const Y& d)
489 for (
int i=0; i<_n; i++) {
499 virtual void post ([[maybe_unused]] X& x)
531 template<
class M,
class X,
class Y,
int l=1>
562 :
SeqILU( A, 0, w, resort )
581 :
SeqILU(A->getmat(), configuration)
598 SeqILU(
const M& A,
const ParameterTree& config)
601 config.
get(
"resort", false))
618 wNotIdentity_([w]{
using std::abs;
return abs(w -
real_field_type(1)) > 1e-15;}() )
623 ILU_.reset(
new matrix_type( A ) );
630 ILU_.reset(
new matrix_type( A.
N(), A.
M(), matrix_type::row_wise) );
648 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
656 virtual void apply (X& v,
const Y& d)
678 virtual void post ([[maybe_unused]] X& x)
689 std::unique_ptr< matrix_type >
ILU_;
694 std::vector< block_type, typename matrix_type::allocator_type >
inv_;
712 template<
class X,
class Y>
755 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
763 virtual void apply (X& v,
const Y& d)
774 virtual void post ([[maybe_unused]] X& x)
788 using D =
typename Dune::TypeListElement<1,
decltype(tl)>::type;
789 using R =
typename Dune::TypeListElement<2,
decltype(tl)>::type;
790 return std::make_shared<Richardson<D,R>>(config);
804 template<
class M,
class X,
class Y >
838 :
SeqILDL(A->getmat(), configuration)
866 : decomposition_( A.N(), A.M(),
matrix_type::random ),
870 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
872 const auto &A_i = *i;
873 const auto ij = A_i.find( i.index() );
874 if( ij != A_i.end() )
875 decomposition_.setrowsize( i.index(), ij.offset()+1 );
877 DUNE_THROW(
ISTLError,
"diagonal entry missing" );
879 decomposition_.endrowsizes();
882 for(
auto i = A.begin(), iend = A.end(); i != iend; ++i )
884 const auto &A_i = *i;
885 for(
auto ij = A_i.begin(); ij.index() < i.index() ; ++ij )
886 decomposition_.addindex( i.index(), ij.index() );
887 decomposition_.addindex( i.index(), i.index() );
889 decomposition_.endindices();
893 for(
auto row = decomposition_.begin(), rowend = decomposition_.end(); row != rowend; ++row, ++i )
895 auto ij = i->begin();
896 for(
auto col = row->begin(), colend = row->end();
col != colend; ++
col, ++ij )
905 void pre ([[maybe_unused]] X &x, [[maybe_unused]] Y &b)
override
909 void apply ( X &v,
const Y &d )
override
916 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:16
Col col
Definition matrixmatrix.hh:351
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition gsetc.hh:646
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition gsetc.hh:658
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition gsetc.hh:634
Definition allocator.hh:11
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition ildl.hh:88
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:293
void bildl_backsolve(const Matrix &A, X &v, const Y &d, bool isLowerTriangular=false)
Definition ildl.hh:149
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:307
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition ilu.hh:94
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition ilu.hh:33
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition ilu.hh:167
compile-time parameter for block recursion depth
Definition gsetc.hh:45
a simple compressed row storage matrix class
Definition ilu.hh:259
derive error class from the base class in common
Definition istlexception.hh:19
size_type M() const
Return the number of columns.
Definition matrix.hh:700
size_type N() const
Return the number of rows.
Definition matrix.hh:695
static void check(const Matrix &mat)
Check whether the a matrix has diagonal values on blocklevel recursion levels.
Definition matrixutils.hh:53
A linear operator exporting itself in matrix form.
Definition operators.hh:109
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:32
Turns an InverseOperator into a Preconditioner.
Definition preconditioners.hh:75
O::range_type range_type
The range type of the preconditioner.
Definition preconditioners.hh:80
O::domain_type domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:78
virtual void post(domain_type &)
Clean up.
Definition preconditioners.hh:111
range_type::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:82
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:86
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:84
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:115
virtual void pre(domain_type &, range_type &)
Prepare the preconditioner.
Definition preconditioners.hh:101
InverseOperator2Preconditioner(InverseOperator &inverse_operator)
Construct the preconditioner from the solver.
Definition preconditioners.hh:94
O InverseOperator
type of the wrapped inverse operator
Definition preconditioners.hh:88
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:104
Sequential SSOR preconditioner.
Definition preconditioners.hh:141
virtual void post(X &x)
Clean up.
Definition preconditioners.hh:229
SeqSSOR(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:182
SeqSSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:199
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:233
X::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:150
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:152
X domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:146
M matrix_type
The matrix type the preconditioner is for.
Definition preconditioners.hh:144
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition preconditioners.hh:216
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition preconditioners.hh:208
Y range_type
The range type of the preconditioner.
Definition preconditioners.hh:148
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:154
SeqSSOR(const M &A, int n, real_field_type w)
Constructor.
Definition preconditioners.hh:163
Sequential SOR preconditioner.
Definition preconditioners.hh:261
SeqSOR(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:302
M matrix_type
The matrix type the preconditioner is for.
Definition preconditioners.hh:264
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:274
void apply(X &v, const Y &d)
Apply the preconditioner in a special direction.
Definition preconditioners.hh:350
X domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:266
virtual void post(X &x)
Clean up.
Definition preconditioners.hh:367
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition preconditioners.hh:328
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:272
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:371
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition preconditioners.hh:336
Y range_type
The range type of the preconditioner.
Definition preconditioners.hh:268
SeqSOR(const M &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:319
X::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:270
SeqSOR(const M &A, int n, real_field_type w)
Constructor.
Definition preconditioners.hh:283
The sequential jacobian preconditioner.
Definition preconditioners.hh:412
virtual void post(X &x)
Clean up.
Definition preconditioners.hh:499
SeqJac(const M &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:470
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition preconditioners.hh:487
M matrix_type
The matrix type the preconditioner is for.
Definition preconditioners.hh:415
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:423
SeqJac(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:453
X::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:421
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition preconditioners.hh:479
X domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:417
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:425
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:503
SeqJac(const M &A, int n, real_field_type w)
Constructor.
Definition preconditioners.hh:434
Y range_type
The range type of the preconditioner.
Definition preconditioners.hh:419
Sequential ILU preconditioner.
Definition preconditioners.hh:532
virtual void post(X &x)
Clean up.
Definition preconditioners.hh:678
SeqILU(const M &A, int n, real_field_type w, const bool resort=false)
Constructor.
Definition preconditioners.hh:612
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition preconditioners.hh:648
virtual void apply(X &v, const Y &d)
Apply the preconditioner.
Definition preconditioners.hh:656
ILU::CRS< block_type, typename M::allocator_type > CRS
type of ILU storage
Definition preconditioners.hh:552
Y range_type
The range type of the preconditioner.
Definition preconditioners.hh:541
CRS lower_
The ILU(n) decomposition of the matrix. As storage a CRS structure is used.
Definition preconditioners.hh:692
const bool wNotIdentity_
true if w != 1.0
Definition preconditioners.hh:699
SeqILU(const M &A, const ParameterTree &config)
Constructor.
Definition preconditioners.hh:598
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition preconditioners.hh:535
matrix_type::block_type block_type
block type of matrix
Definition preconditioners.hh:537
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:549
X::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:544
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:682
SeqILU(const M &A, real_field_type w, const bool resort=false)
Constructor.
Definition preconditioners.hh:561
const real_field_type w_
The relaxation factor to use.
Definition preconditioners.hh:697
SeqILU(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:580
X domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:539
std::vector< block_type, typename matrix_type::allocator_type > inv_
Definition preconditioners.hh:694
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:547
std::unique_ptr< matrix_type > ILU_
The ILU(n) decomposition of the matrix. As storage a BCRSMatrix is used.
Definition preconditioners.hh:689
CRS upper_
Definition preconditioners.hh:693
Richardson preconditioner.
Definition preconditioners.hh:713
X::field_type field_type
The field type of the preconditioner.
Definition preconditioners.hh:720
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:778
Y range_type
The range type of the preconditioner.
Definition preconditioners.hh:718
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition preconditioners.hh:755
Richardson(real_field_type w=1.0)
Constructor.
Definition preconditioners.hh:731
virtual void post(X &x)
Clean up.
Definition preconditioners.hh:774
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:724
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:722
Richardson(const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:746
X domain_type
The domain type of the preconditioner.
Definition preconditioners.hh:716
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition preconditioners.hh:763
sequential ILDL preconditioner
Definition preconditioners.hh:807
SeqILDL(const matrix_type &A, const ParameterTree &config)
Constructor.
Definition preconditioners.hh:853
SeqILDL(const matrix_type &A, real_field_type relax=real_field_type(1))
constructor
Definition preconditioners.hh:865
X domain_type
domain type of the preconditioner
Definition preconditioners.hh:815
void post(X &x) override
Clean up.
Definition preconditioners.hh:916
Y range_type
range type of the preconditioner
Definition preconditioners.hh:817
std::remove_const_t< M > matrix_type
type of matrix the preconditioner is for
Definition preconditioners.hh:813
void apply(X &v, const Y &d) override
Apply one step of the preconditioner to the system A(v)=d.
Definition preconditioners.hh:909
FieldTraits< scalar_field_type >::real_type real_field_type
real scalar type underlying the field_type
Definition preconditioners.hh:823
SeqILDL(const std::shared_ptr< const AssembledLinearOperator< M, X, Y > > &A, const ParameterTree &configuration)
Constructor.
Definition preconditioners.hh:837
void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition preconditioners.hh:905
Simd::Scalar< field_type > scalar_field_type
scalar type underlying the field_type
Definition preconditioners.hh:821
X::field_type field_type
field type of the preconditioner
Definition preconditioners.hh:819
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition preconditioners.hh:920
Statistics about the application of an inverse operator.
Definition solver.hh:48
Abstract base class for all solvers.
Definition solver.hh:99
Category
Definition solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25
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:34