Methods for solving GMRES' projected least-squares problem. More...

#include <BelosProjectedLeastSquaresSolver.hpp>

Public Types
typedef Scalar	scalar_type
	The template parameter of this class. More...

typedef Teuchos::ScalarTraits < Scalar >::magnitudeType	magnitude_type
	The type of the magnitude of a `scalar_type` value. More...

typedef Teuchos::SerialDenseMatrix < int, Scalar >	mat_type
	The type of a dense matrix (or vector) of `scalar_type`. More...

Public Member Functions
	ProjectedLeastSquaresSolver (std::ostream &warnStream, const ERobustness defaultRobustness=ROBUSTNESS_NONE)
	Constructor. More...

magnitude_type	updateColumn (ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
	Update column curCol of the projected least-squares problem. More...

magnitude_type	updateColumns (ProjectedLeastSquaresProblem< Scalar > &problem, const int startCol, const int endCol)
	Update columns [startCol,endCol] of the projected least-squares problem. More...

void	solve (ProjectedLeastSquaresProblem< Scalar > &problem, const int curCol)
	Solve the projected least-squares problem. More...

std::pair< int, bool >	solveUpperTriangularSystem (Teuchos::ESide side, mat_type &X, const mat_type &R, const mat_type &B)
	Solve the given square upper triangular linear system(s). More...

std::pair< int, bool >	solveUpperTriangularSystem (Teuchos::ESide side, mat_type &X, const mat_type &R, const mat_type &B, const ERobustness robustness)
	Solve the given square upper triangular linear system(s). More...

std::pair< int, bool >	solveUpperTriangularSystemInPlace (Teuchos::ESide side, mat_type &X, const mat_type &R)
	Solve square upper triangular linear system(s) in place. More...

std::pair< int, bool >	solveUpperTriangularSystemInPlace (Teuchos::ESide side, mat_type &X, const mat_type &R, const ERobustness robustness)
	Solve square upper triangular linear system(s) in place. More...

bool	testGivensRotations (std::ostream &out)
	Test Givens rotations. More...

bool	testUpdateColumn (std::ostream &out, const int numCols, const bool testBlockGivens=false, const bool extraVerbose=false)
	Test update and solve using Givens rotations. More...

bool	testTriangularSolves (std::ostream &out, const int testProblemSize, const ERobustness robustness, const bool verbose=false)
	Test upper triangular solves. More...

Detailed Description

template<class Scalar>
class Belos::details::ProjectedLeastSquaresSolver< Scalar >

Methods for solving GMRES' projected least-squares problem.

Author: Mark Hoemmen

Template Parameters

Scalar The type of the matrix and vector entries in the least-squares problem.

Expected use of this class:

Use a ProjectedLeastSquaresProblem struct instance to store the projected problem in your GMRES solver.
Instantiate a ProjectedLeastSquaresSolver:
ProjectedLeastSquaresSolver<Scalar> solver;
Update the current column(s) of the QR factorization of GMRES' upper Hessenberg matrix:
typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;

// Standard GMRES: update one column at a time.

MT resNorm = solver.updateColumn (problem, endCol);

or
// Some GMRES variants update multiple columns at a time.

// You can also use this feature to recompute the upper

// Hessenberg's QR factorization from scratch.

MT resNorm = solver.updateColumns (problem, startCol, endCol);
Solve for the current GMRES solution update coefficients:
solver.solve (problem, endCol);

You can defer Step 4 as long as you want. Step 4 must always follow Step 3.

Purposes of this class:

Isolate and factor out BLAS and LAPACK dependencies. This makes it easier to write custom replacements for routines for which no Scalar specialization is available.
Encapsulate common functionality of many GMRES-like solvers. Avoid duplicated code and simplify debugging, testing, and implementation of new GMRES-like solvers.
Provide an option for more robust implementations of solvers for the projected least-squares problem.

"Robust" here means regularizing the least-squares solve, so that the solution is well-defined even if the problem is ill-conditioned. Many distributed-memory iterative solvers, including those in Belos, currently solve the projected least-squares problem redundantly on different processes. If those processes are heterogeneous or implement the BLAS and LAPACK themselves in parallel (via multithreading, for example), then different BLAS or LAPACK calls on different processes may result in different answers. The answers may be significantly different if the projected problem is singular or ill-conditioned. This is bad because GMRES variants use the projected problem's solution as the coefficients for the solution update. The solution update coefficients must be (almost nearly) the same on all processes. Regularizing the projected problem is one way to ensure that different processes compute (almost) the same solution.

Definition at line 931 of file BelosProjectedLeastSquaresSolver.hpp.

Member Typedef Documentation

template<class Scalar >

Belos::details::ProjectedLeastSquaresSolver< Scalar >::scalar_type

The template parameter of this class.

The type of the matrix and vector entries in the projected least-squares problem to solve, and the type of the resulting GMRES solution update coefficients.

Definition at line 939 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

Belos::details::ProjectedLeastSquaresSolver< Scalar >::magnitude_type

The type of the magnitude of a scalar_type value.

If scalar_type is complex-valued, then magnitude_type is real.

Definition at line 944 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

Belos::details::ProjectedLeastSquaresSolver< Scalar >::mat_type

The type of a dense matrix (or vector) of scalar_type.

Definition at line 947 of file BelosProjectedLeastSquaresSolver.hpp.

Constructor & Destructor Documentation

template<class Scalar >

Belos::details::ProjectedLeastSquaresSolver< Scalar >::ProjectedLeastSquaresSolver	(	std::ostream &	warnStream,
		const ERobustness	defaultRobustness = `ROBUSTNESS_NONE`
	)

inline

Constructor.

Parameters

warnStream	[out] Stream to which to output warnings. Set to a Teuchos::oblackholestream if you don't want to display warnings.
defaultRobustness	[in] Default robustness level for operations like triangular solves. For example, at a low robustness level, triangular solves might fail if the triangular matrix is singular or ill-conditioned in a particular way. At a high robustness level, triangular solves will always succeed, but may only be solved in a least-squares sense.

Definition at line 969 of file BelosProjectedLeastSquaresSolver.hpp.

Member Function Documentation

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::updateColumn	(	ProjectedLeastSquaresProblem< Scalar > &	problem,
		const int	curCol
	)

inline

Update column curCol of the projected least-squares problem.

The upper Hessenberg matrix H is read but not touched. The R factor, the cosines and sines, and the right-hand side z are updated. This method does not compute the solution of the least-squares problem; call solve() for that.

Parameters

problem	[in/out] The projected least-squares problem.
curCol	[in] Zero-based index of the current column to update.

Returns: 2-norm of the absolute residual of the projected least-squares problem.

Definition at line 989 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::updateColumns	(	ProjectedLeastSquaresProblem< Scalar > &	problem,
		const int	startCol,
		const int	endCol
	)

inline

Update columns [startCol,endCol] of the projected least-squares problem.

The upper Hessenberg matrix H is read but not touched. The R factor, the cosines and sines, and the right-hand side z are updated. This method does not compute the solution of the least-squares problem; call solve() for that.

Parameters

problem	[in/out] The projected least-squares problem.
startCol	[in] Zero-based index of the first column to update.
endCol	[in] Zero-based index of the last column (inclusive) to update.

Returns: 2-norm of the absolute residual of the projected least-squares problem.

Definition at line 1011 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

void Belos::details::ProjectedLeastSquaresSolver< Scalar >::solve	(	ProjectedLeastSquaresProblem< Scalar > &	problem,
		const int	curCol
	)

inline

Solve the projected least-squares problem.

Call this method only after calling updateColumn() or updateColumns(). If you call updateColumn(), use the same column index when calling this method as you did when calling updateColumn(). If you call updateColumns(), use its endCol argument as the column index for calling this method.

Parameters

problem	[in/out] The projected least-squares problem.
curCol	[in] Zero-based index of the most recently updated column of the least-squares problem.

Definition at line 1034 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

std::pair<int, bool> Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveUpperTriangularSystem	(	Teuchos::ESide	side,
		mat_type &	X,
		const mat_type &	R,
		const mat_type &	B
	)

inline

Solve the given square upper triangular linear system(s).

This method does the same thing as the five-argument overload, but uses the default robustness level.

Definition at line 1045 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

std::pair<int, bool> Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveUpperTriangularSystem	(	Teuchos::ESide	side,
		mat_type &	X,
		const mat_type &	R,
		const mat_type &	B,
		const ERobustness	robustness
	)

inline

Solve the given square upper triangular linear system(s).

Parameters

side	[in] Whether to solve with the triangular matrix R on the left side or right side of the matrix X of unknowns. side = Teuchos::LEFT_SIDE or Teuchos::RIGHT_SIDE.
X	[out] The matrix of unknowns. (We consider a matrix with one column a vector.)
R	[in] The square upper triangular matrix. If R has more rows than columns, this method uses the number of columns of R as the dimension of the square matrix, and ignores the "extra" rows of R. If R has more columns than rows, this method throws an exception.
B	[in] The matrix of right-hand sides. B must have at least as many columns as X, otherwise this method will throw an exception. If B has more columns than X, those columns will be ignored.
robustness	[in] The robustness level to use for the solve. ROBUSTNESS_NONE is fastest but does not attempt to detect rank deficiency in R. ROBUSTNESS_LOTS is slowest, but can "solve" the problem (in a least-squares minimum-norm-solution sense) even if R is singular. ROBUSTNESS_SOME may do something in between.

Depending on the specified robustness level, this method may also attempt to compute the numerical rank and detect rank deficiency in R. The first return value is the computed numerical rank, and the second return value is a Boolean indicating whether any rank deficiency was detected. If the method does not attempt to compute the numerical rank or detect rank deficiency, it will return (N, false), where N is the number of columns in R.

Returns: (detectedRank, foundRankDeficiency).

Definition at line 1092 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

std::pair<int, bool> Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveUpperTriangularSystemInPlace	(	Teuchos::ESide	side,
		mat_type &	X,
		const mat_type &	R
	)

inline

Solve square upper triangular linear system(s) in place.

This method does the same thing as its four-argument overload, but uses the default robustness level.

Definition at line 1128 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

std::pair<int, bool> Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveUpperTriangularSystemInPlace	(	Teuchos::ESide	side,
		mat_type &	X,
		const mat_type &	R,
		const ERobustness	robustness
	)

inline

Solve square upper triangular linear system(s) in place.

This is the "in-place" version of solveUpperTriangularSystem(). The difference between that method and this one is that this method assumes that the right-hand side(s) is/are already stored in X on input. It then overwrites X with the solution.

Definition at line 1143 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

bool Belos::details::ProjectedLeastSquaresSolver< Scalar >::testGivensRotations ( std::ostream & out )

inline

Test Givens rotations.

This routine tests both computing Givens rotations (via computeGivensRotation()) and applying them.

Parameters

out	[out] Stream to which to write test output.

Returns: true if the test succeeded, else false.

Definition at line 1294 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

bool Belos::details::ProjectedLeastSquaresSolver< Scalar >::testUpdateColumn	(	std::ostream &	out,
		const int	numCols,
		const bool	testBlockGivens = `false`,
		const bool	extraVerbose = `false`
	)

inline

Test update and solve using Givens rotations.

Parameters

out	[out] Output stream to which to print results.
numCols	[in] Number of columns in the projected least-squares problem to test. (The number of rows is one plus the number of columns.)
testBlockGivens	[in] Whether to test the "block" (i.e., panel) version of the Givens rotations update.
extraVerbose	[in] Whether to print extra-verbose output (e.g., the test problem and results).

Returns: Whether the test succeeded, meaning that none of the solves reported failure and the least-squares solution error was within the expected bound.

Test updating and solving the least-squares problem using Givens rotations by comparison against LAPACK's least-squares solver. First generate a random least-squares problem that looks like it comes from GMRES. The matrix is upper Hessenberg, and the right-hand side starts out with the first entry being nonzero with nonnegative real part and zero imaginary part, and all the other entries being zero. Then compare the results of updateColumnGivens() (applied to each column in turn) and solveGivens() (applied at the end) with the results of solveLapack() (applied at the end). Print out the results to the given output stream.

Definition at line 1353 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

bool Belos::details::ProjectedLeastSquaresSolver< Scalar >::testTriangularSolves	(	std::ostream &	out,
		const int	testProblemSize,
		const ERobustness	robustness,
		const bool	verbose = `false`
	)

inline

Test upper triangular solves.

Parameters

out	[out] Stream to which this method writes output.
testProblemSize	[in] Dimension (number of rows and columns) of the upper triangular matrices tested by this method.
robustness	[in] Robustness level of the upper triangular solves.
verbose	[in] Whether to print more verbose output.

Returns: True if test passed, else false.

Warning: Before calling this routine, seed the random number generator via Teuchos::ScalarTraits<Scalar>::seedrandom (seed).

Definition at line 1594 of file BelosProjectedLeastSquaresSolver.hpp.

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

BelosProjectedLeastSquaresSolver.hpp

Public Types

Public Member Functions

Detailed Description

template<class Scalar> class Belos::details::ProjectedLeastSquaresSolver< Scalar >

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

template<class Scalar>
class Belos::details::ProjectedLeastSquaresSolver< Scalar >