Left-scaled Chebyshev iteration. More...

#include <Ifpack2_Details_Chebyshev_decl.hpp>

Inheritance diagram for Ifpack2::Details::Chebyshev< ScalarType, MV >:

Public Types
Typedefs
typedef ScalarType	ST
	The type of entries in the matrix and vectors. More...

typedef Teuchos::ScalarTraits < ScalarType >	STS
	Traits class for ST. More...

typedef STS::magnitudeType	MT
	The type of the absolute value of a ScalarType. More...

typedef Tpetra::Operator < typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type >	op_type
	Specialization of Tpetra::Operator. More...

typedef Tpetra::RowMatrix < typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type >	row_matrix_type
	Specialization of Tpetra::RowMatrix. More...

typedef Tpetra::Vector < typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type >	V
	Specialization of Tpetra::Vector. More...

typedef Tpetra::Map< typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type >	map_type
	Specialization of Tpetra::Map. More...

Public Member Functions
	Chebyshev (Teuchos::RCP< const row_matrix_type > A)

	Chebyshev (Teuchos::RCP< const row_matrix_type > A, Teuchos::ParameterList &params)

void	setParameters (Teuchos::ParameterList &plist)
	Set (or reset) parameters. More...

void	compute ()
	(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A. More...

MT	apply (const MV &B, MV &X)
	Solve Ax=b for x with Chebyshev iteration with left diagonal scaling. More...

Teuchos::RCP< const row_matrix_type >	getMatrix () const
	Get the matrix given to the constructor. More...

void	setMatrix (const Teuchos::RCP< const row_matrix_type > &A)
	Set the matrix. More...

bool	hasTransposeApply () const
	Whether it's possible to apply the transpose of this operator. More...

void	print (std::ostream &out)
	Print instance data to the given output stream. More...

Implementation of Teuchos::Describable
std::string	description () const
	A single-line description of the Chebyshev solver. More...

void	describe (Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
	Print a description of the Chebyshev solver to `out`. More...

Detailed Description

template<class ScalarType, class MV>
class Ifpack2::Details::Chebyshev< ScalarType, MV >

Left-scaled Chebyshev iteration.

Template Parameters

ScalarType	The type of entries in the matrix and vectors.
MV	Specialization of Tpetra::MultiVector.

Warning: This class is NOT for users; it is an implementation detail of Ifpack2. Users should use Ifpack2::Chebyshev instead.

This class implements two variants of Chebyshev iteration:

A direct imitation of Ifpack's implementation
A textbook version of the algorithm
Chebyshev polynomials of the 4th kind, using optimal coefficients

All implemented variants use the diagonal of the matrix to precondition the linear system on the left. Diagonal entries less than machine precision are replaced with machine precision.

The first version imitates Ifpack_Chebyshev, both in how it sets parameters and in the actual iteration (ApplyInverse()). The "textbook" in variant #2 above is "Templates for the Solution of Linear Systems," 2nd edition. Experiments show that the Ifpack imitation is much less sensitive to the eigenvalue bounds than the textbook version, so users should prefer it. (In fact, it is the default.) Variant #3 is an experimental implementation of Chebyshev polynomials of the 4th kind with optimal coefficients, from https://arxiv.org/pdf/2202.08830.pdf.

We require that the matrix A be real valued and symmetric positive definite. If users could provide the ellipse parameters ("d" and "c" in the literature, where d is the real-valued center of the ellipse, and d-c and d+c the two foci), the iteration itself would work fine with nonsymmetric real-valued A, as long as the eigenvalues of A can be bounded in an ellipse that is entirely to the right of the origin.

There is also dead code for imitating ML's Chebyshev implementation (ML_Cheby(), in packages/ml/src/Smoother/ml_smoother.c). I couldn't get it to converge in time to be useful for testing, so it is disabled.

Member Typedef Documentation

template<class ScalarType, class MV>

typedef ScalarType Ifpack2::Details::Chebyshev< ScalarType, MV >::ST

The type of entries in the matrix and vectors.

template<class ScalarType, class MV>

typedef Teuchos::ScalarTraits<ScalarType> Ifpack2::Details::Chebyshev< ScalarType, MV >::STS

Traits class for ST.

template<class ScalarType, class MV>

typedef STS::magnitudeType Ifpack2::Details::Chebyshev< ScalarType, MV >::MT

The type of the absolute value of a ScalarType.

template<class ScalarType, class MV>

typedef Tpetra::Operator<typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type> Ifpack2::Details::Chebyshev< ScalarType, MV >::op_type

Specialization of Tpetra::Operator.

template<class ScalarType, class MV>

typedef Tpetra::RowMatrix<typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type> Ifpack2::Details::Chebyshev< ScalarType, MV >::row_matrix_type

Specialization of Tpetra::RowMatrix.

template<class ScalarType, class MV>

typedef Tpetra::Vector<typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type> Ifpack2::Details::Chebyshev< ScalarType, MV >::V

Specialization of Tpetra::Vector.

template<class ScalarType, class MV>

typedef Tpetra::Map<typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type> Ifpack2::Details::Chebyshev< ScalarType, MV >::map_type

Specialization of Tpetra::Map.

Constructor & Destructor Documentation

template<class ScalarType , class MV >

Ifpack2::Details::Chebyshev< ScalarType, MV >::Chebyshev ( Teuchos::RCP< const row_matrix_type > A )

Constructor that takes a sparse matrix and sets default parameters.

Parameters

A	[in] The matrix A in the linear system to solve. If A is nonnull, it must be real-valued and symmetric positive definite. The input A may be null. In that case, you must call setMatrix() with a nonnull input before you may call compute() or apply().

template<class ScalarType , class MV >

Ifpack2::Details::Chebyshev< ScalarType, MV >::Chebyshev	(	Teuchos::RCP< const row_matrix_type >	A,
		Teuchos::ParameterList &	params
	)

Constructor that takes a sparse matrix and sets the user's parameters.

Parameters

A	[in] The matrix A in the linear system to solve. If A is nonnull, it must be real-valued and symmetric positive definite. The input A may be null. In that case, you must call setMatrix() with a nonnull input before you may call compute() or apply().
params	[in/out] On input: the parameters. On output: filled with the current parameter settings.

Member Function Documentation

template<class ScalarType , class MV >

void Ifpack2::Details::Chebyshev< ScalarType, MV >::setParameters ( Teuchos::ParameterList & plist )

Set (or reset) parameters.

This method fills in the input ParameterList with missing parameters set to their default values. You may call this method as many times as you want. On each call, the input ParameterList is treated as a complete list of the desired parameters, not as a "delta" or change list from the current set of parameters. (That is, if you remove parameters from the list that were there in the last call to setParameters() and call setParameters() again with the revised list, this method will use default values for the removed parameters, rather than letting the current settings remain.) However, since the method fills in missing parameters, you may keep calling it with the ParameterList used in the previous call in order to get the same behavior as before.

List of parameters

Parameters that govern spectral bounds of the matrix:

"chebyshev: max eigenvalue" (ScalarType): lambdaMax, an upper bound of the bounding ellipse of the eigenvalues of the matrix A. If you do not set this parameter, we will compute an approximation. See "Parameters that govern eigenvalue analysis" to control this approximation process.
"chebyshev: ratio eigenvalue" (ScalarType): eigRatio, the ratio of lambdaMax to the lower bound of the bounding ellipse of the eigenvalues of A. We use lambdaMax and eigRatio to determine the Chebyshev iteration coefficients. This parameter is optional and defaults to 30.
"chebyshev: min eigenvalue" (ScalarType): lambdaMin, a lower bound of real part of bounding ellipse of eigenvalues of the matrix A. This parameter is optional and only used for a quick check if the matrix is the identity matrix (if lambdaMax == lambdaMin == 1).

Parameters that govern the number of Chebyshev iterations:

"chebyshev: degree" (int): numIters, the number of iterations. This overrides "relaxation: sweeps" and "smoother: sweeps" (see below).
"relaxation: sweeps" (int): numIters, the number of iterations. We include this for compatibility with Ifpack. This overrides "smoother: sweeps" (see below).
"smoother: sweeps" (int): numIters, as above. We include this for compatibility with ML.

Parameters that govern eigenvalue analysis:

"chebyshev: eigenvalue max iterations" (int): eigMaxIters, the number of power method iterations used to compute the maximum eigenvalue. This overrides "eigen-analysis: iterations" (see below).
"eigen-analysis: iterations" (int): eigMaxIters, as above. We include this parameter for compatibility with ML.
"eigen-analysis: type" (std::string): The algorithm to use for estimating the max eigenvalue. This parameter is optional. Currently, we only support "power-method" (or "power method"), which is what Ifpack::Chebyshev uses for eigenanalysis. We include this parameter for compatibility with ML.

Parameters that govern other algorithmic details:

"chebyshev: assume matrix does not change": Whether compute() should always assume that the matrix has not changed since the last call to compute(). The default is false. If true, compute() will not recompute the inverse diagonal or the estimates of the max and min eigenvalues. compute() will always compute any quantity which the user did not provide and which we have not yet computed before.
"chebyshev: operator inv diagonal" (RCP<const V> or const V*): If nonnull, we will use a deep copy of this vector for left scaling as the inverse diagonal of the matrix A, instead of computing the inverse diagonal ourselves. We will make a copy every time you call setParameters(). If you ever call setParameters() without this parameter, we will clear our copy and compute the inverse diagonal ourselves again. If you choose to provide this parameter, you are responsible for updating this if the matrix has changed.
"chebyshev: min diagonal value" (ST): minDiagVal. If any entry of the diagonal of the matrix is less than this in magnitude, it will be replaced with this value in the inverse diagonal used for left scaling.
"chebyshev: zero starting solution" (bool): If true, then always use the zero vector(s) as the initial guess(es). If false, then apply() will use X on input as the initial guess(es).
"chebyshev: compute spectral radius" (bool): If true, the power method will compute the spectral radius of the operator. If false, it will compute the dominant eigenvalue.

Parameters that govern backwards compatibility:

"chebyshev: textbook algorithm" (bool): If true, use the textbook version of Chebyshev iteration. We recommend against this, since the default algorithm is less sensitive to the quality of the eigenvalue bounds.
"chebyshev: compute max residual norm" (bool): If true, apply() will compute and return the max (absolute) residual norm. Otherwise, apply() returns 0. This defaults to false.

The above compatibility parameters are not exposed in the documentation of Ifpack2::Chebyshev, because they are more useful to Ifpack2 developers than to users.

Precondition: lambdaMin, lambdaMax, and eigRatio are real; 0 < lambdaMin <= lambdaMax; numIters >= 0; eigMaxIters >= 0

Default settings for parameters relating to spectral bounds come from Ifpack.

template<class ScalarType , class MV >

void Ifpack2::Details::Chebyshev< ScalarType, MV >::compute ( )

(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A.

You must call this method before calling apply(),

if you have not yet called this method,
if the matrix (either its values or its structure) has changed, or
any time after you call setParameters().

The input matrix must be nonnull before you may call this method. If the input matrix is null, you must first call setMatrix() with a nonnull input matrix before you may call this method.

Users have the option to supply the left scaling vector D_inv and estimates of the min and max eigenvalues of D_inv * A as parameters to setParameters(). If users did not supply a left scaling, then this method will compute it by default (if assumeMatrixUnchanged is false). Likewise, if users did not supply at least an estimate of the max eigenvalue, this method will estimate it by default. If estimation of the eigenvalues is required, this method may take as long as several Chebyshev iterations.

Advanced users may avoid recomputing the left scaling vector and eigenvalue estimates by setting the "chebyshev: assume matrix does not change" parameter of setParameters() to true. The left scaling vector and eigenvalue estimates will always be computed if the user did not provide them and we have not yet computed them. Any changes to parameters that affect computation of the inverse diagonal or estimation of the eigenvalue bounds will not affect subsequent apply() operations, until the "chebyshev: assume matrix does not change" parameter is set back to false (its default value).

template<class ScalarType , class MV>

Chebyshev< ScalarType, MV >::MT Ifpack2::Details::Chebyshev< ScalarType, MV >::apply	(	const MV &	B,
		MV &	X
	)

Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.

Parameters

B	[in] Right-hand side(s) in the linear system to solve.
X	[in] Initial guess(es) for the linear system to solve.

If the "chebyshev: compute max residual norm" parameter is true (not the default), then this method returns the maximum (over all columns) absolute residual 2-norm after iterating. Otherwise, it returns zero.

Warning: If you did not set the "chebyshev: zero starting solution" parameter to true, then this method will use X as the starting guess for Chebyshev iteration. If you did not initialize X before calling this method, then the resulting solution will be undefined, since it will be computed using uninitialized data.

template<class ScalarType , class MV >

Teuchos::RCP< const typename Chebyshev< ScalarType, MV >::row_matrix_type > Ifpack2::Details::Chebyshev< ScalarType, MV >::getMatrix ( ) const

Get the matrix given to the constructor.

template<class ScalarType , class MV >

void Ifpack2::Details::Chebyshev< ScalarType, MV >::setMatrix ( const Teuchos::RCP< const row_matrix_type > & A )

Set the matrix.

It's legal to call this method with a null input. In that case, one must then call this method with a nonnull input before one may call compute() or apply().

template<class ScalarType , class MV >

bool Ifpack2::Details::Chebyshev< ScalarType, MV >::hasTransposeApply ( ) const

Whether it's possible to apply the transpose of this operator.

template<class ScalarType , class MV >

void Ifpack2::Details::Chebyshev< ScalarType, MV >::print ( std::ostream & out )

Print instance data to the given output stream.

template<class ScalarType , class MV >

std::string Ifpack2::Details::Chebyshev< ScalarType, MV >::description ( ) const

virtual

A single-line description of the Chebyshev solver.

Reimplemented from Teuchos::Describable.

template<class ScalarType , class MV >

void Ifpack2::Details::Chebyshev< ScalarType, MV >::describe	(	Teuchos::FancyOStream &	out,
		const Teuchos::EVerbosityLevel	verbLevel = `Teuchos::Describable::verbLevel_default`
	)		const

virtual

Print a description of the Chebyshev solver to out.

Reimplemented from Teuchos::Describable.

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

Public Types

Public Member Functions

Detailed Description

template<class ScalarType, class MV> class Ifpack2::Details::Chebyshev< ScalarType, MV >

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

List of parameters

template<class ScalarType, class MV>
class Ifpack2::Details::Chebyshev< ScalarType, MV >