Preconditioners and Relaxation Operators

Collaboration diagram for Preconditioners and Relaxation Operators:

Topics
	Exceptions and assertions

Namespaces
namespace	internal

Classes
class	PreconditionLU< number >
class	PreconditionIdentity
class	PreconditionRichardson
class	PreconditionUseMatrix< MatrixType, VectorType >
class	PreconditionRelaxation< MatrixType, PreconditionerType >
class	PreconditionJacobi< MatrixType >
class	PreconditionSOR< MatrixType >
class	PreconditionSSOR< MatrixType >
class	PreconditionPSOR< MatrixType >
class	PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >
class	PreconditionBlock< MatrixType, inverse_type >
class	PreconditionBlockJacobi< MatrixType, inverse_type >
class	PreconditionBlockSOR< MatrixType, inverse_type >
class	PreconditionBlockSSOR< MatrixType, inverse_type >
class	PreconditionSelector< MatrixType, VectorType >
class	RelaxationBlock< MatrixType, InverseNumberType, VectorType >
class	RelaxationBlockJacobi< MatrixType, InverseNumberType, VectorType >
class	RelaxationBlockSOR< MatrixType, InverseNumberType, VectorType >
class	RelaxationBlockSSOR< MatrixType, InverseNumberType, VectorType >
class	AdditionalData
class	SparseDirectUMFPACK
class	SparseDirectMUMPS
class	SparseILU< number >
class	SparseMIC< number >
class	SparseVanka< number >
class	SparseBlockVanka< number >
class	TrilinosWrappers::PreconditionBase
class	TrilinosWrappers::PreconditionJacobi
class	TrilinosWrappers::PreconditionSSOR
class	TrilinosWrappers::PreconditionSOR
class	TrilinosWrappers::PreconditionBlockJacobi
class	TrilinosWrappers::PreconditionBlockSSOR
class	TrilinosWrappers::PreconditionBlockSOR
class	TrilinosWrappers::PreconditionIC
class	TrilinosWrappers::PreconditionILU
class	TrilinosWrappers::PreconditionILUT
class	TrilinosWrappers::PreconditionBlockwiseDirect
class	TrilinosWrappers::PreconditionChebyshev
class	TrilinosWrappers::PreconditionAMG
class	TrilinosWrappers::PreconditionAMGMueLu
class	TrilinosWrappers::PreconditionIdentity
class	LinearAlgebra::TpetraWrappers::PreconditionBase< Number, MemorySpace >
class	LinearAlgebra::TpetraWrappers::PreconditionIdentity< Number, MemorySpace >
	Wrapper class for the IdentitySolver preconditioner of Ifpack2. More...
class	LinearAlgebra::TpetraWrappers::PreconditionIfpackBase< Number, MemorySpace >
	The base class for all Ifpack2 preconditioners which are handled through its Factory. More...
class	LinearAlgebra::TpetraWrappers::PreconditionIfpack< Number, MemorySpace >
	Wrapper to create custom Ifpack2 preconditioners. More...
class	LinearAlgebra::TpetraWrappers::PreconditionJacobi< Number, MemorySpace >
	The classical Jacobi preconditioner within Ifpack2. More...
class	LinearAlgebra::TpetraWrappers::PreconditionL1Jacobi< Number, MemorySpace >
	l1 variant of the Jacobi preconditioner. More...
class	LinearAlgebra::TpetraWrappers::PreconditionL1GaussSeidel< Number, MemorySpace >
	l1 variant of the Gauss-Seidel preconditioner. More...
class	LinearAlgebra::TpetraWrappers::PreconditionSOR< Number, MemorySpace >
	The class for the SOR preconditioner within Ifpack2. More...
class	LinearAlgebra::TpetraWrappers::PreconditionSSOR< Number, MemorySpace >
class	LinearAlgebra::TpetraWrappers::PreconditionChebyshev< Number, MemorySpace >
class	LinearAlgebra::TpetraWrappers::PreconditionILU< Number, MemorySpace >
	The ILU/ILU(K)/RILU(K) preconditioner. More...
class	LinearAlgebra::TpetraWrappers::PreconditionILUT< Number, MemorySpace >
class	LinearAlgebra::TpetraWrappers::PreconditionBlockJacobi< Number, MemorySpace >
class	LinearAlgebra::TpetraWrappers::PreconditionBlockSOR< Number, MemorySpace >
	The class for the Block SOR preconditioner within Ifpack2. More...
class	LinearAlgebra::TpetraWrappers::PreconditionBlockSSOR< Number, MemorySpace >
	The class for the Block SSOR preconditioner within Ifpack2. More...

Typedefs
using	size_type = typename SparseMatrix<number>::size_type

Functions
	SparseLUDecomposition ()
virtual	~SparseLUDecomposition () override=0
virtual void	clear () override
template<typename somenumber>
void	initialize (const SparseMatrix< somenumber > &matrix, const AdditionalData parameters)
bool	empty () const
size_type	m () const
size_type	n () const
template<class OutVector, class InVector>
void	vmult_add (OutVector &dst, const InVector &src) const
template<class OutVector, class InVector>
void	Tvmult_add (OutVector &dst, const InVector &src) const
virtual std::size_t	memory_consumption () const
template<typename somenumber>
void	copy_from (const SparseMatrix< somenumber > &matrix)
virtual void	strengthen_diagonal_impl ()
virtual number	get_strengthen_diagonal (const number rowsum, const size_type row) const
void	prebuild_lower_bound ()

Variables
double	strengthen_diagonal
std::vector< const size_type * >	prebuilt_lower_bound
SparsityPattern *	own_sparsity

Detailed Description

Preconditioners

Preconditioners are used to accelerate the iterative solution of linear systems. Typical preconditioners are Jacobi, Gauss-Seidel, or SSOR, but the library also supports more complex ones such as Vanka or incomplete LU decompositions (ILU). In addition, sparse direct solvers can be used as preconditioners when available.

Broadly speaking, preconditioners are operators, which are multiplied with a matrix to improve conditioning. The idea is, that the preconditioned system P^-1Ax = P^-1b is much easier to solve than the original system Ax = b. What this means exactly depends on the structure of the matrix and cannot be discussed here in generality. For symmetric, positive definite matrices A and P, it means that the spectral condition number (the quotient of greatest and smallest eigenvalue) of P^-1A is much smaller than the one of A.

At hand of the simplest example, Richardson iteration, implemented in SolverRichardson, the preconditioned iteration looks like

Accordingly, preconditioning amounts to applying a linear operator to the residual, and consequently, the action of the preconditioner P^-1 is implemented as vmult(). Templates in deal.II that require a preconditioner indicate the requirement with the PreconditionerType concept. In practice, one can usually treat any matrix-like object which defines vmult() and Tvmult() as a preconditioner. All preconditioner classes in this group implement this interface.

When used in Krylov space methods, it is up to the method, whether it simply replaces multiplications with A by those with P^-1A (for instance SolverBicgstab), or does more sophisticated things. SolverCG for instance uses P^-1 to define an inner product, which is the reason why it requires a symmetric, positive definite operator P.

Relaxation methods

Many preconditioners rely on an additive splitting A = P - N into two matrices. In this case, the iteration step of the Richardson method above can be simplified to

thus avoiding multiplication with A completely. We call operators mapping the previous iterate x^k to the next iterate in this way relaxation operators. Their generic interface is given by the RelaxationType concept. The classes with names starting with Relaxation in this group implement this interface, as well as the preconditioners PreconditionJacobi, PreconditionSOR, PreconditionBlockJacobi, PreconditionBlockSOR, and PreconditionBlockSSOR.

The interface

In this section, we discuss the interface preconditioners usually have to provide to work inside the deal.II library.

Initialization

In order to be able to be stored in containers, all preconditioners have a constructor with no arguments. Since this will typically produce a useless object, all preconditioners have a function

void initialize (...)

This function receives the matrix to be preconditioned as well as additional required parameters and sets up the internal structures of the preconditioner.

Relaxation methods

Some preconditioners, like SOR and Jacobi, were used as iterative solvers long before they were used as preconditioners. Thus, they satisfy both MatrixType and RelaxationType concepts.

Typedef Documentation

size_type

using size_type = typename SparseMatrix<number>::size_type

Declare type for container size.

Definition at line 122 of file sparse_decomposition.h.

Function Documentation

SparseLUDecomposition()

SparseLUDecomposition

(

)

Abstract base class for incomplete decompositions of a sparse matrix into sparse factors. This class can't be used by itself, but only as the base class of derived classes that actually implement particular decompositions such as SparseILU or SparseMIC.

The decomposition is stored as a sparse matrix which is why this class is derived from the SparseMatrix. Since it is not a matrix in the usual sense (the stored entries are not those of a matrix, but of the two factors of the original matrix), the derivation is protected rather than public.

Fill-in

Sparse decompositions are frequently used with additional fill-in, i.e., the sparsity structure of the decomposition is denser than that of the matrix to be decomposed. The initialize() function of this class allows this fill-in via the AdditionalData object as long as all entries present in the original matrix are present in the decomposition also, i.e. the sparsity pattern of the decomposition is a superset of the sparsity pattern in the original matrix.

Such fill-in can be accomplished by various ways, one of which is the copy-constructor of the SparsityPattern class that allows the addition of side-diagonals to a given sparsity structure.

Unified use of preconditioners

While objects of this class can not be used directly (this class is only a base class for others implementing actual decompositions), derived classes such as SparseILU and SparseMIC can be used in the usual form as preconditioners. For example, this works:

SparseILU<double> ilu;
ilu.initialize(matrix, SparseILU<double>::AdditionalData(...));
 
somesolver.solve (A, x, f, ilu);

Through the AdditionalData object it is possible to specify additional parameters of the LU decomposition.

1/ The matrix diagonal can be strengthened by adding strengthen_diagonal times the sum of the absolute row entries of each row to the respective diagonal entries. By default no strengthening is performed.

2/ By default, each initialize() function call creates its own sparsity. For that, it copies the sparsity of matrix and adds a specific number of extra off diagonal entries specified by extra_off_diagonals.

3/ By setting use_previous_sparsity=true the sparsity is not recreated but the sparsity of the previous initialize() call is reused (recycled). This might be useful when several linear problems on the same sparsity need to solved, as for example several Newton iteration steps on the same triangulation. The default is false.

4/ It is possible to give a user defined sparsity to use_this_sparsity. Then, no sparsity is created but *use_this_sparsity is used to store the decomposed matrix. For restrictions on the sparsity see section `Fill-in' above).

Particular implementations

It is enough to override the initialize() and vmult() methods to implement particular LU decompositions, like the true LU, or the Cholesky decomposition. Additionally, if that decomposition needs fine tuned diagonal strengthening on a per row basis, it may override the get_strengthen_diagonal() method. */ template <typename number> class SparseLUDecomposition : protected SparseMatrix<number>, public virtual EnableObserverPointer { protected: /** Constructor. Does nothing.

Call the initialize() function before using this object as preconditioner (vmult()).

Definition at line 1770 of file sparse_matrix.h.

~SparseLUDecomposition()

virtual ~SparseLUDecomposition ( )

overridepure virtual

Destruction. Mark the destructor pure to ensure that this class isn't used directly, but only its derived classes.

clear()

virtual void clear ( )

overridevirtual

Deletes all member variables. Leaves the class in the state that it had directly after calling the constructor

initialize()

template<typename somenumber>

void initialize	(	const SparseMatrix< somenumber > &	matrix,
		const AdditionalData	parameters )

This function needs to be called before an object of this class is used as preconditioner.

For more detail about possible parameters, see the class documentation and the documentation of the SparseLUDecomposition::AdditionalData class.

According to the parameters, this function creates a new SparsityPattern or keeps the previous sparsity or takes the sparsity given by the user to data. Then, this function performs the LU decomposition.

After this function is called the preconditioner is ready to be used (using the vmult function of derived classes).

empty()

bool empty

(

)

const

Return whether the object is empty. It calls the inherited SparseMatrix::empty() function.

m()

size_type m

(

)

const

Return the dimension of the codomain (or range) space. It calls the inherited SparseMatrix::m() function. Note that the matrix is of dimension .

n()

size_type n

(

)

const

Return the dimension of the domain space. It calls the inherited SparseMatrix::n() function. Note that the matrix is of dimension .

vmult_add()

template<class OutVector, class InVector>

void vmult_add	(	OutVector &	dst,
		const InVector &	src ) const

Adding Matrix-vector multiplication. Add M*src on dst with M being this matrix.

Source and destination must not be the same vector.

Tvmult_add()

template<class OutVector, class InVector>

void Tvmult_add	(	OutVector &	dst,
		const InVector &	src ) const

Adding Matrix-vector multiplication. Add M^T*src to dst with M being this matrix. This function does the same as vmult_add() but takes the transposed matrix.

Source and destination must not be the same vector.

memory_consumption()

virtual std::size_t memory_consumption ( ) const

virtual

Determine an estimate for the memory consumption (in bytes) of this object.

copy_from()

template<typename somenumber>

void copy_from ( const SparseMatrix< somenumber > & matrix )

protected

Copies the passed SparseMatrix onto this object. This object's sparsity pattern remains unchanged.

strengthen_diagonal_impl()

virtual void strengthen_diagonal_impl ( )

protectedvirtual

Performs the strengthening loop. For each row calculates the sum of absolute values of its elements, determines the strengthening factor (through get_strengthen_diagonal()) sf and multiplies the diagonal entry with sf+1.

get_strengthen_diagonal()

virtual number get_strengthen_diagonal ( const number rowsum, const size_type row ) const

protectedvirtual

In the decomposition phase, computes a strengthening factor for the diagonal entry in row row with sum of absolute values of its elements rowsum.

Note

The default implementation in SparseLUDecomposition returns strengthen_diagonal's value. This variable is set to a nonzero value in several of the derived classes.

prebuild_lower_bound()

void prebuild_lower_bound ( )

protected

Fills the prebuilt_lower_bound array.

Variable Documentation

strengthen_diagonal

double strengthen_diagonal

protected

The default strengthening value, returned by get_strengthen_diagonal().

Definition at line 311 of file sparse_decomposition.h.

prebuilt_lower_bound

std::vector<const size_type *> prebuilt_lower_bound

protected

For every row in the underlying SparsityPattern, this array contains a pointer to the row's first afterdiagonal entry. Becomes available after invocation of prebuild_lower_bound().

Definition at line 318 of file sparse_decomposition.h.

own_sparsity

SparsityPattern* own_sparsity

private

In general this pointer is zero except for the case that no SparsityPattern is given to this class. Then, a SparsityPattern is created and is passed down to the SparseMatrix base class.

Nevertheless, the SparseLUDecomposition needs to keep ownership of this sparsity. It keeps this pointer to it enabling it to delete this sparsity at destruction time.

Definition at line 336 of file sparse_decomposition.h.