Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Classes | Public Types | Public Member Functions | Static Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > Class Template Reference

#include <deal.II/lac/precondition.h>

Inheritance diagram for PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >:
Inheritance graph
[legend]

Classes

struct  AdditionalData
 

Public Types

using size_type = types::global_dof_index
 
using EigenvalueInformation = internal::EigenvalueInformation
 

Public Member Functions

 PreconditionChebyshev ()
 
void initialize (const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
 
void vmult (VectorType &dst, const VectorType &src) const
 
void Tvmult (VectorType &dst, const VectorType &src) const
 
void step (VectorType &dst, const VectorType &src) const
 
void Tstep (VectorType &dst, const VectorType &src) const
 
void clear ()
 
size_type m () const
 
size_type n () const
 
EigenvalueInformation estimate_eigenvalues (const VectorType &src) const
 
template<class Archive >
void serialize (Archive &ar, const unsigned int version)
 
Subscriptor functionality

Classes derived from Subscriptor provide a facility to subscribe to this object. This is mostly used by the SmartPointer class.

void subscribe (std::atomic< bool > *const validity, const std::string &identifier="") const
 
void unsubscribe (std::atomic< bool > *const validity, const std::string &identifier="") const
 
unsigned int n_subscriptions () const
 
template<typename StreamType >
void list_subscribers (StreamType &stream) const
 
void list_subscribers () const
 

Static Public Member Functions

static ::ExceptionBaseExcInUse (int arg1, std::string arg2, std::string arg3)
 
static ::ExceptionBaseExcNoSubscriber (std::string arg1, std::string arg2)
 

Private Types

using map_value_type = decltype(counter_map)::value_type
 
using map_iterator = decltype(counter_map)::iterator
 

Private Member Functions

void check_no_subscribers () const noexcept
 

Private Attributes

SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
 
VectorType solution_old
 
VectorType temp_vector1
 
VectorType temp_vector2
 
AdditionalData data
 
double theta
 
double delta
 
bool eigenvalues_are_initialized
 
Threads::Mutex mutex
 
std::atomic< unsigned intcounter
 
std::map< std::string, unsigned intcounter_map
 
std::vector< std::atomic< bool > * > validity_pointers
 
const std::type_info * object_info
 

Detailed Description

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
class PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >

Preconditioning with a Chebyshev polynomial for symmetric positive definite matrices. This preconditioner is based on an iteration of an inner preconditioner of type PreconditionerType with coefficients that are adapted to optimally cover an eigenvalue range between the largest eigenvalue \(\lambda_{\max{}}\) down to a given lower eigenvalue \(\lambda_{\min{}}\) specified by the optional parameter smoothing_range. The algorithm is based on the following three-term recurrence:

\[ x^{n+1} = x^{n} + \alpha^n_0 (x^{n} - x^{n-1}) + \alpha^n_1 P^{-1} (b-Ax^n) \quad\text{with}\quad \alpha^0_0 := 0,\; \alpha^0_1 := \frac{2\rho_0}{\lambda_{\max}-\lambda_{\min}}\; \alpha^n_0 := \rho_n \rho_{n-1},\;\text{and}\; \alpha^n_1 := \frac{4\rho_n}{\lambda_{\max}-\lambda_{\min}}, \]

where the parameter \(\rho_0\) is set to \(\rho_0 = \frac{\lambda_{\max{}}-\lambda_{\min{}}}{\lambda_{\max{}}+\lambda_{\min{}}}\) for the maximal eigenvalue \(\lambda_{\max{}}\) and updated via \(\rho_n = \left(2\frac{\lambda_{\max{}}+\lambda_{\min{}}} {\lambda_{\max{}}-\lambda_{\min{}}} - \rho_{n-1}\right)^{-1}\). The Chebyshev polynomial is constructed to strongly damp the eigenvalue range between \(\lambda_{\min{}}\) and \(\lambda_{\max{}}\) and is visualized e.g. in Utilities::LinearAlgebra::chebyshev_filter().

The typical use case for the preconditioner is a Jacobi preconditioner specified through DiagonalMatrix, which is also the default value for the preconditioner. Note that if the degree variable is set to one, the Chebyshev iteration corresponds to a Jacobi preconditioner (or the underlying preconditioner type) with relaxation parameter according to the specified smoothing range.

Besides the default choice of a pointwise Jacobi preconditioner, this class also allows for more advanced types of preconditioners, for example iterating block-Jacobi preconditioners in DG methods.

Apart from the inner preconditioner object, this iteration does not need access to matrix entries, which makes it an ideal ingredient for matrix-free computations. In that context, this class can be used as a multigrid smoother that is trivially parallel (assuming that matrix-vector products are parallel and the inner preconditioner is parallel). Its use is demonstrated in the step-37 and step-59 tutorial programs.

Estimation of the eigenvalues

The Chebyshev method relies on an estimate of the eigenvalues of the matrix which are computed during the first invocation of vmult(). This class offers several algorithms to this end, see PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm. The default algorithm invokes the Lanczos method via the SolverCG class, which requires symmetry and positive definiteness of the (preconditioned) matrix system are required. Also note that deal.II needs to be configured with LAPACK support to use this option. The eigenvalue algorithm can be controlled by PreconditionChebyshev::AdditionalData::eig_cg_n_iterations specifying how many iterations should be performed. For all algorithms, the iterative process is started from an initial vector that depends on the vector type. For the classes Vector or LinearAlgebra::distributed::Vector, which have fast element access, it is a vector with entries (-5.5, -4.5, -3.5, -2.5, ..., 3.5, 4.5, 5.5) with appropriate epilogue and adjusted such that its mean is always zero, which works well for the Laplacian. This setup is stable in parallel in the sense that for a different number of processors but the same ordering of unknowns, the same initial vector and thus eigenvalue distribution will be computed, apart from roundoff errors. For other vector types, the initial vector contains all ones, scaled by the length of the vector, except for the very first entry that is zero, triggering high-frequency content again.

The computation of eigenvalues happens the first time one of the vmult(), Tvmult(), step() or Tstep() functions is called or when estimate_eigenvalues() is called directly. In the latter case, it is necessary to provide a temporary vector of the same layout as the source and destination vectors used during application of the preconditioner.

The estimates for minimum and maximum eigenvalue are taken from the underlying solver or eigenvalue algorithm in the given number of iterations, even if the solver did not converge in the requested number of iterations. Finally, the maximum eigenvalue is multiplied by a safety factor of 1.2.

Due to the cost of the eigenvalue estimate, this class is most appropriate if it is applied repeatedly, e.g. in a smoother for a geometric multigrid solver, that can in turn be used to solve several linear systems.

Bypassing the eigenvalue computation

In some contexts, the automatic eigenvalue computation of this class may result in a bad quality, e.g. when the polynomial basis or numbering of unknowns is such that the initial vector described above is a bad choice. It is possible to bypass the automatic eigenvalue computation by setting AdditionalData::eig_cg_n_iterations to zero, and provide the variable AdditionalData::max_eigenvalue instead. The minimal eigenvalue is implicitly specified via max_eigenvalue/smoothing_range.

Using the PreconditionChebyshev as a solver

If the range [max_eigenvalue/smoothing_range, max_eigenvalue] contains all eigenvalues of the preconditioned matrix system and the degree (i.e., number of iterations) is high enough, this class can also be used as a direct solver. For an error estimation of the Chebyshev iteration that can be used to determine the number of iteration, see Varga (2009).

In order to use Chebyshev as a solver, set the degree to numbers::invalid_unsigned_int to force the automatic computation of the number of iterations needed to reach a given target tolerance. Note that this currently only works for symmetric positive definite matrices with the eigenvalue algorithm set to the conjugate gradient algorithm. In this case, the target tolerance is read from the variable PreconditionChebyshev::AdditionalData::smoothing_range (it needs to be a number less than one to force any iterations obviously).

For details on the algorithm, see section 5.1 of

@book{Varga2009,
Title = {Matrix iterative analysis},
Author = {Varga, R. S.},
Publisher = {Springer},
Address = {Berlin},
Edition = {2nd},
Year = {2009},
}

Requirements on the templated classes

The class MatrixType must be derived from Subscriptor because a SmartPointer to MatrixType is held in the class. In particular, this means that the matrix object needs to persist during the lifetime of PreconditionChebyshev. The preconditioner is held in a shared_ptr that is copied into the AdditionalData member variable of the class, so the variable used for initialization can safely be discarded after calling initialize(). Both the matrix and the preconditioner need to provide vmult() functions for the matrix-vector product and m() functions for accessing the number of rows in the (square) matrix. Furthermore, the matrix must provide el(i,i) methods for accessing the matrix diagonal in case the preconditioner type is DiagonalMatrix. Even though it is highly recommended to pass the inverse diagonal entries inside a separate preconditioner object for implementing the Jacobi method (which is the only possible way to operate this class when computing in parallel with MPI because there is no knowledge about the locally stored range of entries that would be needed from the matrix alone), there is a backward compatibility function that can extract the diagonal in case of a serial computation.

Optimized operations with specific MatrixType argument

This class enables to embed the vector updates into the matrix-vector product in case the MatrixType supports this. To this end, the VectorType needs to be of type LinearAlgebra::distributed::Vector, the PreconditionerType needs to be DiagonalMatrix, and the class MatrixType needs to provide a function with the signature

void MatrixType::vmult(
VectorType &,
const VectorType &,
const std::function<void(const unsigned int, const unsigned int)> &,
const std::function<void(const unsigned int, const unsigned int)> &) const

where the two given functions run before and after the matrix-vector product, respectively. They take as arguments a sub-range among the locally owned elements of the vector, defined as half-open intervals. The intervals are designed to be scheduled close to the time the matrix-vector product touches upon the entries in the src and dst vectors, respectively, with the requirement that

The motivation for this function is to increase data locality and hence cache usage. For the example of a class similar to the one in the step-37 tutorial program, the implementation is

void
const std::function<void(const unsigned int, const unsigned int)>
&operation_before_matrix_vector_product,
const std::function<void(const unsigned int, const unsigned int)>
&operation_after_matrix_vector_product) const
{
data.cell_loop(&LaplaceOperator::local_apply,
this,
dst,
src,
operation_before_matrix_vector_product,
operation_after_matrix_vector_product);
}
void vmult(VectorType &dst, const VectorType &src) const

In terms of the Chebyshev iteration, the operation before the loop will set dst to zero, whereas the operation after the loop performs the iteration leading to \(x^{n+1}\) described above, modifying the dst and src vectors.

Definition at line 2101 of file precondition.h.

Member Typedef Documentation

◆ size_type

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
using PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::size_type = types::global_dof_index

Declare type for container size.

Definition at line 2107 of file precondition.h.

◆ EigenvalueInformation

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
using PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::EigenvalueInformation = internal::EigenvalueInformation

Definition at line 2234 of file precondition.h.

◆ map_value_type

using Subscriptor::map_value_type = decltype(counter_map)::value_type
privateinherited

The data type used in counter_map.

Definition at line 229 of file subscriptor.h.

◆ map_iterator

using Subscriptor::map_iterator = decltype(counter_map)::iterator
privateinherited

The iterator type used in counter_map.

Definition at line 234 of file subscriptor.h.

Constructor & Destructor Documentation

◆ PreconditionChebyshev()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::PreconditionChebyshev ( )

Constructor.

Member Function Documentation

◆ initialize()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::initialize ( const MatrixType &  matrix,
const AdditionalData additional_data = AdditionalData() 
)

Initialize function. Takes the matrix which is used to form the preconditioner, and additional flags if there are any. This function works only if the input matrix has an operator el(i,i) for accessing all the elements in the diagonal. Alternatively, the diagonal can be supplied with the help of the AdditionalData field.

This function calculates an estimate of the eigenvalue range of the matrix weighted by its diagonal using a modified CG iteration in case the given number of iterations is positive.

◆ vmult()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::vmult ( VectorType &  dst,
const VectorType &  src 
) const

Compute the action of the preconditioner on src, storing the result in dst.

◆ Tvmult()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::Tvmult ( VectorType &  dst,
const VectorType &  src 
) const

Compute the action of the transposed preconditioner on src, storing the result in dst.

◆ step()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::step ( VectorType &  dst,
const VectorType &  src 
) const

Perform one step of the preconditioned Richardson iteration.

◆ Tstep()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::Tstep ( VectorType &  dst,
const VectorType &  src 
) const

Perform one transposed step of the preconditioned Richardson iteration.

◆ clear()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
void PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::clear ( )

Resets the preconditioner.

◆ m()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
size_type PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::m ( ) const

Return the dimension of the codomain (or range) space. Note that the matrix is of dimension \(m \times n\).

◆ n()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
size_type PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::n ( ) const

Return the dimension of the domain space. Note that the matrix is of dimension \(m \times n\).

◆ estimate_eigenvalues()

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
EigenvalueInformation PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::estimate_eigenvalues ( const VectorType &  src) const

Compute eigenvalue estimates required for the preconditioner.

This function is called automatically on first use of the preconditioner if it is not called by the user. The layout of the vector src is used to create internal temporary vectors and its content does not matter.

Initializes the factors theta and delta based on an eigenvalue computation. If the user set provided values for the largest eigenvalue in AdditionalData, no computation is performed and the information given by the user is used.

◆ subscribe()

void Subscriptor::subscribe ( std::atomic< bool > *const  validity,
const std::string &  identifier = "" 
) const
inherited

Subscribes a user of the object by storing the pointer validity. The subscriber may be identified by text supplied as identifier.

Definition at line 135 of file subscriptor.cc.

◆ unsubscribe()

void Subscriptor::unsubscribe ( std::atomic< bool > *const  validity,
const std::string &  identifier = "" 
) const
inherited

Unsubscribes a user from the object.

Note
The identifier and the validity pointer must be the same as the one supplied to subscribe().

Definition at line 155 of file subscriptor.cc.

◆ n_subscriptions()

unsigned int Subscriptor::n_subscriptions ( ) const
inlineinherited

Return the present number of subscriptions to this object. This allows to use this class for reference counted lifetime determination where the last one to unsubscribe also deletes the object.

Definition at line 300 of file subscriptor.h.

◆ list_subscribers() [1/2]

template<typename StreamType >
void Subscriptor::list_subscribers ( StreamType &  stream) const
inlineinherited

List the subscribers to the input stream.

Definition at line 317 of file subscriptor.h.

◆ list_subscribers() [2/2]

void Subscriptor::list_subscribers ( ) const
inherited

List the subscribers to deallog.

Definition at line 203 of file subscriptor.cc.

◆ serialize()

template<class Archive >
void Subscriptor::serialize ( Archive &  ar,
const unsigned int  version 
)
inlineinherited

Read or write the data of this object to or from a stream for the purpose of serialization using the BOOST serialization library.

This function does not actually serialize any of the member variables of this class. The reason is that what this class stores is only who subscribes to this object, but who does so at the time of storing the contents of this object does not necessarily have anything to do with who subscribes to the object when it is restored. Consequently, we do not want to overwrite the subscribers at the time of restoring, and then there is no reason to write the subscribers out in the first place.

Definition at line 309 of file subscriptor.h.

◆ check_no_subscribers()

void Subscriptor::check_no_subscribers ( ) const
privatenoexceptinherited

Check that there are no objects subscribing to this object. If this check passes then it is safe to destroy the current object. It this check fails then this function will either abort or print an error message to deallog (by using the AssertNothrow mechanism), but will not throw an exception.

Note
Since this function is just a consistency check it does nothing in release mode.
If this function is called when there is an uncaught exception then, rather than aborting, this function prints an error message to the standard error stream and returns.

Definition at line 52 of file subscriptor.cc.

Member Data Documentation

◆ matrix_ptr

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
SmartPointer< const MatrixType, PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> > PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::matrix_ptr
private

A pointer to the underlying matrix.

Definition at line 2258 of file precondition.h.

◆ solution_old

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
VectorType PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::solution_old
mutableprivate

Internal vector used for the vmult operation.

Definition at line 2263 of file precondition.h.

◆ temp_vector1

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
VectorType PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::temp_vector1
mutableprivate

Internal vector used for the vmult operation.

Definition at line 2268 of file precondition.h.

◆ temp_vector2

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
VectorType PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::temp_vector2
mutableprivate

Internal vector used for the vmult operation.

Definition at line 2273 of file precondition.h.

◆ data

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
AdditionalData PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::data
private

Stores the additional data passed to the initialize function, obtained through a copy operation.

Definition at line 2279 of file precondition.h.

◆ theta

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
double PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::theta
private

Average of the largest and smallest eigenvalue under consideration.

Definition at line 2284 of file precondition.h.

◆ delta

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
double PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::delta
private

Half the interval length between the largest and smallest eigenvalue under consideration.

Definition at line 2290 of file precondition.h.

◆ eigenvalues_are_initialized

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
bool PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::eigenvalues_are_initialized
private

Stores whether the preconditioner has been set up and eigenvalues have been computed.

Definition at line 2296 of file precondition.h.

◆ mutex

template<typename MatrixType = SparseMatrix<double>, typename VectorType = Vector<double>, typename PreconditionerType = DiagonalMatrix<VectorType>>
Threads::Mutex PreconditionChebyshev< MatrixType, VectorType, PreconditionerType >::mutex
mutableprivate

A mutex to avoid that multiple vmult() invocations by different threads overwrite the temporary vectors.

Definition at line 2302 of file precondition.h.

◆ counter

std::atomic<unsigned int> Subscriptor::counter
mutableprivateinherited

Store the number of objects which subscribed to this object. Initially, this number is zero, and upon destruction it shall be zero again (i.e. all objects which subscribed should have unsubscribed again).

The creator (and owner) of an object is counted in the map below if HE manages to supply identification.

We use the mutable keyword in order to allow subscription to constant objects also.

This counter may be read from and written to concurrently in multithreaded code: hence we use the std::atomic class template.

Definition at line 218 of file subscriptor.h.

◆ counter_map

std::map<std::string, unsigned int> Subscriptor::counter_map
mutableprivateinherited

In this map, we count subscriptions for each different identification string supplied to subscribe().

Definition at line 224 of file subscriptor.h.

◆ validity_pointers

std::vector<std::atomic<bool> *> Subscriptor::validity_pointers
mutableprivateinherited

In this vector, we store pointers to the validity bool in the SmartPointer objects that subscribe to this class.

Definition at line 240 of file subscriptor.h.

◆ object_info

const std::type_info* Subscriptor::object_info
mutableprivateinherited

Pointer to the typeinfo object of this object, from which we can later deduce the class name. Since this information on the derived class is neither available in the destructor, nor in the constructor, we obtain it in between and store it here.

Definition at line 248 of file subscriptor.h.


The documentation for this class was generated from the following file: