|
Ifpack2 Templated Preconditioning Package
Version 1.0
|
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
|
Ifpack2 Templated Preconditioning Package
Version 1.0
|
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices. More...
#include <Ifpack2_Relaxation_decl.hpp>
Public Types | |
Typedefs | |
| typedef MatrixType::scalar_type | scalar_type |
| The type of the entries of the input MatrixType. More... | |
| typedef MatrixType::local_ordinal_type | local_ordinal_type |
| The type of local indices in the input MatrixType. More... | |
| typedef MatrixType::global_ordinal_type | global_ordinal_type |
| The type of global indices in the input MatrixType. More... | |
| typedef MatrixType::node_type | node_type |
| The Node type used by the input MatrixType. More... | |
| typedef MatrixType::node_type::device_type | device_type |
| The Kokkos device type used by the input MatrixType. More... | |
| typedef Teuchos::ScalarTraits < scalar_type >::magnitudeType | magnitude_type |
| The type of the magnitude (absolute value) of a matrix entry. More... | |
| typedef Tpetra::RowMatrix < scalar_type, local_ordinal_type, global_ordinal_type, node_type > | row_matrix_type |
| Tpetra::RowMatrix specialization used by this class. More... | |
| typedef Tpetra::Operator < scalar_type, local_ordinal_type, global_ordinal_type, node_type > | op_type |
| Tpetra::Operator specialization used by this class. More... | |
 Public Types inherited from Ifpack2::Preconditioner< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > | |
| typedef Teuchos::ScalarTraits < MatrixType::scalar_type > ::magnitudeType | magnitude_type |
| The type of the magnitude (absolute value) of a matrix entry. More... | |
Public Member Functions | |
Constructors and destructors | |
| Relaxation (const Teuchos::RCP< const row_matrix_type > &A) | |
| Constructor. More... | |
| virtual | ~Relaxation ()=default |
| Destructor. More... | |
Preconditioner computation methods | |
| void | setParameters (const Teuchos::ParameterList ¶ms) |
| Set the relaxation / preconditioner parameters. More... | |
| bool | supportsZeroStartingSolution () |
| void | setZeroStartingSolution (bool zeroStartingSolution) |
| Set this preconditioner's parameters. More... | |
| Teuchos::RCP< const Teuchos::ParameterList > | getValidParameters () const |
| Return a list of all the parameters that this class accepts. More... | |
| void | initialize () |
| Initialize the preconditioner ("symbolic setup"). More... | |
| bool | isInitialized () const |
Returns true if the preconditioner has been successfully initialized. More... | |
| void | compute () |
| Compute the preconditioner ("numeric setup");. More... | |
| bool | isComputed () const |
| Return true if compute() has been called. More... | |
Implementation of Ifpack2::Details::CanChangeMatrix | |
| virtual void | setMatrix (const Teuchos::RCP< const row_matrix_type > &A) |
| Change the matrix to be preconditioned. More... | |
Implementation of the Tpetra::Operator interface | |
| void | apply (const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const |
| Apply the preconditioner to X, returning the result in Y. More... | |
| Teuchos::RCP< const Tpetra::Map < local_ordinal_type, global_ordinal_type, node_type > > | getDomainMap () const |
| Returns the Tpetra::Map object associated with the domain of this operator. More... | |
| Teuchos::RCP< const Tpetra::Map < local_ordinal_type, global_ordinal_type, node_type > > | getRangeMap () const |
| Returns the Tpetra::Map object associated with the range of this operator. More... | |
| bool | hasTransposeApply () const |
| Whether apply() and applyMat() let you apply the transpose or conjugate transpose. More... | |
| void | applyMat (const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const |
| Apply the input matrix to X, returning the result in Y. More... | |
Attribute accessor methods | |
| Teuchos::RCP< const Teuchos::Comm< int > > | getComm () const |
| The communicator over which the matrix and vectors are distributed. More... | |
| Teuchos::RCP< const row_matrix_type > | getMatrix () const |
| The matrix to be preconditioned. More... | |
| double | getComputeFlops () const |
| Total number of floating-point operations over all calls to compute(). More... | |
| double | getApplyFlops () const |
| Total number of floating-point operations over all calls to apply(). More... | |
| int | getNumInitialize () const |
| Total number of calls to initialize(). More... | |
| int | getNumCompute () const |
| Total number of calls to compute(). More... | |
| int | getNumApply () const |
| Total number of calls to apply(). More... | |
| double | getInitializeTime () const |
| Total time in seconds spent in all calls to initialize(). More... | |
| double | getComputeTime () const |
| Total time in seconds spent in all calls to compute(). More... | |
| double | getApplyTime () const |
| Total time in seconds spent in all calls to apply(). More... | |
| size_t | getNodeSmootherComplexity () const |
| Get a rough estimate of cost per iteration. More... | |
Implementation of Teuchos::Describable interface | |
| std::string | description () const |
| A simple one-line description of this object. More... | |
| void | describe (Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const |
| Print the object's attributes to the given output stream. More... | |
| Public Member Functions inherited from Ifpack2::Preconditioner< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > | |
| virtual | ~Preconditioner () |
| Destructor. More... | |
| Public Member Functions inherited from Ifpack2::Details::CanChangeMatrix< Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > | |
| virtual void | setMatrix (const Teuchos::RCP< const Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > &A)=0 |
| Set the new matrix. More... | |
| virtual | ~CanChangeMatrix () |
| Destructor. More... | |
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
| MatrixType | A specialization of Tpetra::RowMatrix. |
This class implements several different relaxation preconditioners for Tpetra::RowMatrix or its subclass Tpetra::CrsMatrix. Relaxation is derived from Preconditioner, which is itself derived from Tpetra::Operator. Therefore this object may be used as a preconditioner for Belos linear solvers, and for any linear solver that treats preconditioners as instances of Tpetra::Operator.
This class implements the following relaxation methods:
All methods allow you to set an optional damping parameter. The "Gauss-Seidel" methods technically only perform Gauss-Seidel within an MPI process, but Jacobi between processes. To compensate, these methods include an "L1" option, which can improve convergence by weighting contributions near process boundaries differently. For more details, please refer to the following publication:
A. H. Baker, R. D. Falgout, T. V. Kolev, and U. M. Yang. "Multigrid Smoothers for Ultraparallel Computing." SIAM J. Sci. Comput., Vol. 33, No. 5. (2011), pp. 2864-2887.
Richardson and Jacobi will always use your matrix's native sparse matrix-vector multiply kernel. This should give good performance, since we have spent a lot of effort tuning Tpetra's kernels. Depending on the Node type of your Tpetra matrix, it may also exploit threads for additional parallelism within each MPI process. In contrast, Gauss-Seidel and symmetric Gauss-Seidel are intrinsically sequential methods within an MPI process. This prevents us from exposing more parallelism via threads. The difference should become more apparent as your code moves away from a "one MPI process per core" model, to a "one MPI process per socket or node" model, assuming that you are using a thread-parallel Node type.
Relaxation works with any Tpetra::RowMatrix input. If your Tpetra::RowMatrix happens to be a Tpetra::CrsMatrix, the Gauss-Seidel and symmetric Gauss-Seidel relaxations may be able to exploit this for better performance. You don't have to do anything to figure this out (we test via dynamic_cast).
The following code snippet shows how to create a Relaxation preconditioner.
We now briefly describe the relaxation algorithms this class implements. Consider the linear system \(Ax=b\), where \(A\) is a square matrix, and \(x\) and \(b\) are two vectors of compatible dimensions. Suppose that \(x^{(0)}\) is the starting vector and \(x^{(k)}\) is the approximate solution for \(x\) computed by iteration $k+1$ of whatever relaxation method we are using. Here, \(x^{(k)}_i\) is the $i$-th element of vector \(x^{(k)}\).
The Richardson method computes
\[ x^{(k+1)}_i = x_^{(k)}_i + alpha ( b_i - \sum_{j} A_{ij} x^{(k)}_j ). \]
The Jacobi method computes
\[ x^{(k+1)}_i = A_{ii}^{-1} ( b_i - \sum_{j \neq i} A_{ij} x^{(k)}_j ). \]
The "damped" Jacobi method generalizes Jacobi. It introduces a damping parameter \(\omega \), and computes
\[ x^{(k+1)}_i = (1 - \omega) x^{(k)}_i + \omega A_{ii}^{-1} ( b_i - \sum_{j \neq i} A_{ij} x^{(k)}_j ). \]
The "damped Gauss-Seidel method" is actually successive over-relaxation (SOR), with Gauss-Seidel as a special case when the damping parameter \(\omega = 1\). We implement has two different sweep directions: Forward and Backward. The Forward sweep direction computes
\[ x^{(k+1)}_i = (1 - \omega) x^{(k)}_i + \omega A_{ii}^{-1} ( b_i - \sum_{j < i} A_{ij} x^{(k+1)}_j - \sum_{j > i} A_{ij} x^{(k)}_j ), \]
and the Backward sweep direction computes
\[ x^{(k+1)}_i = (1 - \omega) x^{(k)}_i + \omega A_{ii}^{-1} ( b_i - \sum_{j > i} A_{ij} x^{(k+1)}_j - \sum_{j < i} A_{ij} x^{(k)}_j ), \]
Users may set the sweep direction via the "relaxation: backward mode" option. See the documentation of setParameters() for details.
Gauss-Seidel / SOR also comes in a symmetric version. This method first does a Forward sweep, then a Backward sweep. Only the symmetric version of this preconditioner is guaranteed to be symmetric (or Hermitian, if the matrix data are complex).
Users may set the relaxation method via the "relaxation: type" parameter. For all relaxation methods, users may specify the number of sweeps per call to apply() and the damping parameter \(\omega \). For a list of all supported parameters, please refer to the documentation of the setParameters() method. For advice on picking \(\omega \) for a preconditioner, please refer to the following book: "Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition," R. Barrett et al., SIAM, 1994.
\[ x^{(k+1)}_i = x^{(k)}_i + \omega b_i - \frac{\omega}{A_{ii}} ( \sum_{j \geq i} A_{ij} x^{(k)}_j + \sum_{j < i} x^{(k+1)}_j ). \]
Executing this expression in a forward sweep does not require distinguishing between the lower and upper triangle of A. The same thing holds for the backward sweep.| typedef MatrixType::scalar_type Ifpack2::Relaxation< MatrixType >::scalar_type |
The type of the entries of the input MatrixType.
| typedef MatrixType::local_ordinal_type Ifpack2::Relaxation< MatrixType >::local_ordinal_type |
The type of local indices in the input MatrixType.
| typedef MatrixType::global_ordinal_type Ifpack2::Relaxation< MatrixType >::global_ordinal_type |
The type of global indices in the input MatrixType.
| typedef MatrixType::node_type Ifpack2::Relaxation< MatrixType >::node_type |
The Node type used by the input MatrixType.
| typedef MatrixType::node_type::device_type Ifpack2::Relaxation< MatrixType >::device_type |
The Kokkos device type used by the input MatrixType.
| typedef Teuchos::ScalarTraits<scalar_type>::magnitudeType Ifpack2::Relaxation< MatrixType >::magnitude_type |
The type of the magnitude (absolute value) of a matrix entry.
| typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> Ifpack2::Relaxation< MatrixType >::row_matrix_type |
Tpetra::RowMatrix specialization used by this class.
| typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> Ifpack2::Relaxation< MatrixType >::op_type |
Tpetra::Operator specialization used by this class.
|
explicit |
Constructor.
| A | [in] The matrix for which to make the constructor. Tpetra::RowMatrix is the base class of Tpetra::CrsMatrix, so you may give either a Tpetra::RowMatrix or a Tpetra::CrsMatrix here. |
The results of apply() are undefined if you change the diagonal entries of the sparse matrix after invoking this constructor, without first calling compute(). In particular, the compute() method may extract the diagonal entries and precompute their inverses, in order to speed up any of the relaxation methods that this class implements.
The "explicit" keyword just means that you must invoke the Relaxation constructor explicitly; you aren't allowed to use it as an implicit conversion ("cast"). For example, you may do this (namespaces and Tpetra template parameters omitted for brevity):
but you may not do this:
|
virtualdefault |
Destructor.
|
virtual |
Set the relaxation / preconditioner parameters.
The "relaxation: type" (string) parameter sets the relaxation / preconditioner method you want to use. It currently accepts the following values (the default is "Jacobi"):
The "relaxation: sweeps" (int) parameter sets the number of sweeps, that is, the number of times to apply the relaxation on each invocation of apply(). The default number of sweeps is 1.
The "relaxation: damping factor" (scalar_type – the type of the entries of the matrix) parameter is the value of the damping factor \(\omega \). The main documentation of this class explains how we use this value. The default value is 1.0.
The "relaxation: zero starting solution" (bool) parameter governs whether or not we use the existing values in the output multivector Y when applying the relaxation. Its default value is true, meaning that we fill Y with zeros before applying relaxation sweeps. If false, we use the existing values in Y.
If the "relaxation: backward mode" (bool) parameter is true, we perform Gauss-Seidel in reverse mode. The default value is false, meaning that we do forward-mode Gauss-Seidel. This only affects standard Gauss-Seidel, not symmetric Gauss-Seidel.
The "relaxation: fix tiny diagonal entries" (bool) parameter defaults to false. If true, the compute() method will do extra work (computation only, no MPI communication) to "fix" diagonal entries that are less than or equal to the threshold given by the (magnitude of the) "relaxation: min diagonal value" parameter. The default behavior imitates that of Aztec, which does not do any special modification of the diagonal.
The "relaxation: min diagonal value" (scalar_type) parameter only matters if "relaxation: fix tiny diagonal entries" (see above) is true. This parameter limits how close to zero the diagonal elements of the matrix are allowed to be. If the magnitude of a diagonal element of the matrix is less than the magnitude of this value, then we set that diagonal element to this value. (We don't actually modify the matrix; we just remember the diagonal values.) The use of magnitude rather than the value itself makes this well defined if scalar_type is complex. The default value of this parameter is zero, in which case we will replace diagonal entries that are exactly equal to zero with a small nonzero value (machine precision for the given Scalar type) before inverting them. Note that if "relaxation: fix tiny diagonal entries" is false, the default value, this parameter does nothing.)
The "relaxation: check diagonal entries" (bool) parameter defaults to false. If true, the compute() method will do extra work (both computation and communication) to count diagonal entries that are zero, have negative real part, or are small in magnitude. The describe() method will then print this information for you. You may find this useful for checking whether your input matrix has issues that make Jacobi or Gauss-Seidel a poor choice of preconditioner.
The last two parameters govern the L1 variant of Gauss-Seidel. The "relaxation: use l1" (bool) parameter, if true, turns on the L1 variant. (In "l1", the first character is a lower-case L, and the second character is the numeral 1 (one).) This parameter's value is false by default. The "relaxation: l1 eta" (magnitude_type) parameter is the \(\eta \) parameter associated with that method; its default value is 1.5. Recall that "magnitude_type" is the type of the absolute value of a scalar_type value. This is the same as scalar_type for real-valued floating-point types (like float and double). If scalar_type is std::complex<T> for some type T, then magnitude_type is T.
|
inlinevirtual |
Set this preconditioner's parameters.
Reimplemented from Ifpack2::Preconditioner< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type >.
| Teuchos::RCP< const Teuchos::ParameterList > Ifpack2::Relaxation< MatrixType >::getValidParameters | ( | ) | const |
Return a list of all the parameters that this class accepts.
|
virtual |
Initialize the preconditioner ("symbolic setup").
|
inlinevirtual |
Returns true if the preconditioner has been successfully initialized.
|
virtual |
Compute the preconditioner ("numeric setup");.
|
inlinevirtual |
Return true if compute() has been called.
|
virtual |
Change the matrix to be preconditioned.
| A | [in] The new matrix. |
! isInitialized () ! isComputed ()Calling this method with a matrix different than the current matrix resets the preconditioner's state. After calling this method with a nonnull input, you must first call initialize() and compute() (in that order) before you may call apply().
You may call this method with a null input. If A is null, then you may not call initialize() or compute() until you first call this method again with a nonnull input. This method invalidates any previous factorization whether or not A is null, so calling setMatrix() with a null input is one way to clear the preconditioner's state (and free any memory that it may be using).
The new matrix A need not necessarily have the same Maps or even the same communicator as the original matrix.
|
virtual |
Apply the preconditioner to X, returning the result in Y.
This method computes Y = beta*Y + alpha*M*X, where M*X represents the action of the preconditioner on the input multivector X.
| X | [in] The multivector input of the preconditioner. |
| Y | [in/out] The multivector output of the preconditioner. |
| mode | [in] Whether to apply the transpose or conjugate transpose of the preconditioner. Not all preconditioners support options other than the default (no transpose); please call hasTransposeApply() to determine whether nondefault options are supported. |
| alpha | [in] Scaling factor for the preconditioned input. |
| beta | [in] Scaling factor for the output. |
|
virtual |
Returns the Tpetra::Map object associated with the domain of this operator.
|
virtual |
Returns the Tpetra::Map object associated with the range of this operator.
| bool Ifpack2::Relaxation< MatrixType >::hasTransposeApply | ( | ) | const |
Whether apply() and applyMat() let you apply the transpose or conjugate transpose.
| void Ifpack2::Relaxation< MatrixType >::applyMat | ( | const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > & | X, |
| Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > & | Y, | ||
| Teuchos::ETransp | mode = Teuchos::NO_TRANS |
||
| ) | const |
Apply the input matrix to X, returning the result in Y.
| X | [in] The multivector input of the preconditioner. |
| Y | [in/out] The multivector output of the preconditioner. |
| mode | [in] Whether to apply the transpose or conjugate transpose of the matrix. |
| Teuchos::RCP< const Teuchos::Comm< int > > Ifpack2::Relaxation< MatrixType >::getComm | ( | ) | const |
The communicator over which the matrix and vectors are distributed.
|
virtual |
The matrix to be preconditioned.
| double Ifpack2::Relaxation< MatrixType >::getComputeFlops | ( | ) | const |
Total number of floating-point operations over all calls to compute().
| double Ifpack2::Relaxation< MatrixType >::getApplyFlops | ( | ) | const |
Total number of floating-point operations over all calls to apply().
|
virtual |
Total number of calls to initialize().
|
virtual |
Total number of calls to compute().
|
virtual |
Total number of calls to apply().
|
virtual |
Total time in seconds spent in all calls to initialize().
|
virtual |
Total time in seconds spent in all calls to compute().
|
virtual |
Total time in seconds spent in all calls to apply().
| size_t Ifpack2::Relaxation< MatrixType >::getNodeSmootherComplexity | ( | ) | const |
Get a rough estimate of cost per iteration.
| std::string Ifpack2::Relaxation< MatrixType >::description | ( | ) | const |
A simple one-line description of this object.
Be aware that this will print a very long line, because some users really like to see all the attributes of the object in a single line. If you prefer multiple lines of output, you should call describe() instead.
| void Ifpack2::Relaxation< MatrixType >::describe | ( | Teuchos::FancyOStream & | out, |
| const Teuchos::EVerbosityLevel | verbLevel = Teuchos::Describable::verbLevel_default |
||
| ) | const |
Print the object's attributes to the given output stream.
This method will print a constant amount of information (not proportional to the matrix's dimensions or number of entries) on Process 0 of the communicator over which this object is distributed.
You may wrap an std::ostream in a Teuchos::FancyOStream by including "Teuchos_FancyOStream.hpp" and calling Teuchos::getFancyOStream(). For example:
1.8.5