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> 71 template <
class ScalarType,
class MV,
class OP>
76 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 "StatusTestGenResSubNorm::StatusTestGenResSubNorm(): StatusTestGenResSubNorm only available for blocked operators (e.g., Thyra).");
115 int defineResForm(
NormType TypeOfNorm) {
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& os,
int indent = 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 266 template <
class ScalarType>
267 class 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<ScalarType>::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<ScalarType>::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);
797 typedef MultiVecTraits<ScalarType, MV>
MVT;
802 void MvScalingRatio(
const MV& mv,
size_t block,
MagnitudeType& lengthRatio) {
803 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
805 typedef typename Thyra::ProductMultiVectorBase<ScalarType> TPMVB;
806 Teuchos::RCP<const TPMVB> thyProdVec = Teuchos::rcp_dynamic_cast<
const TPMVB>(input);
808 TEUCHOS_TEST_FOR_EXCEPTION(thyProdVec == Teuchos::null, std::invalid_argument,
809 "Belos::StatusTestGenResSubNorm::MvScalingRatio (Thyra specialization): " 810 "mv must be a Thyra::ProductMultiVector, but is of type " << thyProdVec);
812 Teuchos::RCP<const MV> thySubVec = thyProdVec->getMultiVectorBlock(block);
814 lengthRatio = Teuchos::as<MagnitudeType>(thySubVec->range()->dim()) / Teuchos::as<MagnitudeType>(thyProdVec->range()->dim());
832 bool showMaxResNormOnly_;
847 std::vector<MagnitudeType> scalevector_;
850 std::vector<MagnitudeType> resvector_;
853 std::vector<MagnitudeType> testvector_;
856 std::vector<int> ind_;
859 Teuchos::RCP<MV> curSoln_;
871 std::vector<int> curLSIdx_;
880 bool firstcallCheckStatus_;
883 bool firstcallDefineResForm_;
886 bool firstcallDefineScaleForm_;
892 #endif // HAVE_BELOS_THYRA
virtual std::vector< int > convIndices()=0
Returns the std::vector containing the indices of the residuals that passed the test.
ScaleType
The type of scaling to use on the residual norm value.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
virtual int setQuorum(int quorum)=0
Sets the number of residuals that must pass the convergence test before Passed is returned...
virtual int getQuorum() const =0
Returns the number of residuals that must pass the convergence test before Passed is returned...
An implementation of StatusTestResNorm using a family of norms of subvectors of the residual vectors...
Declaration of basic traits for the multivector type.
An abstract class of StatusTest for stopping criteria using residual norms.
virtual MagnitudeType getTolerance() const =0
Returns the value of the tolerance, , set in the constructor.
virtual StatusType checkStatus(Iteration< ScalarType, MV, OP > *iSolver)=0
Check convergence status: Unconverged, Converged, Failed.
Class which defines basic traits for the operator type.
StatusType
Whether the StatusTest wants iteration to stop.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
virtual void print(std::ostream &os, int indent=0) const =0
Output formatted description of stopping test to output stream.
virtual int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue=Teuchos::ScalarTraits< MagnitudeType >::one())=0
Define the form of the scaling for the residual.
MultiVecTraits< ScalarType, MV > MVT
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
virtual void printStatus(std::ostream &os, StatusType type) const
Output the result of the most recent CheckStatus call.
virtual Teuchos::RCP< MV > getSolution()=0
Returns the current solution estimate that was computed for the most recent residual test...
virtual int setShowMaxResNormOnly(bool showMaxResNormOnly)=0
Set whether the only maximum residual norm is displayed when the print() method is called...
virtual bool getShowMaxResNormOnly()=0
Returns whether the only maximum residual norm is displayed when the print() method is called...
Class which describes the linear problem to be solved by the iterative solver.
SCT::magnitudeType MagnitudeType
virtual const std::vector< MagnitudeType > * getTestValue() const =0
Returns the test value, , computed in most recent call to CheckStatus.
Belos::StatusTest abstract class for specifying a residual norm stopping criteria.
Teuchos::ScalarTraits< ScalarType > SCT
NormType
The type of vector norm to compute.
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 ...
virtual void reset()=0
Informs the convergence test that it should reset its internal configuration to the initialized state...
virtual bool getLOADetected() const =0
Returns a boolean indicating a loss of accuracy has been detected in computing the residual...
virtual StatusType getStatus() const =0
Return the result of the most recent CheckStatus call.
virtual int setTolerance(MagnitudeType tolerance)=0
Set the value of the tolerance.