#include <Thyra_LinearOpWithSolveBaseDecl.hpp>
Inheritance diagram for Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >:
Pure virtual functions that must be overridden in subclasses | |
virtual void | solve (const EConj conj, const MultiVectorBase< RangeScalar > &B, MultiVectorBase< DomainScalar > *X, const int numBlocks=0, const BlockSolveCriteria< PromotedScalar > blockSolveCriteria[]=NULL, SolveStatus< PromotedScalar > blockSolveStatus[]=NULL) const =0 |
Request the forward solution of a block system with different targeted solution criteria. | |
Virtual functions with default implementations | |
virtual bool | solveSupportsConj (EConj conj) const |
Return if solve() supports the argument conj . | |
virtual bool | solveTransposeSupportsConj (EConj conj) const |
Return if solveTranspose() supports the argument conj . | |
virtual bool | solveSupportsSolveMeasureType (EConj conj, const SolveMeasureType &solveMeasureType) const |
Return if solve() supports the given the solve measure type. | |
virtual bool | solveTransposeSupportsSolveMeasureType (EConj conj, const SolveMeasureType &solveMeasureType) const |
Return if solveTranspose() supports the given the solve measure type. | |
virtual void | solveTranspose (const EConj conj, const MultiVectorBase< DomainScalar > &B, MultiVectorBase< RangeScalar > *X, const int numBlocks=0, const BlockSolveCriteria< PromotedScalar > blockSolveCriteria[]=NULL, SolveStatus< PromotedScalar > blockSolveStatus[]=NULL) const |
Request the transpose (or adjoint) solution of a block system with different targeted solution criteria. | |
Public Types | |
typedef Teuchos::PromotionTraits< RangeScalar, DomainScalar >::promote | PromotedScalar |
Local typedef for promoted scalar type. |
solve()
) of the form:
and/or a transpose solve operation (using solveTranspose()
) of the form:
and/or an adjoint solve operation (using solveTranspose()
) of the form:
where is *this
linear operator, is an appropriate RHS multi-vector with columns, and is a LHS multi-vector with columns that is computed by this interface. Note that if the underlying operator has real-valued entries then the transpose and the adjoint are the same.
Note that this interface does not assume that the linear operator itself is nonsingular or invertible in the classic sense (i.e. an inverse operator may not exist).
Let signify either the forward operator , the transpose operator or the adjoint operator . What this interface assumes is that for any appropriately selected consistent multi-vector RHS B
that a solve of will yield an approximate solution LHS multi-vector X
such that A X == B
. Note that this interface does not assume that a solution can be computed for any random RHS multi-vector . Solutions for any random RHS can on be expected for relatively well conditioned non-singular operators.
Note: It is recommended that clients use the non-member helper functions defined here rather than call these member functions directly as they support a number of other simpler use cases.
*this
operator to exclude support certain types of solve measures (see the functions solveSupportsSolveMeasureType()
and solveTransposeSupportsSolveMeasureType()
). Also, this interface assumes that all implementations can support a "default" solve criteria that is determined internally to *this
.This interface is meant to support direct and iterative linear solvers as well as combinations of the two in a variety of configurations. Because of the almost infinite number of types of linear solver configurations possible, this interface tries not to specify any particular solver-specific solution control options. The one exception is a maximum number of iterations which is totally implementation defined. These types of control options are better specified in lower lever implementations and should be kept out of an interface such as this.
The functions solve()
and solveTranspose()
both take the arguments:
,const int numBlocks = 0 ,const BlockSolveCriteria<PromotedScalar> blockSolveCriteria[] = NULL ,SolveStatus<PromotedScalar> blockSolveStatus[] = NULL
The array arguments blockSolveCriteria[]
and blockSolveStatus[]
specify different blocks of solution criteria and the corresponding solve return statuses for a partitioned set of linear systems. Assuming that the client passes in arrays of dimensions numBlocks
, these tolerances define the solution criteria for the block systems:
where the column indexes are given by , for .
The solve criteria for the block system
is given by blockSolveCriteria[j]
(if blockSolveCriteria!=NULL
) and the solution status after return is given by blockSolveStatus[j]
(if blockSolveStatus!=NULL
).
By specifying solution criteria in blocks and then only requesting basic tolerances, we allow linear solver implementations every opportunity to perform as many optimizations as possible in solving the linear systems. For example, SVD could be performed on each block of RHSs and then a reduced set of linear systems could be solved.
For the remainder of this discussion we will focus on how the solution criteria for single block of linear systems specified by a single BlockSolveCriteria
object is specified and how the status of this block linear solve is reported in a SolveStatus
object.
The struct BlockSolveCriteria
contains a SolveCriteria
member and a number of RHSs that it applies to. It is the SolveCriteria
object that determines the type and tolerance of a block solve request.
Let solveCriteria
be a SolveCriteria
object setup by a client to be passed into a solve operation. This object can be set up in a variety of ways to support several different use cases which are described below:
solveCriteria.solveMeasureType.useDefault()==true
]: In this mode, criteria for the solution tolerance is determined internally by *this
object. Usually, it would be assumed that *this
would solve the linear systems to a sufficient tolerance for the given ANA client. This is the mode that many ANAs not designed to control inexactness would work with. In this mode, the solution criteria will be determined in the end by the application or the user in an "appropriate" manner. In this case, the value of solveCriteria.requestedTol
is ignored and no meaningful solve status can be returned to the client.
solveCriteria.solveMeasureType.numerator==SOLVE_MEASURE_NORM_RESIDUAL && solveCriteria.solveMeasureType.denominator==SOLVE_MEASURE_NORM_RHS && solveCriteria.requestedTol!=SolveCriteriaunspecifiedTolerance()
]: In this mode, the solution algorithm would be requested to solve the block system to the relative residual tolerance of
for , where solveCriteria.requestedTol
and where the norm is given by the natural norm defined by the range space of and computed from Thyra::norm()
. Many linear solvers should be able to monitor this tolerance and be able to achieve it, or on failure be able to report the actual tolerance achieved.
solve()
and solveTranspose()
functions return, the client can optionally get back a solution status for each block of linear systems for of block solve criteria. Specifically, for each block of linear systems
whose solution criteria is specified by a SolveCriteria
object, a SolveStatus
object can optionally be returned that lets the client know the status of the linear solve.
A note about direct solvers is in order. The "inexact" solve features of this interface are primarily designed to support "loose" solve tolerances that exploit the properties of iterative linear solvers. With that said, any decent direct solver can assume that it has met the convergence criteria as requested by the client but does not have to return an estimate of the actual tolerance achieved.
If solveStatus
is a SolveStatus
object returned for the above block linear system the the following return status are significant:
solveStatus.solveStatus==SOLVE_STATUS_CONVERGED
]: This status is returned by the linear solver if the solution criteria was likely achieved (within some acceptable cushion cased by round-off etc.). This should almost always be the return value for a direct linear solver. The maximum actual tolerance achieved may or may not be returned in the field solveStatus.achievedTol
. The two sub-cases are:
solveStatus.achievedTol!=SolveStatusunknownTolerance()
] : The linear solver knows the approximate order-of-magnitude estimate of the maximum tolerance achieved. An order-of-magnitude (or so) estimate of the achieved tolerance would likely be known by any iterative linear solver when solveCriteria.solveMeasureType.numerator==SOLVE_MEASURE_NORM_RESIDUAL&&solveCriteria.solveMeasureType.denominator==SOLVE_MEASURE_NORM_RHS
. Most direct linear solvers will not return a known tolerance except for the case where solveCriteria.solveMeasureType.numerator==SOLVE_MEASURE_NORM_RESIDUAL&&solveCriteria.solveMeasureType.denominator==SOLVE_MEASURE_NORM_RHS
and where iterative refinement is used.
solveStatus.achievedTol==SolveStatusunknownTolerance()
] : The linear solver does not know the tolerance that was achieved but the achieved tolerance should be very close to the requested tolerance. This would be the most likely return status for a direct linear solver.
solveStatus.solveStatus==SOLVE_STATUS_UNCONVERGED
]: The linear solver was most likely not able to achieve the requested tolerance. A direct linear solver should almost never return this status except for in extreme cases (e.g. highly ill conditioned matrix and a tight requested tolerance). The linear solver may not be able to return the actual tolerance achieved and the same two cases as for the <it>unconverged</it> case are possible and the two subdcases are:
solveStatus.achievedTol!===SolveStatusunknownTolerance()0
] : The linear solver knows the approximate order-of-magnitude estimate of the maximum tolerance achieved.
solveStatus.achievedTol==SolveStatusunknownTolerance()
] : The linear solver does not know the tolerance that was achieved but the achieved tolerance is most likely significantly greater than the requested tolerance.
solveStatus.solveStatus==SOLVE_STATUS_UNKNOWN
]: The linear solver does not know if the solution status was achieved or not. In this case, the value of solveStatus.achievedTol==SolveStatusunknownTolerance()
may or many not be returned. This may be the return value when there is no reasonable way that the linear solver algorithm can know now to compute or estimate the requested tolerance. This will also always be the return status when solveCriteria.solveMeasureType.useDefault()==true
since the client would have no way to interpret this tolerance. The value of solveStatus.achievedTol!=SolveStatusunknownTolerance()
in this case should only be returned when solveCriteria.solveMeasureType.useDefault()==true
and therefore the client would have no way to interpret this tolerance as a residual or an solution norm.
The implementation of the function accumulateSolveStatus()
defines how to accumulate the individual solve status for each RHS in a block into the overall solve status for an entire block returned by blockSolveStatus[]
.
ToDo: Finish documentation!
solve()
must be overridden. See LinearOpBase
for what other virtual functions must be overridden to define a concrete subclass.
Definition at line 338 of file Thyra_LinearOpWithSolveBaseDecl.hpp.
typedef Teuchos::PromotionTraits<RangeScalar,DomainScalar>::promote Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::PromotedScalar |
Local typedef for promoted scalar type.
Definition at line 345 of file Thyra_LinearOpWithSolveBaseDecl.hpp.
virtual void Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solve | ( | const EConj | conj, | |
const MultiVectorBase< RangeScalar > & | B, | |||
MultiVectorBase< DomainScalar > * | X, | |||
const int | numBlocks = 0 , |
|||
const BlockSolveCriteria< PromotedScalar > | blockSolveCriteria[] = NULL , |
|||
SolveStatus< PromotedScalar > | blockSolveStatus[] = NULL | |||
) | const [pure virtual] |
Request the forward solution of a block system with different targeted solution criteria.
conj | [in] Determines if the elements are non-conjugate (NONCONJ_ELE ) or conjugate (CONJ_ELE ). For real valued operator, this argument is meaningless. Most ANAs will request NONCONJ_ELE . | |
B | [in] The RHS multi-vector with m = B.domain()->dim() columns. | |
X | [in/out] The LHS multi-vector with with m = X->domain()->dim() columns. On input, contains the initial guess for the solution (only significant for iterative solvers) and on output contains an estimate of the solution. | |
numBlocks | [in] The number of blocks for which solve tolerances will be specified for. If numBlocks==0 then this is a flag that a default set of tolerances should be used for all the linear systems. Default numBlocks=0 . | |
blockSolveCriteria | [in] Array (length numBlocks ) which gives the desired solution criteria for each of the numBlocks blocks of RHS. If numBlocks>0 then this argument must be non-NULL and point to a valid array of numBlocks entries. | |
blockSolveStatus | [out] Array (length numBlocks ) which gives the status of each set of block systems. A value of blockSolveStatus==NULL is allowed and means that the client is not interested in solution status of the linear systems. |
this->solveSupportsConj(conj)==true
X!=NULL
this->range()->isCompatible(*B.range())==true
this->domain()->isCompatible(*X->range())==true
B->domain()->isCompatible(*X->domain())==true
numBlocks >= 0
numBlocks == 0
] blockSolveCriteria==NULL && blockSolveStatus==NULL
. numBlocks > 0
] blockSolveCriteria!=NULL
and points to an array of length at least numBlocks
. numBlocks > 0
] this->solveSupportsSolveMeasureType(conj,blockSolveCriteria[k].solveMeasureType)==true
, for k = 0...numBlocks-1
. blockSolveStatus!=NULL
] blockSolveStatus
and points to an array of length at least numBlocks
. Postconditions:
CatastrophicSolveFailure
will be thrown. If the function returns normally, then the following postconditions apply. blockSolveStatus!=NULL
then blockSolveStatus[k]
gives the solution status of the block of linear systems specified by blockSolveCriteria[k]
where, for k=0...numBlocks-1
.
See the above introduction for a more complete description of how this function behaves and the meaning of the arguments blockSolveCriteria[]
and blockSolveStatus[]
.
bool Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solveSupportsConj | ( | EConj | conj | ) | const [virtual] |
Return if solve()
supports the argument conj
.
The default implementation returns true
for real valued scalar types or when conj==NONCONJ_ELE
for complex valued types.
Reimplemented in Thyra::SingleScalarLinearOpWithSolveBase< Scalar >, RealComplexFFTLinearOp< RealScalar >, and Thyra::SingleScalarLinearOpWithSolveBase< std::complex< RealScalar > >.
Definition at line 38 of file Thyra_LinearOpWithSolveBase.hpp.
bool Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solveTransposeSupportsConj | ( | EConj | conj | ) | const [virtual] |
Return if solveTranspose()
supports the argument conj
.
The default implementation returns false
.
Reimplemented in Thyra::SingleScalarLinearOpWithSolveBase< Scalar >, RealComplexFFTLinearOp< RealScalar >, and Thyra::SingleScalarLinearOpWithSolveBase< std::complex< RealScalar > >.
Definition at line 44 of file Thyra_LinearOpWithSolveBase.hpp.
bool Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solveSupportsSolveMeasureType | ( | EConj | conj, | |
const SolveMeasureType & | solveMeasureType | |||
) | const [virtual] |
Return if solve()
supports the given the solve measure type.
The default implementation returns true
for solveMeasureType.inNone()
.
Reimplemented in Thyra::SingleScalarLinearOpWithSolveBase< Scalar >, and Thyra::SingleScalarLinearOpWithSolveBase< std::complex< RealScalar > >.
Definition at line 50 of file Thyra_LinearOpWithSolveBase.hpp.
bool Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solveTransposeSupportsSolveMeasureType | ( | EConj | conj, | |
const SolveMeasureType & | solveMeasureType | |||
) | const [virtual] |
Return if solveTranspose()
supports the given the solve measure type.
The default implementation returns true
for solveMeasureType.inNone()
.
Reimplemented in Thyra::SingleScalarLinearOpWithSolveBase< Scalar >, and Thyra::SingleScalarLinearOpWithSolveBase< std::complex< RealScalar > >.
Definition at line 56 of file Thyra_LinearOpWithSolveBase.hpp.
void Thyra::LinearOpWithSolveBase< RangeScalar, DomainScalar >::solveTranspose | ( | const EConj | conj, | |
const MultiVectorBase< DomainScalar > & | B, | |||
MultiVectorBase< RangeScalar > * | X, | |||
const int | numBlocks = 0 , |
|||
const BlockSolveCriteria< PromotedScalar > | blockSolveCriteria[] = NULL , |
|||
SolveStatus< PromotedScalar > | blockSolveStatus[] = NULL | |||
) | const [virtual] |
Request the transpose (or adjoint) solution of a block system with different targeted solution criteria.
conj | [in] Determines if the elements are non-conjugate (NONCONJ_ELE ) or conjugate (CONJ_ELE ). For real valued operator, this argument is meaningless. The transpose solve is requested with conj==NONCONJ_ELE and the adjoint solve is requested with conj==CONJ_ELE . | |
B | [in] The RHS multi-vector with m = B.domain()->dim() columns. | |
X | [in/out] The LHS multi-vector with with m = X->domain()->dim() columns. On input, contains the initial guess for the solution (only significant for iterative solvers) and on output contains an estimate of the solution. | |
numBlocks | [in] The number of blocks for which solve tolerances will be specified for. If numBlocks==0 then this is a flag that a default set of tolerances should be used for all the linear systems. Default numBlocks=0 . | |
blockSolveCriteria | [in] Array (length numBlocks ) which gives the desired solution criteria for each of the numBlocks blocks of RHS. If numBlocks>0 then this argument must be non-NULL and point to a valid array of numBlocks entries. | |
blockSolveStatus | [out] Array (length numBlocks ) which gives the status of each set of block systems. A value of blockSolveStatus==NULL is allowed and means that the client is not interested in solution status of the linear systems. |
this->solveTransposeSupportsConj(conj)==true
X!=NULL
this->domain()->isCompatible(*B.range())==true
this->range()->isCompatible(*X->range())==true
B->domain()->isCompatible(*X->domain())==true
numBlocks >= 0
numBlocks == 0
] blockSolveCriteria==NULL && blockSolveStatus==NULL
. numBlocks > 0
] blockSolveCriteria!=NULL
and points to an array of length at least numBlocks
. numBlocks > 0
] this->solveTransposeSupportsSolveMeasureType(conj,blockSolveCriteria[k].solveMeasureType)==true
, for k = 0...numBlocks-1
. blockSolveStatus!=NULL
] blockSolveStatus
and points to an array of length at least numBlocks
. Postconditions:
CatastrophicSolveFailure
will be thrown. If the function returns normally, then the following postconditions apply. blockSolveStatus!=NULL
then blockSolveStatus[k]
gives the solution status of the block of linear systems specified by blockSolveCriteria[k]
where, for k=0...numBlocks-1
.
See the above introduction for a more complete description of how this function behaves and the meaning of the arguments blockSolveCriteria[]
and blockSolveStatus[]
.
Definition at line 62 of file Thyra_LinearOpWithSolveBase.hpp.