 |
Reference documentation for deal.II version 9.4.1
|
Loading...
Searching...
No Matches
|
Reference documentation for deal.II version 9.4.1
|
#include <deal.II/lac/precondition.h>
Classes | |
| struct | AdditionalData |
| struct | EigenvalueInformation |
Public Types | |
| using | size_type = types::global_dof_index |
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 |
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 |
Subscriptor functionality | |
Classes derived from Subscriptor provide a facility to subscribe to this object. This is mostly used by the SmartPointer class. | |
| std::atomic< unsigned int > | counter |
| std::map< std::string, unsigned int > | counter_map |
| std::vector< std::atomic< bool > * > | validity_pointers |
| const std::type_info * | object_info |
| 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 |
| template<class Archive > | |
| void | serialize (Archive &ar, const unsigned int version) |
| using | map_value_type = decltype(counter_map)::value_type |
| using | map_iterator = decltype(counter_map)::iterator |
| static ::ExceptionBase & | ExcInUse (int arg1, std::string arg2, std::string arg3) |
| static ::ExceptionBase & | ExcNoSubscriber (std::string arg1, std::string arg2) |
| void | check_no_subscribers () const noexcept |
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.
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.
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.
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
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.
MatrixType argumentThis 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
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
src or dst once the operation_before_matrix_vector_product has been run on that vector entry; operation_after_matrix_vector_product may run on a range of entries [i,j) once the matrix-vector product does not access the entries [i,j) in src and dst any more. 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
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 1594 of file precondition.h.
1.9.5