Bordered system solver strategy based on Householder transformations. More...

#include <LOCA_BorderedSolver_TpetraHouseholder.hpp>

Inheritance diagram for LOCA::BorderedSolver::TpetraHouseholder:

Collaboration diagram for LOCA::BorderedSolver::TpetraHouseholder:

Public Member Functions
	TpetraHouseholder (const Teuchos::RCP< LOCA::GlobalData > &global_data, const Teuchos::RCP< LOCA::Parameter::SublistParser > &topParams, const Teuchos::RCP< Teuchos::ParameterList > &solverParams)
	Constructor. More...

virtual	~TpetraHouseholder ()
	Destructor.

virtual void	setMatrixBlocks (const Teuchos::RCP< const LOCA::BorderedSolver::AbstractOperator > &op, const Teuchos::RCP< const NOX::Abstract::MultiVector > &blockA, const Teuchos::RCP< const LOCA::MultiContinuation::ConstraintInterface > &blockB, const Teuchos::RCP< const NOX::Abstract::MultiVector::DenseMatrix > &blockC)
	Set blocks. More...

virtual NOX::Abstract::Group::ReturnType	initForSolve ()
	Intialize solver for a solve. More...

virtual NOX::Abstract::Group::ReturnType	initForTransposeSolve ()
	Intialize solver for a transpose solve. More...

virtual NOX::Abstract::Group::ReturnType	apply (const NOX::Abstract::MultiVector &X, const NOX::Abstract::MultiVector::DenseMatrix &Y, NOX::Abstract::MultiVector &U, NOX::Abstract::MultiVector::DenseMatrix &V) const
	Computed extended matrix-multivector product. More...

virtual NOX::Abstract::Group::ReturnType	applyTranspose (const NOX::Abstract::MultiVector &X, const NOX::Abstract::MultiVector::DenseMatrix &Y, NOX::Abstract::MultiVector &U, NOX::Abstract::MultiVector::DenseMatrix &V) const
	Computed extended matrix transpose-multivector product. More...

virtual NOX::Abstract::Group::ReturnType	applyInverse (Teuchos::ParameterList &params, const NOX::Abstract::MultiVector F, const NOX::Abstract::MultiVector::DenseMatrix G, NOX::Abstract::MultiVector &X, NOX::Abstract::MultiVector::DenseMatrix &Y) const
	Solves the extended system using the technique described above. More...

virtual NOX::Abstract::Group::ReturnType	applyInverseTranspose (Teuchos::ParameterList &params, const NOX::Abstract::MultiVector F, const NOX::Abstract::MultiVector::DenseMatrix G, NOX::Abstract::MultiVector &X, NOX::Abstract::MultiVector::DenseMatrix &Y) const
	Solves the transpose of the extended system as defined above. More...

void	updateCrsMatrixForPreconditioner (const NOX::Abstract::MultiVector &U, const NOX::Abstract::MultiVector &V, NOX::TCrsMatrix &mat) const
	Overwrites the Jacobian with for computing the preconditioner of . More...

Public Member Functions inherited from LOCA::BorderedSolver::AbstractStrategy
	AbstractStrategy ()
	Constructor.

virtual	~AbstractStrategy ()
	Destructor.

virtual void	setMatrixBlocksMultiVecConstraint (const Teuchos::RCP< const LOCA::BorderedSolver::AbstractOperator > &op, const Teuchos::RCP< const NOX::Abstract::MultiVector > &blockA, const Teuchos::RCP< const NOX::Abstract::MultiVector > &blockB, const Teuchos::RCP< const NOX::Abstract::MultiVector::DenseMatrix > &blockC)
	Set blocks with multivector constraint. More...

Protected Types
enum	PRECONDITIONER_METHOD { JACOBIAN, SMW }
	Enumerated type indicating preconditioner method.

Protected Member Functions
virtual NOX::Abstract::Group::ReturnType	solve (Teuchos::ParameterList &params, const NOX::Abstract::MultiVector F, const NOX::Abstract::MultiVector::DenseMatrix G, NOX::Abstract::MultiVector &X, NOX::Abstract::MultiVector::DenseMatrix &Y) const
	Solves the extended system using the technique described above.

virtual NOX::Abstract::Group::ReturnType	solveTranspose (Teuchos::ParameterList &params, const NOX::Abstract::MultiVector F, const NOX::Abstract::MultiVector::DenseMatrix G, NOX::Abstract::MultiVector &X, NOX::Abstract::MultiVector::DenseMatrix &Y) const
	Solves the transpose of the extended system as defined above.

NOX::Abstract::Group::ReturnType	computeUV (const NOX::Abstract::MultiVector::DenseMatrix &Y1, const NOX::Abstract::MultiVector &Y2, const NOX::Abstract::MultiVector::DenseMatrix &T, const NOX::Abstract::MultiVector &A, NOX::Abstract::MultiVector &U, NOX::Abstract::MultiVector &V, bool use_jac_transpose)
	Compute and multivectors in .

Teuchos::RCP < NOX::Abstract::MultiVector >	createBlockMV (const NOX::Abstract::MultiVector &v) const

void	setBlockMV (const NOX::Abstract::MultiVector &bv, NOX::Abstract::MultiVector &v) const

Protected Attributes
Teuchos::RCP< LOCA::GlobalData >	globalData
	Global data object.

Teuchos::RCP < Teuchos::ParameterList >	solverParams
	Solver parameters.

Teuchos::RCP< LOCA::Thyra::Group >	grp
	Pointer to group storing J.

Teuchos::RCP< const LOCA::BorderedSolver::AbstractOperator >	op

Teuchos::RCP< const NOX::Abstract::MultiVector >	A
	Pointer to A block.

Teuchos::RCP< const NOX::Abstract::MultiVector >	B
	Pointer to B block.

Teuchos::RCP< const NOX::Abstract::MultiVector::DenseMatrix >	C
	Pointer to C block.

Teuchos::RCP< const LOCA::MultiContinuation::ConstraintInterfaceMVDX >	constraints
	Pointer to constraint interface.

LOCA::BorderedSolver::HouseholderQR	qrFact
	QR Factorization object.

Teuchos::RCP < NOX::Abstract::MultiVector >	house_x
	Solution component of Householder multivec.

NOX::Abstract::MultiVector::DenseMatrix	house_p
	Parameter component of Householder multivec.

NOX::Abstract::MultiVector::DenseMatrix	T
	T matrix in compact WY representation.

NOX::Abstract::MultiVector::DenseMatrix	R
	R matrix in QR factorization.

Teuchos::RCP < NOX::Abstract::MultiVector >	U
	U matrix in low-rank update form P = J + U*V^T.

Teuchos::RCP < NOX::Abstract::MultiVector >	V
	V matrix in low-rank update form P = J + U*V^T.

Teuchos::RCP < NOX::Abstract::MultiVector >	house_x_trans
	Solution component of Householder multivec for transposed system.

NOX::Abstract::MultiVector::DenseMatrix	house_p_trans
	Parameter component of Householder multivec for transposed system.

NOX::Abstract::MultiVector::DenseMatrix	T_trans
	T matrix in compact WY representation for transposed system.

NOX::Abstract::MultiVector::DenseMatrix	R_trans
	R matrix in QR factorization for transposed system.

Teuchos::RCP < NOX::Abstract::MultiVector >	U_trans
	U matrix in low-rank update form P = J + U*V^T for transposed system.

Teuchos::RCP < NOX::Abstract::MultiVector >	V_trans
	V matrix in low-rank update form P = J + U*V^T for transposed system.

Teuchos::RCP< const NOX::Abstract::MultiVector >	Ablock
	Pointer to A block as an Tpetra multivector.

Teuchos::RCP< const NOX::Abstract::MultiVector >	Bblock
	Pointer to B block as an Tpetra multivector.

Teuchos::RCP < NOX::Abstract::MultiVector >	Ascaled
	Pointer to scaled A block.

Teuchos::RCP < NOX::Abstract::MultiVector >	Bscaled
	Pointer to scaled B block.

Teuchos::RCP < NOX::Abstract::MultiVector::DenseMatrix >	Cscaled
	Pointer to scaled C block.

Teuchos::RCP< NOX::TOperator >	tpetraOp
	Pointer to Tpetra J operator.

Teuchos::RCP< NOX::TCrsMatrix >	tpetraPrecMatrix
	Pointer to Tpetra Preconditioner operator.

Teuchos::RCP <::Thyra::DefaultLinearOpSource < double > >	prec_losb
	Thyra wrapped preconditioner matrix (tpetraPrecMatrix) for when includeUV is true and use_P_for_Prec is false.

int	numConstraints
	Number of constraint equations.

bool	isZeroA
	flag indicating whether A block is zero

bool	isZeroB
	flag indicating whether B block is zero

bool	isZeroC
	flag indicating whether C block is zero

bool	isValidForSolve
	Flag indicating whether constraint factorization for solve has been computed.

bool	isValidForTransposeSolve
	Flag indicating whether constraint factorization for transpoe solve has been computed.

Teuchos::BLAS< int, double >	dblas
	BLAS Wrappers.

bool	scale_rows
	Whether we should scale augmented rows to have unit 2-norm.

std::vector< double >	scale_vals
	Scale values for each row.

PRECONDITIONER_METHOD	precMethod
	Preconditioner method.

bool	includeUV
	Flag indicating whether to include U*V^T terms in preconditioner.

bool	use_P_For_Prec
	Flag indicating whether to use P = J + U*V^T in preconditioner.

bool	isComplex
	Flag indicating whether we are doing a complex solve.

double	omega
	Frequency for complex systems.

Detailed Description

Bordered system solver strategy based on Householder transformations.

This class solves the extended system of equations

$\begin{bmatrix} J & A \\ B^T & C \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} F \\ G \end{bmatrix}$

using Householder tranformations. The algorithm works as follows: First consider a slightly rearranged version of the extended system of equations:

$\begin{bmatrix} C & B^T \\ A & J \end{bmatrix} \begin{bmatrix} Y \\ X \end{bmatrix} = \begin{bmatrix} G \\ F \end{bmatrix}.$

Let

$Q^T \begin{bmatrix} C^T \\ B \end{bmatrix} = \begin{bmatrix} R \\ 0 \end{bmatrix}$

be the QR decomposition of the constraints matrix where $Q\in\Re^{n+m\times n+m}$ and $R\in\Re^{m\times m}$ . Define

$\begin{bmatrix} Z_Y \\ Z_X \end{bmatrix} = Q^T \begin{bmatrix} Y \\ X \end{bmatrix},$

then the extended system of equations is equivalent to

$\begin{bmatrix} R^T & 0 \\ [A & J] Q \end{bmatrix} \begin{bmatrix} Z_Y \\ Z_X \end{bmatrix} = \begin{bmatrix} G \\ F \end{bmatrix}$

and hence

$\begin{split} Z_Y &= R^{-T} G \\ [A \;\; J] Q \begin{bmatrix} 0 \\ Z_X \end{bmatrix} &= F - [A \;\; J] Q \begin{bmatrix} Z_Y \\ 0 \end{bmatrix}. \end{split}$

This last equation equation can be written

$P Z_X = \tilde{F}$

where $P\in\Re^{n\times n}$ is given by

$P Z_X = [A \;\; J] Q \begin{bmatrix} 0 \\ Z_X \end{bmatrix}$

and

$\tilde{F} = F - [A \;\; J] Q \begin{bmatrix} Z_Y \\ 0 \end{bmatrix}.$

We then recover $X$ and $Y$ by

$\begin{bmatrix} Y \\ X \end{bmatrix} = Q \begin{bmatrix} Z_Y \\ Z_X \end{bmatrix}.$

It can be further shown that the $P$ operator above can be written

$P = J + U V^T$

where $U = A*Y_1 + J*Y_2$ , $V = Y_2*T^T$ and $Y = [Y_1 ; Y_2]$ . The equation $P Z_X = \tilde{F}$ is solved using an iterative solver using the definition of $P Z_X$ above, in this case AztecOO. The system is preconditioned using the preconditioner for $J$ . The operator $Q$ is generated using the standard Householder QR algorithm (Algorithm 5.2.1, G. Golub and C. Van Loan, "Matrix Computations," 3rd Edition, Johns Hopkins, Baltimore, 1996) and is stored using the compact WY representation: $Q = I + Y T Y^T$ (see R. Schreiber and C. Van Loan, "A Storage-Efficient WY Representation for Products of Householder Transformations," SIAM J. Sci. Stat. Comput., Vol. 10, No. 1, pp. 53-57, January 1989).

The operator representing $P$ is encapsulated in the class LOCA::Tpetra::LowRankUpdateRowMatrix if $J$ is an Tpetra::RowMatrix and LOCA::Tpetra::LowRankUpdateOp otherwise. If the row matrix version is available $P$ can be scaled and also used to construct a preconditioner. If "Include UV In Preconditioner" is true as discussed below, the $U$ and $V$ terms will be included when computing this preconditioner, which can help stability when $J$ is nearly singular.

The class is intialized via the solverParams parameter list argument to the constructor. The parameters this class recognizes are:

"Preconditioner Method" – [string] (default: "Jacobian") - Method for preconditioning the operator. Choices are:
- "Jacobian" (default) – Use the preconditioner for
- "SMW" – Use the Sherman-Morrison-Woodbury formula for the inverse of , replacing the inverse of with the preconditioner for .
"Scale Augmented Rows" – [bool] (default: true) - Scale augmented rows to unit 2-norm before computing QR factorization.
"Include UV In Preconditioner" – [bool] (default: false) - Flag indicating whether to use the and terms in the preconditioner for when using the "Jacobian" preconditioner method.
"Use P For Preconditioner" – [bool] (default: false) - Flag indicating whether to use the representation of as a LOCA::Tpetra::LowRankUpdateRowMatrix for computing the preconditioner when using the "Jacobian" preconditioner method. This is valid only for preconditioners that accept an Tpetra::RowMatrix interface.
"Transpose Solver Method" – [string] (default: "Transpose Preconditioner") Method for preconditioning the transpose linear system. See LOCA::Tpetra::TransposeLinearSystem::Factory for available choices.

Constructor & Destructor Documentation

LOCA::BorderedSolver::TpetraHouseholder::TpetraHouseholder	(	const Teuchos::RCP< LOCA::GlobalData > &	global_data,
		const Teuchos::RCP< LOCA::Parameter::SublistParser > &	topParams,
		const Teuchos::RCP< Teuchos::ParameterList > &	solverParams
	)

Constructor.

Parameters

global_data	[in] Global data object
topParams	[in] Parsed top-level parameter list
solverParams	[in] Bordered solver parameters as described above

References Teuchos::ParameterList::get(), globalData, includeUV, LOCA::GlobalData::locaErrorCheck, precMethod, rcp(), scale_rows, Teuchos::ParameterList::set(), solverParams, use_P_For_Prec, and Teuchos::ParameterList::validateParametersAndSetDefaults().

Member Function Documentation

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::apply	(	const NOX::Abstract::MultiVector &	X,
		const NOX::Abstract::MultiVector::DenseMatrix &	Y,
		NOX::Abstract::MultiVector &	U,
		NOX::Abstract::MultiVector::DenseMatrix &	V
	)		const

virtual

Computed extended matrix-multivector product.

Computes

$\begin{bmatrix} U \\ V \end{bmatrix} = \begin{bmatrix} J & A \\ B^T & C \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} J*X + A*Y \\ B^T*X + C*Y \end{bmatrix}.$

Implements LOCA::BorderedSolver::AbstractStrategy.

References NOX::Abstract::Group::Failed, Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::multiply(), Teuchos::NO_TRANS, and NOX::Abstract::MultiVector::update().

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::applyInverse	(	Teuchos::ParameterList &	params,
		const NOX::Abstract::MultiVector *	F,
		const NOX::Abstract::MultiVector::DenseMatrix *	G,
		NOX::Abstract::MultiVector &	X,
		NOX::Abstract::MultiVector::DenseMatrix &	Y
	)		const

virtual

Solves the extended system using the technique described above.

The params argument is the linear solver parameters. If isZeroF or isZeroG is true, than the corresponding F or G pointers may be NULL.

Note that if either the A or B blocks are zero, the system is solved using a simple block elimination scheme instead of the Householder scheme.

Implements LOCA::BorderedSolver::AbstractStrategy.

References Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::assign(), Teuchos::ParameterList::get(), Teuchos::RCP< T >::get(), NOX::Abstract::MultiVector::init(), Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::numCols(), Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::numRows(), NOX::Abstract::Group::Ok, Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::putScalar(), Teuchos::rcp(), LOCA::BorderedSolver::UpperTriangularBlockElimination::solve(), and LOCA::BorderedSolver::LowerTriangularBlockElimination::solve().

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::applyInverseTranspose	(	Teuchos::ParameterList &	params,
		const NOX::Abstract::MultiVector *	F,
		const NOX::Abstract::MultiVector::DenseMatrix *	G,
		NOX::Abstract::MultiVector &	X,
		NOX::Abstract::MultiVector::DenseMatrix &	Y
	)		const

virtual

Solves the transpose of the extended system as defined above.

The params argument is the linear solver parameters.

Implements LOCA::BorderedSolver::AbstractStrategy.

References Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::assign(), Teuchos::ParameterList::get(), Teuchos::RCP< T >::get(), NOX::Abstract::MultiVector::init(), Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::numCols(), Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::numRows(), NOX::Abstract::Group::Ok, Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::putScalar(), Teuchos::rcp(), LOCA::BorderedSolver::UpperTriangularBlockElimination::solveTranspose(), and LOCA::BorderedSolver::LowerTriangularBlockElimination::solveTranspose().

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::applyTranspose	(	const NOX::Abstract::MultiVector &	X,
		const NOX::Abstract::MultiVector::DenseMatrix &	Y,
		NOX::Abstract::MultiVector &	U,
		NOX::Abstract::MultiVector::DenseMatrix &	V
	)		const

virtual

Computed extended matrix transpose-multivector product.

Computes

$\begin{bmatrix} U \\ V \end{bmatrix} = \begin{bmatrix} J^T & B \\ A^T & C \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} J^T*X + B*Y \\ A^T*X + C^T*Y \end{bmatrix}.$

Implements LOCA::BorderedSolver::AbstractStrategy.

References NOX::Abstract::Group::Failed, Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::multiply(), NOX::Abstract::MultiVector::multiply(), Teuchos::NO_TRANS, and Teuchos::TRANS.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::initForSolve ( )

virtual

Intialize solver for a solve.

This should be called after setMatrixBlocks(), but before applyInverse().

Implements LOCA::BorderedSolver::AbstractStrategy.

References NOX::Abstract::Group::Ok, and NOX::ShapeCopy.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::TpetraHouseholder::initForTransposeSolve ( )

virtual

Intialize solver for a transpose solve.

This should be called after setMatrixBlocks(), but before applyInverseTranspose().

Implements LOCA::BorderedSolver::AbstractStrategy.

References NOX::Abstract::Group::Ok, and NOX::ShapeCopy.

void LOCA::BorderedSolver::TpetraHouseholder::setMatrixBlocks	(	const Teuchos::RCP< const LOCA::BorderedSolver::AbstractOperator > &	op,
		const Teuchos::RCP< const NOX::Abstract::MultiVector > &	blockA,
		const Teuchos::RCP< const LOCA::MultiContinuation::ConstraintInterface > &	blockB,
		const Teuchos::RCP< const NOX::Abstract::MultiVector::DenseMatrix > &	blockC
	)

virtual

Set blocks.

The blockA or blockC pointer may be null if either is zero. Whether block B is zero will be determined by querying blockB via ConstraintInterface::isConstraintDerivativesXZero.

Implements LOCA::BorderedSolver::AbstractStrategy.

References Teuchos::RCP< T >::get(), LOCA::BorderedSolver::JacobianOperator::getGroup(), Teuchos::SerialDenseMatrix< OrdinalType, ScalarType >::putScalar(), Teuchos::rcp(), and TEUCHOS_TEST_FOR_EXCEPTION.

void LOCA::BorderedSolver::TpetraHouseholder::updateCrsMatrixForPreconditioner	(	const NOX::Abstract::MultiVector &	U,
		const NOX::Abstract::MultiVector &	V,
		NOX::TCrsMatrix &	mat
	)		const

Overwrites the Jacobian $J$ with $J + U V^T$ for computing the preconditioner of $P$ .

NOTE: This should be a protected method, but cuda lambda forces this to be public!

References TEUCHOS_ASSERT.

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

LOCA_BorderedSolver_TpetraHouseholder.hpp
LOCA_BorderedSolver_TpetraHouseholder.cpp

Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation