42#ifndef BELOS_STATUS_TEST_GEN_RESSUBNORM_H
43#define BELOS_STATUS_TEST_GEN_RESSUBNORM_H
55#ifdef HAVE_BELOS_THYRA
56#include <Thyra_MultiVectorBase.hpp>
57#include <Thyra_MultiVectorStdOps.hpp>
58#include <Thyra_ProductMultiVectorBase.hpp>
71template <
class ScalarType,
class MV,
class OP>
76 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 "StatusTestGenResSubNorm::StatusTestGenResSubNorm(): StatusTestGenResSubNorm only available for blocked operators (e.g., Thyra).");
116 TEUCHOS_TEST_FOR_EXCEPTION(
true,StatusTestError,
117 "StatusTestGenResSubNorm::defineResForm(): StatusTestGenResSubNorm only available for blocked operators (e.g., Thyra).");
118 TEUCHOS_UNREACHABLE_RETURN(0);
143 TEUCHOS_TEST_FOR_EXCEPTION(
true,StatusTestError,
144 "StatusTestGenResSubNorm::defineScaleForm(): StatusTestGenResSubNorm only available for blocked operators (e.g., Thyra).");
145 TEUCHOS_UNREACHABLE_RETURN(0);
157 int setSubIdx (
size_t subIdx ) {
return 0;}
195 void print(std::ostream& ,
int = 0)
const { }
207 Teuchos::RCP<MV>
getSolution() {
return Teuchos::null; }
214 size_t getSubIdx()
const {
return 0; }
220 std::vector<int>
convIndices() {
return std::vector<int>(0); }
226 const std::vector<MagnitudeType>*
getTestValue()
const {
return NULL;};
229 const std::vector<MagnitudeType>* getResNormValue()
const {
return NULL;};
232 const std::vector<MagnitudeType>* getScaledNormValue()
const {
return NULL;};
249 StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver) {
258 std::string description()
const
259 {
return std::string(
""); }
263#ifdef HAVE_BELOS_THYRA
266template <
class ScalarType>
267class StatusTestGenResSubNorm<ScalarType,Thyra::MultiVectorBase<ScalarType>,Thyra::LinearOpBase<ScalarType> >
268 :
public StatusTestResNorm<ScalarType,Thyra::MultiVectorBase<ScalarType>,Thyra::LinearOpBase<ScalarType> > {
272 typedef Thyra::MultiVectorBase<ScalarType> MV;
273 typedef Thyra::LinearOpBase<ScalarType> OP;
275 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
277 typedef MultiVecTraits<ScalarType,MV>
MVT;
278 typedef OperatorTraits<ScalarType,MV,OP> OT;
295 StatusTestGenResSubNorm(
MagnitudeType Tolerance,
size_t subIdx,
int quorum = -1,
bool showMaxResNormOnly =
false )
296 : tolerance_(Tolerance),
299 showMaxResNormOnly_(showMaxResNormOnly),
309 firstcallCheckStatus_(true),
310 firstcallDefineResForm_(true),
311 firstcallDefineScaleForm_(true) { }
314 virtual ~StatusTestGenResSubNorm() { };
327 int defineResForm(
NormType TypeOfNorm) {
328 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==
false,StatusTestError,
329 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
330 firstcallDefineResForm_ =
false;
332 resnormtype_ = TypeOfNorm;
359 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==
false,StatusTestError,
360 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
361 firstcallDefineScaleForm_ =
false;
363 scaletype_ = TypeOfScaling;
364 scalenormtype_ = TypeOfNorm;
365 scalevalue_ = ScaleValue;
379 int setSubIdx (
size_t subIdx ) { subIdx_ = subIdx;
return(0);}
383 int setQuorum(
int quorum) {quorum_ = quorum;
return(0);}
386 int setShowMaxResNormOnly(
bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly;
return(0);}
400 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
401 const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
403 if (firstcallCheckStatus_) {
404 StatusType status = firstCallCheckStatusSetup(iSolver);
414 if ( curLSNum_ != lp.getLSNumber() ) {
418 curLSNum_ = lp.getLSNumber();
419 curLSIdx_ = lp.getLSIndex();
420 curBlksz_ = (int)curLSIdx_.size();
422 for (
int i=0; i<curBlksz_; ++i) {
423 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
426 curNumRHS_ = validLS;
427 curSoln_ = Teuchos::null;
433 if (status_==
Passed) {
return status_; }
439 Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
440 curSoln_ = lp.updateSolution( cur_update );
442 lp.computeCurrResVec( &*cur_res, &*curSoln_ );
444 MvSubNorm( *cur_res, subIdx_, tmp_resvector, resnormtype_ );
446 typename std::vector<int>::iterator pp = curLSIdx_.begin();
447 for (
int i=0; pp<curLSIdx_.end(); ++pp, ++i) {
450 resvector_[*pp] = tmp_resvector[i];
457 if ( scalevector_.size() > 0 ) {
458 typename std::vector<int>::iterator p = curLSIdx_.begin();
459 for (; p<curLSIdx_.end(); ++p) {
463 if ( scalevector_[ *p ] != zero ) {
465 testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
467 testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
473 typename std::vector<int>::iterator ppp = curLSIdx_.begin();
474 for (; ppp<curLSIdx_.end(); ++ppp) {
477 testvector_[ *ppp ] = resvector_[ *ppp ] / scalevalue_;
482 ind_.resize( curLSIdx_.size() );
483 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
484 for (; p2<curLSIdx_.end(); ++p2) {
488 if (testvector_[ *p2 ] > tolerance_) {
490 }
else if (testvector_[ *p2 ] <= tolerance_) {
496 TEUCHOS_TEST_FOR_EXCEPTION(
true,StatusTestError,
"StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
501 int need = (quorum_ == -1) ? curNumRHS_: quorum_;
522 firstcallCheckStatus_ =
true;
523 curSoln_ = Teuchos::null;
532 void print(std::ostream& os,
int indent = 0)
const {
533 os.setf(std::ios_base::scientific);
534 for (
int j = 0; j < indent; j ++)
539 os <<
", tol = " << tolerance_ << std::endl;
542 if(showMaxResNormOnly_ && curBlksz_ > 1) {
544 testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
546 for (
int j = 0; j < indent + 13; j ++)
548 os <<
"max{residual["<<curLSIdx_[0]<<
"..."<<curLSIdx_[curBlksz_-1]<<
"]} = " << maxRelRes
549 << ( maxRelRes <= tolerance_ ?
" <= " :
" > " ) << tolerance_ << std::endl;
552 for (
int i=0; i<numrhs_; i++ ) {
553 for (
int j = 0; j < indent + 13; j ++)
555 os <<
"residual [ " << i <<
" ] = " << testvector_[ i ];
556 os << ((testvector_[i]<tolerance_) ?
" < " : (testvector_[i]==tolerance_) ?
" == " : (testvector_[i]>tolerance_) ?
" > " :
" " ) << tolerance_ << std::endl;
565 os << std::left << std::setw(13) << std::setfill(
'.');
578 os << std::left << std::setfill(
' ');
587 Teuchos::RCP<MV>
getSolution() {
return curSoln_; }
591 int getQuorum()
const {
return quorum_; }
594 size_t getSubIdx()
const {
return subIdx_; }
606 const std::vector<MagnitudeType>*
getTestValue()
const {
return(&testvector_);};
609 const std::vector<MagnitudeType>* getResNormValue()
const {
return(&resvector_);};
612 const std::vector<MagnitudeType>* getScaledNormValue()
const {
return(&scalevector_);};
629 StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver) {
631 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
632 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
633 const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
635 if (firstcallCheckStatus_) {
639 firstcallCheckStatus_ =
false;
642 Teuchos::RCP<const MV> rhs = lp.getRHS();
644 scalevector_.resize( numrhs_ );
645 MvSubNorm( *rhs, subIdx_, scalevector_, scalenormtype_ );
648 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
650 scalevector_.resize( numrhs_ );
651 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
654 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
656 scalevector_.resize( numrhs_ );
657 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
660 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
662 scalevector_.resize( numrhs_ );
663 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
664 scalevalue_ = Teuchos::ScalarTraits<MagnitudeType>::one();
667 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
669 scalevector_.resize( numrhs_ );
670 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
671 scalevalue_ = Teuchos::ScalarTraits<MagnitudeType>::one();
674 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
676 scalevector_.resize( numrhs_ );
677 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
678 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
681 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
683 scalevector_.resize( numrhs_ );
684 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
685 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
691 resvector_.resize( numrhs_ );
692 testvector_.resize( numrhs_ );
694 curLSNum_ = lp.getLSNumber();
695 curLSIdx_ = lp.getLSIndex();
696 curBlksz_ = (int)curLSIdx_.size();
698 for (i=0; i<curBlksz_; ++i) {
699 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
702 curNumRHS_ = validLS;
705 for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
708 if (scalevalue_ == zero) {
720 std::string description()
const
722 std::ostringstream oss;
723 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
724 oss <<
", tol = " << tolerance_;
736 std::string resFormStr()
const
738 std::ostringstream oss;
740 oss << ((resnormtype_==
OneNorm) ?
"1-Norm" : (resnormtype_==
TwoNorm) ?
"2-Norm" :
"Inf-Norm");
742 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
745 if (scaletype_!=
None)
752 oss <<
" (User Scale)";
755 oss << ((scalenormtype_==
OneNorm) ?
"1-Norm" : (resnormtype_==
TwoNorm) ?
"2-Norm" :
"Inf-Norm");
757 oss <<
" Res0 [" << subIdx_ <<
"]";
759 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
761 oss <<
" Full Res0 [" << subIdx_ <<
"]";
763 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
765 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
767 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
769 oss <<
" RHS [" << subIdx_ <<
"]";
785 void MvSubNorm(
const MV& mv,
size_t block, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normVec,
NormType type =
TwoNorm) {
786 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
788 typedef typename Thyra::ProductMultiVectorBase<ScalarType> TPMVB;
789 Teuchos::RCP<const TPMVB> thyProdVec = Teuchos::rcp_dynamic_cast<const TPMVB>(input);
791 TEUCHOS_TEST_FOR_EXCEPTION(thyProdVec == Teuchos::null, std::invalid_argument,
792 "Belos::StatusTestGenResSubNorm::MvSubNorm (Thyra specialization): "
793 "mv must be a Thyra::ProductMultiVector, but is of type " << thyProdVec);
795 Teuchos::RCP<const MV> thySubVec = thyProdVec->getMultiVectorBlock(block);
801 void MvScalingRatio(
const MV& mv,
size_t block,
MagnitudeType& lengthRatio) {
802 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
804 typedef typename Thyra::ProductMultiVectorBase<ScalarType> TPMVB;
805 Teuchos::RCP<const TPMVB> thyProdVec = Teuchos::rcp_dynamic_cast<const TPMVB>(input);
807 TEUCHOS_TEST_FOR_EXCEPTION(thyProdVec == Teuchos::null, std::invalid_argument,
808 "Belos::StatusTestGenResSubNorm::MvScalingRatio (Thyra specialization): "
809 "mv must be a Thyra::ProductMultiVector, but is of type " << thyProdVec);
811 Teuchos::RCP<const MV> thySubVec = thyProdVec->getMultiVectorBlock(block);
813 lengthRatio = Teuchos::as<MagnitudeType>(thySubVec->range()->dim()) / Teuchos::as<MagnitudeType>(thyProdVec->range()->dim());
831 bool showMaxResNormOnly_;
846 std::vector<MagnitudeType> scalevector_;
849 std::vector<MagnitudeType> resvector_;
852 std::vector<MagnitudeType> testvector_;
855 std::vector<int> ind_;
858 Teuchos::RCP<MV> curSoln_;
870 std::vector<int> curLSIdx_;
879 bool firstcallCheckStatus_;
882 bool firstcallDefineResForm_;
885 bool firstcallDefineScaleForm_;
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Belos::StatusTest abstract class for specifying a residual norm stopping criteria.
Traits class which defines basic operations on multivectors.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ,...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
An implementation of StatusTestResNorm using a family of norms of subvectors of the residual vectors.
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
MultiVecTraits< ScalarType, MV > MVT
An abstract class of StatusTest for stopping criteria using residual norms.
virtual const std::vector< MagnitudeType > * getTestValue() const =0
Returns the test value, , computed in most recent call to CheckStatus.
virtual int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue=Teuchos::ScalarTraits< MagnitudeType >::one())=0
Define the form of the scaling for the residual.
virtual MagnitudeType getTolerance() const =0
Returns the value of the tolerance, , set in the constructor.
virtual bool getShowMaxResNormOnly()=0
Returns whether the only maximum residual norm is displayed when the print() method is called.
virtual bool getLOADetected() const =0
Returns a boolean indicating a loss of accuracy has been detected in computing the residual.
virtual int getQuorum() const =0
Returns the number of residuals that must pass the convergence test before Passed is returned.
virtual int setQuorum(int quorum)=0
Sets the number of residuals that must pass the convergence test before Passed is returned.
virtual int setShowMaxResNormOnly(bool showMaxResNormOnly)=0
Set whether the only maximum residual norm is displayed when the print() method is called.
virtual int setTolerance(MagnitudeType tolerance)=0
Set the value of the tolerance.
virtual Teuchos::RCP< MV > getSolution()=0
Returns the current solution estimate that was computed for the most recent residual test.
virtual std::vector< int > convIndices()=0
Returns the std::vector containing the indices of the residuals that passed the test.
virtual StatusType getStatus() const =0
Return the result of the most recent CheckStatus call.
virtual void print(std::ostream &os, int indent=0) const =0
Output formatted description of stopping test to output stream.
virtual StatusType checkStatus(Iteration< ScalarType, MV, OP > *iSolver)=0
Check convergence status: Unconverged, Converged, Failed.
virtual void reset()=0
Informs the convergence test that it should reset its internal configuration to the initialized state...
virtual void printStatus(std::ostream &os, StatusType type) const
Output the result of the most recent CheckStatus call.
NormType
The type of vector norm to compute.
StatusType
Whether the StatusTest wants iteration to stop.
ScaleType
The type of scaling to use on the residual norm value.
@ NormOfFullScaledPrecInitRes
@ NormOfFullScaledInitRes