LOCA::BorderedSolver::EpetraHouseholder Class Reference

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

#include <LOCA_BorderedSolver_EpetraHouseholder.H>

Inheritance diagram for LOCA::BorderedSolver::EpetraHouseholder:

[legend]
Collaboration diagram for LOCA::BorderedSolver::EpetraHouseholder:
[legend]
List of all members.

Public Member Functions

 EpetraHouseholder (const Teuchos::RefCountPtr< LOCA::GlobalData > &global_data, const Teuchos::RefCountPtr< LOCA::Parameter::SublistParser > &topParams, const Teuchos::RefCountPtr< Teuchos::ParameterList > &solverParams)
 Constructor.
virtual ~EpetraHouseholder ()
 Destructor.
virtual void setMatrixBlocks (const Teuchos::RefCountPtr< const NOX::Abstract::Group > &group, const Teuchos::RefCountPtr< const NOX::Abstract::MultiVector > &blockA, const Teuchos::RefCountPtr< const LOCA::MultiContinuation::ConstraintInterface > &blockB, const Teuchos::RefCountPtr< const NOX::Abstract::MultiVector::DenseMatrix > &blockC)
 Set blocks.
virtual NOX::Abstract::Group::ReturnType initForSolve ()
 Intialize solver for a solve.
virtual NOX::Abstract::Group::ReturnType initForTransposeSolve ()
 Intialize solver for a transpose solve.
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.
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.
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.
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.

Protected Member Functions

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)
void updateJacobianForPreconditioner (const NOX::Abstract::MultiVector &U, const NOX::Abstract::MultiVector &V, Epetra_CrsMatrix &jac) const
 Overwrites the Jacobian $J$ with $J + U V^T$ for computing the preconditioner of $P$.

Protected Attributes

Teuchos::RefCountPtr< LOCA::GlobalDataglobalData
 Global data object.
Teuchos::RefCountPtr< Teuchos::ParameterListsolverParams
 Solver parameters.
Teuchos::RefCountPtr< const
LOCA::Epetra::Group
grp
 Pointer to group storing J.
Teuchos::RefCountPtr< const
NOX::Abstract::MultiVector
A
 Pointer to A block.
Teuchos::RefCountPtr< const
NOX::Abstract::MultiVector
B
 Pointer to B block.
Teuchos::RefCountPtr< const
NOX::Abstract::MultiVector::DenseMatrix
C
 Pointer to C block.
Teuchos::RefCountPtr< const
LOCA::MultiContinuation::ConstraintInterfaceMVDX
constraints
 Pointer to constraint interface.
LOCA::BorderedSolver::HouseholderQR qrFact
 QR Factorization object.
Teuchos::RefCountPtr< NOX::Abstract::MultiVectorhouse_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::RefCountPtr< NOX::Abstract::MultiVectorU
 U matrix in low-rank update form P = J + U*V^T.
Teuchos::RefCountPtr< NOX::Abstract::MultiVectorV
 V matrix in low-rank update form P = J + U*V^T.
Teuchos::RefCountPtr< NOX::Abstract::MultiVectorhouse_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::RefCountPtr< NOX::Abstract::MultiVectorU_trans
 U matrix in low-rank update form P = J + U*V^T for transposed system.
Teuchos::RefCountPtr< NOX::Abstract::MultiVectorV_trans
 V matrix in low-rank update form P = J + U*V^T for transposed system.
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 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.

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 T Y$ (see R. Schreiver and C. Van Loan, "A Storage-Efficient WY Represntation 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::Epetra::LowRankUpdateRowMatrix if $J$ is an Epetra_RowMatrix and LOCA::Epetra::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:


Constructor & Destructor Documentation

LOCA::BorderedSolver::EpetraHouseholder::EpetraHouseholder const Teuchos::RefCountPtr< LOCA::GlobalData > &  global_data,
const Teuchos::RefCountPtr< LOCA::Parameter::SublistParser > &  topParams,
const Teuchos::RefCountPtr< 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


Member Function Documentation

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::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.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::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.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::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.

Reimplemented from LOCA::BorderedSolver::AbstractStrategy.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::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.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::initForSolve  )  [virtual]
 

Intialize solver for a solve.

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

Implements LOCA::BorderedSolver::AbstractStrategy.

NOX::Abstract::Group::ReturnType LOCA::BorderedSolver::EpetraHouseholder::initForTransposeSolve  )  [virtual]
 

Intialize solver for a transpose solve.

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

Implements LOCA::BorderedSolver::AbstractStrategy.

void LOCA::BorderedSolver::EpetraHouseholder::setMatrixBlocks const Teuchos::RefCountPtr< const NOX::Abstract::Group > &  group,
const Teuchos::RefCountPtr< const NOX::Abstract::MultiVector > &  blockA,
const Teuchos::RefCountPtr< const LOCA::MultiContinuation::ConstraintInterface > &  blockB,
const Teuchos::RefCountPtr< 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.


The documentation for this class was generated from the following files:
Generated on Thu Sep 18 12:38:30 2008 for NOX by doxygen 1.3.9.1