Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2::Relaxation< MatrixType > Class Template Reference

Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices. More...

#include <Ifpack2_Relaxation_decl.hpp>

Inheritance diagram for Ifpack2::Relaxation< MatrixType >:
[legend]

List of all members.

## Public Member Functions

Constructors and destructors
Relaxation (const Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > &Matrix)
Constructor.
virtual ~Relaxation ()
Destructor.
Preconditioner computation methods
void setParameters (const Teuchos::ParameterList &params)
Set the relaxation / preconditioner parameters.
void initialize ()
Initialize the preconditioner.
bool isInitialized () const
Returns true if the preconditioner has been successfully initialized.
void compute ()
Compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation parameters.
bool isComputed () const
Return true if compute() has been called.
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.
const 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.
const 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.
bool hasTransposeApply () const
Whether apply() and applyMat() let you apply the transpose or conjugate transpose.
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 preconditioner to X, returning the result in Y.
Mathematical functions
magnitude_type computeCondEst (CondestType CT=Cheap, local_ordinal_type MaxIters=1550, magnitude_type Tol=1e-9, const Teuchos::Ptr< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > &matrix=Teuchos::null)
Computes and returns the estimated condition number.
Attribute accessor methods
magnitude_type getCondEst () const
The computed estimated condition number, or -1 if not previously computed.
const Teuchos::RCP< const
Teuchos::Comm< int > > &
getComm () const
The communicator over which the matrix and vectors are distributed.
Teuchos::RCP< const
Tpetra::RowMatrix< scalar_type,
local_ordinal_type,
global_ordinal_type, node_type > >
getMatrix () const
The matrix to be preconditioned.
double getComputeFlops () const
Total number of floating-point operations over all calls to compute().
double getApplyFlops () const
Total number of floating-point operations over all calls to apply().
int getNumInitialize () const
Total number of calls to initialize().
int getNumCompute () const
Total number of calls to compute().
int getNumApply () const
Total number of calls to apply().
double getInitializeTime () const
Total time in seconds spent in all calls to initialize().
double getComputeTime () const
Total time in seconds spent in all calls to compute().
double getApplyTime () const
Total time in seconds spent in all calls to apply().
Implementation of Teuchos::Describable interface
std::string description () const
A simple one-line description of this object.
void describe (Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to a Teuchos::FancyOStream.

## Typedefs

typedef MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
typedef
MatrixType::local_ordinal_type
local_ordinal_type
The type of local indices in the input MatrixType.
typedef
MatrixType::global_ordinal_type
global_ordinal_type
The type of global indices in the input MatrixType.
typedef MatrixType::node_type node_type
The type of the Kokkos Node used by the input MatrixType.
typedef Teuchos::ScalarTraits
< scalar_type >::magnitudeType
magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
TEUCHOS_DEPRECATED typedef
MatrixType::scalar_type
Scalar
Preserved only for backwards compatibility. Please use "scalar_type".
TEUCHOS_DEPRECATED typedef
MatrixType::local_ordinal_type
LocalOrdinal
Preserved only for backwards compatibility. Please use "local_ordinal_type".
TEUCHOS_DEPRECATED typedef
MatrixType::global_ordinal_type
GlobalOrdinal
Preserved only for backwards compatibility. Please use "global_ordinal_type".
TEUCHOS_DEPRECATED typedef
MatrixType::node_type
Node
Preserved only for backwards compatibility. Please use "node_type".
TEUCHOS_DEPRECATED typedef
Teuchos::ScalarTraits
< scalar_type >::magnitudeType
magnitudeType
Preserved only for backwards compatibility. Please use "magnitude_type".

## Detailed Description

### template<class MatrixType> class Ifpack2::Relaxation< MatrixType >

Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.

Template Parameters:
 MatrixType A specialization of Tpetra::CrsMatrix.

## Summary

This class implements several different relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix. Relaxation is derived from Preconditioner, which is itself derived from Tpetra::Operator. Therefore this object can be used as a preconditioner everywhere an apply() method is required in the preconditioning step.

This class implements the following relaxation methods:

• Jacobi
• Gauss-Seidel
• Symmetric Gauss-Seidel

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.

## Performance

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 Kokkos 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 threaded Kokkos Node type.

Relaxation works with any Tpetra::RowMatrix. 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 normally don't have to do anything to figure this out (we test via dynamic_cast), but it may help to use a Tpetra::CrsMatrix specialization as the MatrixType template parameter, rather than a Tpetra::RowMatrix specialization. (This matters if you are using a nondefault value of the fifth template parameter of Tpetra::CrsMatrix.)

## Creating a Relaxation preconditioner

The following code snippet shows how to create a Relaxation preconditioner.

#include "Ifpack2_Relaxation.hpp"

...
using Teuchos::ParameterList;
using Teuchos::RCP;
typedef double ST;
typedef int    LO;
typedef int    GO;
typedef Tpetra::CrsMatrix<ST, LO, GO> crs_matrix_type;
typedef Ifpack2::Preconditioner<ST, LO, GO> precond_type;
...

// Create the sparse matrix A somehow.  It must be fill complete
// before you may create an Ifpack2 preconditioner from it.
RCP<crs_matrix_type> A = ...;

// Create the relaxation.  You could also do this using
// Ifpack2::Factory (the preconditioner factory) if you like.
precond_type prec (A);

// Make the list of relaxation parameters.
Teuchos::ParameterList params;
// Do symmetric SOR / Gauss-Seidel.
params.set ("relaxation: type", "Symmetric Gauss-Seidel");
// Two sweeps (of symmetric SOR / Gauss-Seidel) per apply() call.
params.set ("relaxation: sweeps", 2);
// ... Set any other parameters you want to set ...

// Set parameters.
prec.setParameters (params);

// Prepare the relaxation instance for use.
prec.initialize ();
prec.compute ();

// Now prec may be used as a preconditioner or smoother,
// by calling its apply() method, just like any Tpetra::Operator.


## Algorithms

We now briefly describe the relaxation algorithms this class implements. Consider a linear system of type

$A x = b,$

where $$A$$ is a square matrix, and $$x$$, $$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$. Here, $$x^{(k)}_i$$ is the $i$-th element of vector $$x^{(k)}$$.

The Jacobi method computes

$x^{(k+1)}_i = A_{ii}^{-1} ( b_i - \sum_{j != 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 != 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's 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.

Note that this class does not actually use the formulae above to apply Jacobi or SOR. For example, we optimize the formulae to avoid branches.

## Member Typedef Documentation

template<class MatrixType>
 typedef MatrixType::scalar_type Ifpack2::Relaxation< MatrixType >::scalar_type

The type of the entries of the input MatrixType.

template<class MatrixType>
 typedef MatrixType::local_ordinal_type Ifpack2::Relaxation< MatrixType >::local_ordinal_type

The type of local indices in the input MatrixType.

template<class MatrixType>
 typedef MatrixType::global_ordinal_type Ifpack2::Relaxation< MatrixType >::global_ordinal_type

The type of global indices in the input MatrixType.

template<class MatrixType>
 typedef MatrixType::node_type Ifpack2::Relaxation< MatrixType >::node_type

The type of the Kokkos Node used by the input MatrixType.

template<class MatrixType>
 typedef Teuchos::ScalarTraits::magnitudeType Ifpack2::Relaxation< MatrixType >::magnitude_type

The type of the magnitude (absolute value) of a matrix entry.

## Constructor & Destructor Documentation

template<class MatrixType >
 Ifpack2::Relaxation< MatrixType >::Relaxation ( const Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > & Matrix )  [explicit]

Constructor.

Parameters:
 Matrix [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. In particular, the compute() method may extract the diagonal entries and precompute their inverses, in order to speed up Gauss-Seidel or to implement the L1 version of various relaxation methods. If you plan to change the diagonal entries of the matrix after making a Relaxation instance with that matrix, you must destroy the old Relaxation instance and create a new one after changing the diagonal entries.

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):

 RCP<const CrsMatrix<...> > A = ...;
Relaxation<CrsMatrix<...> > R (A);


but you may not do this:

 void foo (const Relaxation<CrsMatrix<...> >& R);

RCP<const CrsMatrix<...> > A = ...;
foo (A);

template<class MatrixType >
 Ifpack2::Relaxation< MatrixType >::~Relaxation ( )  [virtual]

Destructor.

## Member Function Documentation

template<class MatrixType >
 void Ifpack2::Relaxation< MatrixType >::setParameters ( const Teuchos::ParameterList & params )  [virtual]

Set the relaxation / preconditioner parameters.

Warning:
All parameters are case sensitive. We make no attempt to check the spelling of parameter names in your input list.

The "relaxation: type" (string) parameter sets the relaxation / preconditioner method you want to use. It currently accepts the following values (the default is "Jacobi"):

• "Jacobi"
• "Gauss-Seidel"
• "Symmetric Gauss-Seidel"

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: min diagonal value" (scalar_type) 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, meaning that we do not impose a minimum diagonal value by default.

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 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.

template<class MatrixType >
 void Ifpack2::Relaxation< MatrixType >::initialize ( )  [virtual]

Initialize the preconditioner.

template<class MatrixType>
 bool Ifpack2::Relaxation< MatrixType >::isInitialized ( ) const [inline, virtual]

Returns true if the preconditioner has been successfully initialized.

template<class MatrixType >
 void Ifpack2::Relaxation< MatrixType >::compute ( )  [virtual]

Compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation parameters.

template<class MatrixType>
 bool Ifpack2::Relaxation< MatrixType >::isComputed ( ) const [inline, virtual]

Return true if compute() has been called.

template<class MatrixType>
 void Ifpack2::Relaxation< MatrixType >::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 [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.

Parameters:
 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.
template<class MatrixType >
 const Teuchos::RCP< const Tpetra::Map< typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > & Ifpack2::Relaxation< MatrixType >::getDomainMap ( ) const [virtual]

Returns the Tpetra::Map object associated with the domain of this operator.

template<class MatrixType >
 const Teuchos::RCP< const Tpetra::Map< typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > & Ifpack2::Relaxation< MatrixType >::getRangeMap ( ) const [virtual]

Returns the Tpetra::Map object associated with the range of this operator.

template<class MatrixType >
 bool Ifpack2::Relaxation< MatrixType >::hasTransposeApply ( ) const

Whether apply() and applyMat() let you apply the transpose or conjugate transpose.

template<class MatrixType>
 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 preconditioner to X, returning the result in Y.

This method computes Y = M*X, where M*X represents the action of the preconditioner on the input multivector X.

Parameters:
 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.
template<class MatrixType>
 magnitude_type Ifpack2::Relaxation< MatrixType >::computeCondEst ( CondestType CT = Cheap, local_ordinal_type MaxIters = 1550, magnitude_type Tol = 1e-9, const Teuchos::Ptr< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > & matrix = Teuchos::null )  [virtual]

Computes and returns the estimated condition number.

We use an iterative process to estimate the condition number. You can control the number of iterations we use, the iteration tolerance, and how hard to work at estimating.

Parameters:
 CondestType [in] How hard to work at estimating the condition number. Cheap means not very hard. MaxIters [in] Maximum number of iterations for estimating the condition number. Tol [in] Iteration tolerance.
template<class MatrixType >
 Teuchos::ScalarTraits< typename MatrixType::scalar_type >::magnitudeType Ifpack2::Relaxation< MatrixType >::getCondEst ( ) const [virtual]

The computed estimated condition number, or -1 if not previously computed.

If you haven't yet called computeCondEst(), then this method returns -1.

template<class MatrixType >
 const Teuchos::RCP< const Teuchos::Comm< int > > & Ifpack2::Relaxation< MatrixType >::getComm ( ) const

The communicator over which the matrix and vectors are distributed.

template<class MatrixType >
 Teuchos::RCP< const Tpetra::RowMatrix< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > Ifpack2::Relaxation< MatrixType >::getMatrix ( ) const [virtual]

The matrix to be preconditioned.

template<class MatrixType >
 double Ifpack2::Relaxation< MatrixType >::getComputeFlops ( ) const

Total number of floating-point operations over all calls to compute().

template<class MatrixType >
 double Ifpack2::Relaxation< MatrixType >::getApplyFlops ( ) const

Total number of floating-point operations over all calls to apply().

template<class MatrixType >
 int Ifpack2::Relaxation< MatrixType >::getNumInitialize ( ) const [virtual]

Total number of calls to initialize().

template<class MatrixType >
 int Ifpack2::Relaxation< MatrixType >::getNumCompute ( ) const [virtual]

Total number of calls to compute().

template<class MatrixType >
 int Ifpack2::Relaxation< MatrixType >::getNumApply ( ) const [virtual]

Total number of calls to apply().

template<class MatrixType >
 double Ifpack2::Relaxation< MatrixType >::getInitializeTime ( ) const [virtual]

Total time in seconds spent in all calls to initialize().

template<class MatrixType >
 double Ifpack2::Relaxation< MatrixType >::getComputeTime ( ) const [virtual]

Total time in seconds spent in all calls to compute().

template<class MatrixType >
 double Ifpack2::Relaxation< MatrixType >::getApplyTime ( ) const [virtual]

Total time in seconds spent in all calls to apply().

template<class MatrixType >
 std::string Ifpack2::Relaxation< MatrixType >::description ( ) const

A simple one-line description of this object.

template<class MatrixType >
 void Ifpack2::Relaxation< MatrixType >::describe ( Teuchos::FancyOStream & out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default ) const

Print the object with some verbosity level to a Teuchos::FancyOStream.

## Member Data Documentation

template<class MatrixType>
 TEUCHOS_DEPRECATED typedef MatrixType::scalar_type Ifpack2::Relaxation< MatrixType >::Scalar

Preserved only for backwards compatibility. Please use "scalar_type".

template<class MatrixType>
 TEUCHOS_DEPRECATED typedef MatrixType::local_ordinal_type Ifpack2::Relaxation< MatrixType >::LocalOrdinal

Preserved only for backwards compatibility. Please use "local_ordinal_type".

template<class MatrixType>
 TEUCHOS_DEPRECATED typedef MatrixType::global_ordinal_type Ifpack2::Relaxation< MatrixType >::GlobalOrdinal

Preserved only for backwards compatibility. Please use "global_ordinal_type".

template<class MatrixType>
 TEUCHOS_DEPRECATED typedef MatrixType::node_type Ifpack2::Relaxation< MatrixType >::Node

Preserved only for backwards compatibility. Please use "node_type".

template<class MatrixType>
 TEUCHOS_DEPRECATED typedef Teuchos::ScalarTraits::magnitudeType Ifpack2::Relaxation< MatrixType >::magnitudeType

Preserved only for backwards compatibility. Please use "magnitude_type".

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