Belos Version of the Day
Public Types | Public Member Functions
Belos::GmresBase< Scalar, MV, OP > Class Template Reference

Common state and functionality for new GMRES implementations. More...

#include <BelosGmresBase.hpp>

Inheritance diagram for Belos::GmresBase< Scalar, MV, OP >:
Inheritance graph
[legend]

List of all members.

Public Types

typedef Scalar scalar_type
typedef Teuchos::ScalarTraits
< Scalar >::magnitudeType 
magnitude_type
typedef MV multivector_type
typedef OP operator_type

Public Member Functions

 GmresBase (const Teuchos::RCP< LinearProblem< Scalar, MV, OP > > &problem, const Teuchos::RCP< const OrthoManager< Scalar, MV > > &ortho, const Teuchos::RCP< OutputManager< Scalar > > &outMan, const int maxIterCount, const bool flexible)
 Constructor.
bool canAdvance () const
 Whether it is legal to call advance()
void advance ()
 Perform the next group of iteration(s)
virtual void acceptCandidateBasis (const int newNumVectors=1)
 Accept the candidate basis and prepare for the next iteration.
virtual void rejectCandidateBasis (const int newNumVectors=1)
 Reject the candidate basis.
int getNumIters () const
 The number of completed iteration(s) (>= 0)
int maxNumIters () const
 Maximum number of iterations allowed before restart.
void restart (const int maxIterCount)
 Restart, and optionally reset max iteration count.
void restart ()
 Restart with the same max iteration count as before.
void backOut (const int numIters)
 "Back out" iterations following numIters
Teuchos::RCP< MV > getCurrentUpdate ()
 Compute and return current solution update.
Teuchos::RCP< MV > currentNativeResidualVector ()
 The current "native" residual vector.
Teuchos::ScalarTraits< Scalar >
::magnitudeType 
currentNativeResidualNorm ()
 The current "native" residual norm.
Teuchos::ScalarTraits< Scalar >
::magnitudeType 
arnoldiRelationError () const
 Forward normwise error in the Arnoldi relation.
void updateSolution ()
 Update the current approximate solution.

Protected Member Functions

Subclass-specific implementation details

Subclasses, by implementing all of these methods, implement a specific (F)GMRES algorithm.

virtual bool canExtendBasis () const =0
 Whether it's legal to call extendBasis()
virtual void extendBasis (Teuchos::RCP< MV > &V_cur, Teuchos::RCP< MV > &Z_cur)=0
 Extend the basis by 1 (or more) vector(s)
virtual void orthogonalize (const Teuchos::RCP< MV > &V_cur, const Teuchos::RCP< const MV > &V_prv, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &C_V, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &B_V, const Teuchos::RCP< MV > &Z_cur, const Teuchos::RCP< const MV > &Z_prv, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &C_Z, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &B_Z)=0
 Orthogonalize the candidate basis/es.
virtual void updateUpperHessenbergMatrix (const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &C_V, const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &B_V, const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &C_Z, const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &B_Z)=0
 Update the upper Hessenberg matrix (H_)
virtual bool acceptedCandidateBasis () const =0
 Whether the subclass "accepted" the candidate basis.
virtual bool shouldRetryAdvance ()
 Whether advance() should retry advancing the basis.
Hooks

Hooks allow subclasses to preface and postface restarting, and/or the orthogonalization of new basis vector(s), with additional features.

Typical uses would be implementing Recycling GMRES (GCRO-DR), or otherwise ensuring orthogonality of the computed Krylov basis/es to another given basis. GmresBase provides trivial default implementations of all the hooks.

Note:
In terms of aspect-oriented programming, restart() and orthogonalize() are "join points," and the hooks are "advice." Unfortunately we have to implement the "pointcuts" by hand, since otherwise we would need C++ syntax extensions.
If you know something about Emacs Lisp, you would say these hooks "advice" the restart() and orthogonalize() methods.
Warning:
Subclasses are responsible for invoking (or not invoking) their parent class' hook implementation, and for choosing when to do so. Each child class need only invoke its immediate parent's hook; recursion will take care of the rest. However, if your class inherits from multiple classes that all inherit from a subclass of GmresBase ("diamonds" in the class hierarchy), you'll need to resolve the hooks manually. (Otherwise, recursion will invoke the same hooks multiple times.)
virtual void preRestartHook (const int maxIterCount)
 Hook for subclass features before restarting.
virtual void postRestartHook ()
 Hook for subclass features after restarting.
virtual void preOrthogonalizeHook (const Teuchos::RCP< MV > &V_cur, const Teuchos::RCP< MV > &Z_cur)
 Hook to be invoked in advance() before orthogonalize().
void postOrthogonalizeHook (const Teuchos::RCP< MV > &V_cur, const Teuchos::RCP< mat_type > &C_V, const Teuchos::RCP< mat_type > &B_V, const Teuchos::RCP< MV > &Z_cur, const Teuchos::RCP< mat_type > &C_Z, const Teuchos::RCP< mat_type > &B_Z)
 Hook to be invoked in advance() after orthogonalize().

Protected Attributes

Member data, inherited by subclasses
Teuchos::RCP< LinearProblem
< Scalar, MV, OP > > 
lp_
 The linear problem to solve.
Teuchos::RCP< const
OrthoManager< Scalar, MV > > 
ortho_
 Orthogonalization manager.
Teuchos::RCP< OutputManager
< Scalar > > 
outMan_
 Output manager.
Teuchos::RCP< MV > nativeResVec_
 "Native" residual vector
Teuchos::RCP< MV > initResVec_
 Initial residual vector.
Teuchos::RCP< MV > xUpdate_
 Current solution update.
Teuchos::RCP< MV > V_
 First set of (orthogonalized) Krylov basis vectors.
Teuchos::RCP< MV > Z_
 Second set of (orthogonalized) Krylov basis vectors.
Teuchos::RCP
< details::ProjectedLeastSquaresProblem
< Scalar > > 
projectedProblem_
 The projected least-squares problem.
Teuchos::ScalarTraits< Scalar >
::magnitudeType 
initialResidualNorm_
 The initial residual norm.
Teuchos::ScalarTraits< Scalar >
::magnitudeType 
nativeResidualNorm_
 The "native" residual norm.
int lastUpdatedCol_
 Last column of H_ for which the QR factorization (implicitly stored in theCosines_, theSines_, and R_) has been computed.
int curNumIters_
 Current number of completed iterations.
int maxNumIters_
 Maximum number of iterations allowed.
bool flexible_
 Whether we are running Flexible GMRES.

Detailed Description

template<class Scalar, class MV, class OP>
class Belos::GmresBase< Scalar, MV, OP >

Common state and functionality for new GMRES implementations.

Author:
Mark Hoemmen
Warning:
This is EXPERIMENTAL CODE. DO NOT RELY ON THIS CODE. The interface or implementation may change at any time.

This class includes both state and functionality that are useful for different implementations of the Generalized Minimal Residual (GMRES) method of Saad and Schultz, for iterative solution of nonsymmetric nonsingular linear systems. Both the original GMRES algorithm (Saad and Schultz 1986) and Flexible GMRES (FGMRES) (Saad 1993) are supported. This class does not implement the actual iterations; this is left for subclasses. Furthermore, it does not implement features like recycling, but it does include hooks for subclasses to add in such features.

References:

GmresBase is an implementation class. The Belos solution manager for new GMRES implementations will interact with it mainly through a wrapper, which is a subclass of Belos::Iteration. Interaction with GMRES through an Iteration subclass will allow GMRES implementations to use Belos' stopping criteria (subclasses of StatusTest) and output methods (subclasses of StatusTestOutput).

Definition at line 226 of file BelosGmresBase.hpp.


Member Typedef Documentation

template<class Scalar , class MV , class OP >
typedef Scalar Belos::GmresBase< Scalar, MV, OP >::scalar_type
template<class Scalar , class MV , class OP >
typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Belos::GmresBase< Scalar, MV, OP >::magnitude_type
template<class Scalar , class MV , class OP >
typedef MV Belos::GmresBase< Scalar, MV, OP >::multivector_type
template<class Scalar , class MV , class OP >
typedef OP Belos::GmresBase< Scalar, MV, OP >::operator_type

Constructor & Destructor Documentation

template<class Scalar , class MV , class OP >
Belos::GmresBase< Scalar, MV, OP >::GmresBase ( const Teuchos::RCP< LinearProblem< Scalar, MV, OP > > &  problem,
const Teuchos::RCP< const OrthoManager< Scalar, MV > > &  ortho,
const Teuchos::RCP< OutputManager< Scalar > > &  outMan,
const int  maxIterCount,
const bool  flexible 
)

Constructor.

Construct and initialize a new GmresBase iteration to solve a given linear problem.

Parameters:
problem[in/out] Linear problem. On input, we use the starting guess (x0) and the initial residual (r0) to initialize the iteration. The advance() method does _not_ call updateSolution(); it only iterates on the solution update (the "delta"), not the initial guess. The restart() method _does_ call updateSolution().
ortho[in] Orthogonalization manager
outMan[in/out] Output manager
maxIterCount[in] Maximum number of iterations before restart. The number of vectors' worth of storage this constructor allocates is proportional to this, so choose carefully.
flexible[in] Whether or not to run the Flexible variant of GMRES (FGMRES). This requires twice as much vector storage, but lets the preconditioner change in every iteration. This only works with right preconditioning, not left or split preconditioning.
Note:
Instantiating a new GmresBase subclass instance is the only way to change the linear problem. This ensures that the V (and Z, for Flexible GMRES) spaces (especially their data distributions) are correct. It's allowed to change a LinearProblem instance's operator and right-hand, for example, and afterwards they might not even have the same number of rows, let alone the same data distribution as before. Belos::MultiVecTraits offers a limited interface to multivectors and operators that doesn't let us access the data distribution, so we choose instead to reinitialize from scratch whenever the linear problem is changed. The way to avoid this reinitialization would be to add a check for compatibility of multivectors in MultiVecTraits.

Definition at line 1113 of file BelosGmresBase.hpp.


Member Function Documentation

template<class Scalar , class MV , class OP >
bool Belos::GmresBase< Scalar, MV, OP >::canAdvance ( ) const [inline]

Whether it is legal to call advance()

Definition at line 286 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::advance ( )

Perform the next group of iteration(s)

Perform the next group of iteration(s), without any convergence or termination tests. Throw GmresCantExtendBasis if there is no more room for basis vectors; throw GmresRejectsCandidateBasis if the candidate basis vector(s) can't be accepted.

Note:
Different subclasses may compute different numbers of iteration(s) in this method. Standard (F)GMRES attempts exactly one iteration, whereas CA-(F)GMRES does work equivalent to several iterations of standard (F)GMRES.

Definition at line 951 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::acceptCandidateBasis ( const int  newNumVectors = 1) [inline, virtual]

Accept the candidate basis and prepare for the next iteration.

Note:
Subclasses may override this implementation if they need to do something unusual with the new candidate basis, like checkpoint it.

Reimplemented in Belos::CaGmres< Scalar, MV, OP >.

Definition at line 309 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::rejectCandidateBasis ( const int  newNumVectors = 1) [inline, virtual]

Reject the candidate basis.

Note:
Subclasses that retry rejected iterations should override this implementation appropriately.

Reimplemented in Belos::CaGmres< Scalar, MV, OP >.

Definition at line 325 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
int Belos::GmresBase< Scalar, MV, OP >::getNumIters ( ) const [inline]

The number of completed iteration(s) (>= 0)

Does not include the starting vector.

Definition at line 337 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
int Belos::GmresBase< Scalar, MV, OP >::maxNumIters ( ) const [inline]

Maximum number of iterations allowed before restart.

Does not include the starting vector.

Definition at line 342 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::restart ( const int  maxIterCount)

Restart, and optionally reset max iteration count.

First run preRestartHook(). Then, update the LinearProblem with the current approximate solution, and let it recompute the (exact, not "native") residual. Restart, reallocating space for basis vectors and the projected least-squares problem if necessary. Finally, run postRestartHook(). The hooks let subclasses implement things like subspace recycling.

Parameters:
maxIterCount[in] New maximum iteration count for the next restart cycle. If same as the old iteration count (the default), we reuse existing space for the Krylov basis vectors.
Note:
The restart will be invalid if the LinearProblem changes between restart cycles (except perhaps for the preconditioner if we are running Flexible GMRES). It's impossible for iteration classes to enforce this requirement, since LinearProblem objects are mutable. If you want to solve a different linear problem (e.g., with a different matrix A or right-hand side b), instantiate a new subclass of GmresBase. Changing the preconditioner is allowed if the flexible option is enabled.

Definition at line 1643 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::restart ( ) [inline]

Restart with the same max iteration count as before.

Definition at line 370 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::backOut ( const int  numIters)

"Back out" iterations following numIters

"Back out" iterations after, but not including, numIters. This is permanent, but relatively inexpensive. The main cost is $O(m^2)$ floating-point operations on small dense matrices and vectors, where m = numIters().

Parameters:
numIters[in] 0 <= numIters <= getNumIters().

Definition at line 1619 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP< MV > Belos::GmresBase< Scalar, MV, OP >::getCurrentUpdate ( )

Compute and return current solution update.

Compute and return current solution update. This also updates the value returned by currentNativeResidualNorm().

Definition at line 1392 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP< MV > Belos::GmresBase< Scalar, MV, OP >::currentNativeResidualVector ( )

The current "native" residual vector.

Compute the current "native" residual vector, which is formed from the basis vectors V, the upper Hessenberg matrix H, the current projected least-squares solution y, and the initial residual vector. Computing this does not invoke the operator or preconditioner(s).

Definition at line 1215 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::ScalarTraits< Scalar >::magnitudeType Belos::GmresBase< Scalar, MV, OP >::currentNativeResidualNorm ( )

The current "native" residual norm.

This is computed from the projected least-squares problem involving the upper Hessenberg matrix. Calling this does not invoke the operator or preconditioner(s), nor does it require computations on basis vectors.

Definition at line 1322 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::ScalarTraits< Scalar >::magnitudeType Belos::GmresBase< Scalar, MV, OP >::arnoldiRelationError ( ) const

Forward normwise error in the Arnoldi relation.

Compute and return the forward normwise error in the Arnoldi relation. For a fixed preconditioner (nonflexible GMRES), we compute

$\| \tilde{A} V_m - V_{m+1} \underline{H}_m \|_F$,

where $\tilde{A}$ is the preconditioned matrix. (We use the Frobenius norm to minimize dependency on the factorization codes which would be necessary for computing the 2-norm.) For a varying right preconditioner (Flexible GMRES), we compute

$\| A Z_m - V_{m+1} \underline{H}_m \|_F$.

In these expressions, "m" is the current iteration count, as returned by getNumIters().

Definition at line 1064 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::updateSolution ( ) [inline]

Update the current approximate solution.

Modify the LinearProblem instance's current approximate solution (the multivector returned by getLHS()) by adding in the current solution update (xUpdate_). Tell the LinearProblem to compute a new residual vector as well.

Definition at line 433 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual bool Belos::GmresBase< Scalar, MV, OP >::canExtendBasis ( ) const [protected, pure virtual]

Whether it's legal to call extendBasis()

Whether the basis (bases, if Flexible GMRES) can be extended by computing candidate basis vectors, i.e., whether it's legal to call extendBasis(). Subclasses' implementations of extendBasis() must throw GmresCantExtendBasis if no more candidate basis vectors can be computed.

Implemented in Belos::CaGmres< Scalar, MV, OP >, and Belos::StandardGmres< Scalar, MV, OP >.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::extendBasis ( Teuchos::RCP< MV > &  V_cur,
Teuchos::RCP< MV > &  Z_cur 
) [protected, pure virtual]

Extend the basis by 1 (or more) vector(s)

Extend the basis (bases, if Flexible GMRES) by 1 (or more) vector(s), without orthogonalizing it/them. Standard GMRES only adds one basis vector at a time; CA-GMRES may add more than one at a time.

Parameters:
V_cur[out] New basis vector(s). The number of column(s) gives the number of basis vector(s) added. This may be a view of member data, rather than freshly allocated storage.
Z_cur[out] If running Flexible GMRES, the preconditioned basis vector(s); else, Teuchos::null. The number of column(s) gives the number of basis vector(s) added. This may be a view of member data, rather than freshly allocated storage.
Warning:
Don't call this if canExtendBasis() returns false.
Note:
This method is non-const so that implementations may use previously allocated space for the basis vectors. (They are not _required_ to use previously allocated space, so this method may, if it wishes, allocate a new multivector to return.)

Implemented in Belos::CaGmres< Scalar, MV, OP >, and Belos::StandardGmres< Scalar, MV, OP >.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::orthogonalize ( const Teuchos::RCP< MV > &  V_cur,
const Teuchos::RCP< const MV > &  V_prv,
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  C_V,
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  B_V,
const Teuchos::RCP< MV > &  Z_cur,
const Teuchos::RCP< const MV > &  Z_prv,
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  C_Z,
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  B_Z 
) [protected, pure virtual]

Orthogonalize the candidate basis/es.

Flexible CA-GMRES must orthogonalize the Z basis along with the V basis, and it must keep both sets of orthogonalization coefficients. Flexible standard GMRES, in contrast, only needs to orthogonalize the new V basis vector; it doesn't need to orthogonalize the Z basis vector. This is why we let the subclass implement this functionality, rather than calling the OrthoManager's projectAndNormalize() method directly; the subclass has to decide whether it needs to orthogonalize Z_cur as well as V_cur.

Subclasses may want to save the rank (if applicable) from ortho_->projectAndNormalize() somewhere. In the case of Flexible CA-GMRES, there are _two_ ranks -- one for the V basis, and the other for the Z basis. This is why this method doesn't return the rank, as ortho_->projectAndNormalize() does.

C_V and B_V are the "C" (projection) resp. "B" (normalization) block coefficients from ortho_->projectAndNormalize() on the new V basis vectors. If the new Z basis vectors need to be orthogonalized as well, then C_Z and B_Z are the "C" resp. "B" block coefficients from ortho_->projectAndNormalize() on the new Z basis vectors. Subclasses' implementations of orthogonalize() are responsible for allocating space for C_V and B_V, and C_Z and B_Z if necessary (otherwise they are set to Teuchos::null, as the RCPs are passed by reference).

Parameters:
V_cur[in/out] Candidate basis vector(s) for V basis
V_prv[in] Previously orthogonalized V basis vectors
C_V[out] Projection coefficients (V_prv^* V_cur) of candidate V basis vector(s)
B_V[out] Normalization coefficients (after projection) of candidate V basis vector(s)
Z_cur[in/out] Candidate basis vector(s) for Z basis, or null if there is no Z basis
Z_prv[in] Previously orthogonalized Z basis vectors, or null if there is no Z basis
C_Z[out] Projection coefficients (Z_prv^* Z_cur) of candidate Z basis vector(s), or null if there is no Z basis
B_Z[out] Normalization coefficients (after projection) of candidate Z basis vector(s), or null if there is no Z basis

Implemented in Belos::StandardGmres< Scalar, MV, OP >.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::updateUpperHessenbergMatrix ( const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  C_V,
const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  B_V,
const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  C_Z,
const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > &  B_Z 
) [protected, pure virtual]

Update the upper Hessenberg matrix (H_)

Parameters:
C_V[in] Projection coefficients (V_prv^* V_cur) of candidate V basis vector(s)
B_V[in] Normalization coefficients (after projection) of candidate V basis vector(s)
C_Z[in] Projection coefficients (Z_prv^* Z_cur) of candidate Z basis vector(s), or null if there is no Z basis
B_Z[in] Normalization coefficients (after projection) of candidate Z basis vector(s), or null if there is no Z basis
Note:
For Flexible GMRES, it may be desirable to compute or update a rank-revealing decomposition of the upper square submatrix of H_. This is because FGMRES only promises convergence in exact arithmetic if this submatrix is full rank and the computed residual norm is zero. Any decomposition of H_ should not be done in place; this is because updateProjectedLeastSquaresProblem() depends on H_ being intact. That method does not itself modify H_, so implementations of updateUpperHessenbergMatrix() can also rely on H_ being intact.
For an algorithm for updating a rank-revealing decomposition, see e.g., G. W. Stewart, "Updating a rank-revealing ULV decomposition", SIAM J. Matrix Anal. & Appl., Volume 14, Issue 2, pp. 494-499 (April 1993).

Implemented in Belos::CaGmres< Scalar, MV, OP >, and Belos::StandardGmres< Scalar, MV, OP >.

template<class Scalar , class MV , class OP >
virtual bool Belos::GmresBase< Scalar, MV, OP >::acceptedCandidateBasis ( ) const [protected, pure virtual]

Whether the subclass "accepted" the candidate basis.

Implemented in Belos::CaGmres< Scalar, MV, OP >, and Belos::StandardGmres< Scalar, MV, OP >.

template<class Scalar , class MV , class OP >
virtual bool Belos::GmresBase< Scalar, MV, OP >::shouldRetryAdvance ( ) [inline, protected, virtual]

Whether advance() should retry advancing the basis.

In advance(), if the subclass rejects the candidate basis vectors (acceptedCandidateBasis() returns false), then we have to decide whether to retry the current iteration. Retry is done by calling advance() recursively. Subclasses should set state before advance() calls acceptedCandidateBasis(), to ensure that this recursion terminates either with success, or by throwing GmresRejectsCandidateBasis.

Default behavior is never to retry. This always results in finite termination of advance().

Definition at line 587 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::preRestartHook ( const int  maxIterCount) [inline, protected, virtual]

Hook for subclass features before restarting.

This method is called by restart(), before restart() does anything. The default implementation does nothing. This would be the place to add things like Krylov subspace recycling (the part where the recycling subspace is computed from the existing basis).

Note:
Implementations of this method should not update the linear problem; this happens immediately following the invocation of preRestartHook().

Definition at line 633 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::postRestartHook ( ) [inline, protected, virtual]

Hook for subclass features after restarting.

This method is called by restart(), after restart() finishes everything else. The default implementation does nothing. This would be the place where subclasses may implement things like orthogonalizing the first basis vector (after restart) against the recycled subspace.

Note:
The initial residual has already been set by restart() before this method is called.

Definition at line 645 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
virtual void Belos::GmresBase< Scalar, MV, OP >::preOrthogonalizeHook ( const Teuchos::RCP< MV > &  V_cur,
const Teuchos::RCP< MV > &  Z_cur 
) [inline, protected, virtual]

Hook to be invoked in advance() before orthogonalize().

Hook called in advance() before calling orthogonalize() on the given newly generated basis vectors V_cur (of the V basis) and Z_cur (of the Z basis). The default implementation does nothing.

Note:
Subclasses that perform subspace recycling could use implement orthogonalization against the recycled subspace either here or in postOrthogonalizeHook().

Definition at line 658 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
void Belos::GmresBase< Scalar, MV, OP >::postOrthogonalizeHook ( const Teuchos::RCP< MV > &  V_cur,
const Teuchos::RCP< mat_type > &  C_V,
const Teuchos::RCP< mat_type > &  B_V,
const Teuchos::RCP< MV > &  Z_cur,
const Teuchos::RCP< mat_type > &  C_Z,
const Teuchos::RCP< mat_type > &  B_Z 
) [inline, protected]

Hook to be invoked in advance() after orthogonalize().

Hook called in advance() after calling orthogonalize() on the given newly generated basis vectors V_cur (of the V basis) and Z_cur (of the Z basis). The default implementation does nothing.

Note:
Subclasses that perform subspace recycling could use implement orthogonalization against the recycled subspace either here or in preOrthogonalizeHook().

Definition at line 672 of file BelosGmresBase.hpp.


Member Data Documentation

template<class Scalar , class MV , class OP >
Teuchos::RCP<LinearProblem<Scalar, MV, OP> > Belos::GmresBase< Scalar, MV, OP >::lp_ [protected]

The linear problem to solve.

Definition at line 765 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<const OrthoManager<Scalar, MV> > Belos::GmresBase< Scalar, MV, OP >::ortho_ [protected]

Orthogonalization manager.

Definition at line 768 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<OutputManager<Scalar> > Belos::GmresBase< Scalar, MV, OP >::outMan_ [protected]

Output manager.

Definition at line 771 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<MV> Belos::GmresBase< Scalar, MV, OP >::nativeResVec_ [protected]

"Native" residual vector

"Native" means computed using exact-arithmetic invariants of the iterative method. In this case, if m is the current number of iterations and $r_0 = A x_0 - b$,

\begin{eqnarray*} A x_k - b &=& A (x_0 + Z(1:m) y(1:m)) - b \\ &=& r_0 + A Z(1:m) y(1:m) \\ &=& r_0 + V(1:m+1) H(1:m+1, 1:m) y(1:m). \end{eqnarray*}

Accuracy of the above formula depends only on the Arnoldi relation, which in turn depends mainly on the residual error of the orthogonalization being small. This is true for just about any orthogonalization method, so the above computation should hold just about always.

Storage for the native residual vector is allocated only on demand. The storage is cached and reallocated only when the LinearProblem changes. (This ensures that the dimensions and data distribution are correct, i.e., are the same as the initial residual vector in the linear problem).

Definition at line 796 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<MV> Belos::GmresBase< Scalar, MV, OP >::initResVec_ [protected]

Initial residual vector.

This is the left-preconditioned residual vector if a left preconditioner has been provided, otherwise it is the unpreconditioned residual vector.

Definition at line 803 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<MV> Belos::GmresBase< Scalar, MV, OP >::xUpdate_ [protected]

Current solution update.

The current approximate solution for Flexible GMRES is x0 + Z_[0:numIters-1]y, and for (nonflexible) GMRES is x0 + V_[0:numIters-1]y.

Note:
Invalidated when linear problem changes, etc.

Definition at line 812 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<MV> Belos::GmresBase< Scalar, MV, OP >::V_ [protected]

First set of (orthogonalized) Krylov basis vectors.

Note:
These vectors are always used, whether we are performing Flexible GMRES (left preconditing with a possibly changing preconditioner, so we keep a second set of basis vectors Z_) or standard GMRES.

Definition at line 820 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<MV> Belos::GmresBase< Scalar, MV, OP >::Z_ [protected]

Second set of (orthogonalized) Krylov basis vectors.

Note:
These vectors are only used if we are performing Flexible GMRES. In that case, the solution update coefficients are computed for this basis.

Definition at line 827 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::RCP<details::ProjectedLeastSquaresProblem<Scalar> > Belos::GmresBase< Scalar, MV, OP >::projectedProblem_ [protected]

The projected least-squares problem.

GMRES computes the solution update's coefficients by solving a small dense least-squares problem. The latter is the projection of the full residual norm minimization problem onto the Krylov search space.

Definition at line 835 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::ScalarTraits<Scalar>::magnitudeType Belos::GmresBase< Scalar, MV, OP >::initialResidualNorm_ [protected]

The initial residual norm.

GMRES makes use of the initial residual norm for solving the projected least-squares problem for the solution update coefficients. In that case, it is usually called $\beta$. For left-preconditioned GMRES, this is the preconditioned initial residual norm, else it's the unpreconditioned version.

Definition at line 844 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
Teuchos::ScalarTraits<Scalar>::magnitudeType Belos::GmresBase< Scalar, MV, OP >::nativeResidualNorm_ [protected]

The "native" residual norm.

Definition at line 852 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
int Belos::GmresBase< Scalar, MV, OP >::lastUpdatedCol_ [protected]

Last column of H_ for which the QR factorization (implicitly stored in theCosines_, theSines_, and R_) has been computed.

Should be initialized to -1 (which means the 0th column of H_ will be the first to be updated).

Definition at line 858 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
int Belos::GmresBase< Scalar, MV, OP >::curNumIters_ [protected]

Current number of completed iterations.

This also counts the number of orthogonalized basis vectors in V_. Columns [0, curNumIters_] of V_ are always orthogonalized basis vectors.

Definition at line 865 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
int Belos::GmresBase< Scalar, MV, OP >::maxNumIters_ [protected]

Maximum number of iterations allowed.

This affects the amount of storage consumed by V_ and Z_. Our Arnoldi/GMRES implementation requires that all the basis vectors in a set occupy a contiguous block of storage. (This differs from some other Belos implementations of iterative methods, where each basis vector (or block of vectors) is allocated separately on demand.)

Definition at line 875 of file BelosGmresBase.hpp.

template<class Scalar , class MV , class OP >
bool Belos::GmresBase< Scalar, MV, OP >::flexible_ [protected]

Whether we are running Flexible GMRES.

Flexible GMRES (FGMRES) is a variant of standard GMRES that allows the preconditioner to change in every iteration. It only works for right preconditioning. FGMRES requires keeping two sets of Krylov basis vectors: one for the Krylov subspace

\[ \text{span}\{ r_0, A M^{-1} r_0, \dots, (A M^{-1})^k r_0 \}, \]

where $r_0$ is the unpreconditioned residual, and one for the Krylov subspace

\[ \text{span}\{ M^{-1} r_0, M^{-1} A M^{-1} r_0, \dots, M^{-1} (A M^{-1})^{k-1} r_0 \}. \]

We store the basis for the first subspace in V_, and the basis for the second subspace in Z_. If we are not running FGMRES, we let Z_ be Teuchos::null. FGMRES reduces to standard GMRES in the case of nopreconditioning, and then we also let Z_ be null.

Note:
In the original Flexible GMRES paper, Saad suggested implementing FGMRES and GMRES in a single routine, with a runtime switch to decide whether to keep the Z basis and whether update the solution with the Q or Z basis. This is in fact what we do. The run-time cost per iteration is no more than that of checking this Boolean a small constant number of times.

Definition at line 906 of file BelosGmresBase.hpp.


The documentation for this class was generated from the following file:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines