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

Private Types
typedef Teuchos::ScalarTraits < scalar_type >	STS

typedef Teuchos::ScalarTraits < magnitude_type >	STM

typedef Teuchos::BLAS< int, scalar_type >	blas_type

typedef Teuchos::LAPACK< int, scalar_type >	lapack_type

Private Member Functions
int	solveLeastSquaresUsingSVD (mat_type &A, mat_type &X)
	Solve $\min_X \\|AX - B\\|_2$ using the SVD. More...

void	solveGivens (mat_type &y, mat_type &R, const mat_type &z, const int curCol)
	Solve the projected least-squares problem, assuming Givens rotations updates. More...

void	makeRandomProblem (mat_type &H, mat_type &z)
	Make a random projected least-squares problem. More...

void	computeGivensRotation (const Scalar &x, const Scalar &y, Scalar &theCosine, Scalar &theSine, Scalar &result)
	Compute the Givens rotation corresponding to [x; y]. More...

void	singularValues (const mat_type &A, Teuchos::ArrayView< magnitude_type > sigmas)
	Compute the singular values of A. Store them in the sigmas array. More...

std::pair< magnitude_type, magnitude_type >	extremeSingularValues (const mat_type &A)
	The (largest, smallest) singular values of the given matrix. More...

magnitude_type	leastSquaresConditionNumber (const mat_type &A, const mat_type &b, const magnitude_type &residualNorm)
	Normwise 2-norm condition number of the least-squares problem. More...

magnitude_type	leastSquaresResidualNorm (const mat_type &A, const mat_type &x, const mat_type &b)
	$\\| b - A x \\|_2$ (Frobenius norm if b has more than one column). More...

magnitude_type	solutionError (const mat_type &x_approx, const mat_type &x_exact)
	\|\|x_approx - x_exact\|\|_2 // \|\|x_exact\|\|_2. More...

magnitude_type	updateColumnGivens (const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int curCol)
	Update current column using Givens rotations. More...

magnitude_type	solveLapack (const mat_type &H, mat_type &R, mat_type &y, mat_type &z, const int curCol)
	Solve the least-squares problem using LAPACK. More...

magnitude_type	updateColumnsGivens (const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
	Update columns [startCol,endCol] of the projected least-squares problem. More...

magnitude_type	updateColumnsGivensBlock (const mat_type &H, mat_type &R, mat_type &y, mat_type &z, Teuchos::ArrayView< scalar_type > theCosines, Teuchos::ArrayView< scalar_type > theSines, const int startCol, const int endCol)
	Update columns [startCol,endCol] of the projected least-squares problem. More...

Private Attributes
std::ostream &	warn_
	Stream to which to output warnings. More...

ERobustness	defaultRobustness_
	Default robustness level, for things like 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.

template<class Scalar >

typedef Teuchos::ScalarTraits<scalar_type> Belos::details::ProjectedLeastSquaresSolver< Scalar >::STS

private

Definition at line 950 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

typedef Teuchos::ScalarTraits<magnitude_type> Belos::details::ProjectedLeastSquaresSolver< Scalar >::STM

private

Definition at line 951 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

typedef Teuchos::BLAS<int, scalar_type> Belos::details::ProjectedLeastSquaresSolver< Scalar >::blas_type

private

Definition at line 952 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

typedef Teuchos::LAPACK<int, scalar_type> Belos::details::ProjectedLeastSquaresSolver< Scalar >::lapack_type

private

Definition at line 953 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.

template<class Scalar >

int Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveLeastSquaresUsingSVD	(	mat_type &	A,
		mat_type &	X
	)

inlineprivate

Solve $\min_X \|AX - B\|_2$ using the SVD.

Parameters

X	[in/out] On input: the right-hand side(s) B. On output: the solution to the least-squares problem.
A	[in/out] On input: the matrix in the least-squares problem. On output, A will be overwritten with factorization information.

Returns: Numerical rank of A, computed using Scalar's machine precision as the rank tolerance.

Definition at line 1705 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

void Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveGivens	(	mat_type &	y,
		mat_type &	R,
		const mat_type &	z,
		const int	curCol
	)

inlineprivate

Solve the projected least-squares problem, assuming Givens rotations updates.

Call this method after invoking either updateColumnGivens() with the same curCol, or updateColumnsGivens() with curCol = endCol.

Definition at line 1780 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

void Belos::details::ProjectedLeastSquaresSolver< Scalar >::makeRandomProblem	(	mat_type &	H,
		mat_type &	z
	)

inlineprivate

Make a random projected least-squares problem.

Definition at line 1800 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

void Belos::details::ProjectedLeastSquaresSolver< Scalar >::computeGivensRotation	(	const Scalar &	x,
		const Scalar &	y,
		Scalar &	theCosine,
		Scalar &	theSine,
		Scalar &	result
	)

inlineprivate

Compute the Givens rotation corresponding to [x; y].

The result of applying the rotation is [result; 0].

Definition at line 1845 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

void Belos::details::ProjectedLeastSquaresSolver< Scalar >::singularValues	(	const mat_type &	A,
		Teuchos::ArrayView< magnitude_type >	sigmas
	)

inlineprivate

Compute the singular values of A. Store them in the sigmas array.

Definition at line 1874 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

std::pair<magnitude_type, magnitude_type> Belos::details::ProjectedLeastSquaresSolver< Scalar >::extremeSingularValues ( const mat_type & A )

inlineprivate

The (largest, smallest) singular values of the given matrix.

We use these for computing the 2-norm condition number of the matrix A. We separate out the singular values rather than returning their quotient, so that you can see the value of the largest singular value, even if the smallest singular value is zero.

Definition at line 1934 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::leastSquaresConditionNumber	(	const mat_type &	A,
		const mat_type &	b,
		const magnitude_type &	residualNorm
	)

inlineprivate

Normwise 2-norm condition number of the least-squares problem.

Parameters

A	[in] The matrix in the least-squares problem.
b	[in] The right-hand side in the least-squares problem.
residualNorm	[in] Residual norm from solving the least-squares problem using a known good method.

For details on the condition number formula, see Section 3.3 of J. W. Demmel, "Applied Numerical Linear Algebra," SIAM Press.

Definition at line 1957 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::leastSquaresResidualNorm	(	const mat_type &	A,
		const mat_type &	x,
		const mat_type &	b
	)

inlineprivate

$\| b - A x \|_2$ (Frobenius norm if b has more than one column).

Definition at line 2005 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::solutionError	(	const mat_type &	x_approx,
		const mat_type &	x_exact
	)

inlineprivate

||x_approx - x_exact||_2 // ||x_exact||_2.

Use the Frobenius norm if more than one column. Don't scale if ||x_exact|| == 0.

Definition at line 2023 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::updateColumnGivens	(	const mat_type &	H,
		mat_type &	R,
		mat_type &	y,
		mat_type &	z,
		Teuchos::ArrayView< scalar_type >	theCosines,
		Teuchos::ArrayView< scalar_type >	theSines,
		const int	curCol
	)

inlineprivate

Update current column using Givens rotations.

Update the current column of the projected least-squares problem, using Givens rotations. This updates the QR factorization of the upper Hessenberg matrix H. The resulting R factor is stored in the matrix R. The Q factor is stored implicitly in the list of cosines and sines, representing the Givens rotations applied to the problem. These Givens rotations are also applied to the right-hand side z.

Returns: The residual of the resulting least-squares problem, assuming that the upper triangular system Ry=z can be solved exactly (with zero residual). (This may not be the case if R is singular and the system Ry=z is inconsistent.)

Parameters

H	[in] The upper Hessenberg matrix from GMRES. We only view H( 1:curCol+2, 1:curCol+1 ). It's copied and not overwritten, so as not to disturb any backtracking or other features of GMRES.
R	[in/out] Matrix with the same dimensions as H, used for storing the incrementally computed upper triangular factor from the QR factorization of H.
y	[out] Vector (one-column matrix) with the same number of rows as H. On output: the solution of the projected least-squares problem. The vector should have one more entry than necessary for the solution, because of the way we solve the least-squares problem.
z	[in/out] Vector (one-column matrix) with the same number of rows as H. On input: the current right-hand side of the projected least-squares problem. On output: the updated right-hand side.
theCosines	[in/out] On input: Array of cosines from the previously computed Givens rotations. On output: the same cosines, with the new Givens rotation's cosine appended.
theSines	[in/out] On input: Array of sines from the previously computed Givens rotations. On output: the same sines, with the new Givens rotation's sine appended.
curCol	[in] Index of the current column to update.

Definition at line 2089 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::solveLapack	(	const mat_type &	H,
		mat_type &	R,
		mat_type &	y,
		mat_type &	z,
		const int	curCol
	)

inlineprivate

Solve the least-squares problem using LAPACK.

In particular, this method uses LAPACK's _GELS dense least-squares solver. _GELS assumes that the input matrix has full rank, so this method might fail. This method is inefficient compared to the Givens rotations approach, but useful for testing the latter.

Parameters

H	[in] The upper Hessenberg matrix from GMRES. We only view H( 1:curCol+2, 1:curCol+1 ). It's copied and not overwritten, so as not to disturb any backtracking or other features of GMRES.
R	[out] Matrix with the same dimensions as H. On output: the upper triangle of R contains the R factor in the QR factorization of H.
y	[out] Vector (one-column matrix) with the same number of rows as H. On output: the solution of the projected least-squares problem. The vector should have one more entry than necessary for the solution, because of the way we solve the least-squares problem.
z	[in] Vector (one-column matrix) with the same number of rows as H. On input: the current right-hand side of the projected least-squares problem.
curCol	[in] Index of the current column of H.

Returns: The 2-norm of the residual of the least-squares problem, assuming that it can be solved using LAPACK's _GELS least-squares solver. (This may not be the case if H is singular.)

Definition at line 2226 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::updateColumnsGivens	(	const mat_type &	H,
		mat_type &	R,
		mat_type &	y,
		mat_type &	z,
		Teuchos::ArrayView< scalar_type >	theCosines,
		Teuchos::ArrayView< scalar_type >	theSines,
		const int	startCol,
		const int	endCol
	)

inlineprivate

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

This method implements a "left-looking panel QR factorization" of the upper Hessenberg matrix in the projected least-squares problem. It's "left-looking" because we don't updating anything to the right of columns [startCol, endCol], which is the "panel."

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

Definition at line 2314 of file BelosProjectedLeastSquaresSolver.hpp.

template<class Scalar >

magnitude_type Belos::details::ProjectedLeastSquaresSolver< Scalar >::updateColumnsGivensBlock	(	const mat_type &	H,
		mat_type &	R,
		mat_type &	y,
		mat_type &	z,
		Teuchos::ArrayView< scalar_type >	theCosines,
		Teuchos::ArrayView< scalar_type >	theSines,
		const int	startCol,
		const int	endCol
	)

inlineprivate

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

This method implements a "left-looking panel QR factorization" of the upper Hessenberg matrix in the projected least-squares problem. It's "left-looking" because we don't updating anything to the right of columns [startCol, endCol], which is the "panel."