42 #ifndef BELOS_STATUS_TEST_GEN_IMPRESNORM_MP_VECTOR_HPP
43 #define BELOS_STATUS_TEST_GEN_IMPRESNORM_MP_VECTOR_HPP
54 #include "BelosStatusTestImpResNorm.hpp"
109 template <
class Storage,
class MV,
class OP>
111 public StatusTestResNorm<Sacado::MP::Vector<Storage>,MV,OP> {
116 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
121 typedef Teuchos::ScalarTraits<ScalarType>
STS;
122 typedef Teuchos::ScalarTraits<MagnitudeType>
STM;
123 typedef MultiVecTraits<ScalarType,MV>
MVT;
143 bool showMaxResNormOnly =
false);
189 tolerance_ = tolerance;
202 showMaxResNormOnly_ = showMaxResNormOnly;
235 void print(std::ostream& os,
int indent = 0)
const;
276 return currTolerance_;
280 const std::vector<MagnitudeType>*
getTestValue()
const {
return(&testvector_);}
313 std::ostringstream oss;
314 oss <<
"Belos::StatusTestImpResNorm<>: " << resFormStr();
315 oss <<
", tol = " << tolerance_;
329 std::ostringstream oss;
331 oss << ((resnormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
335 if (scaletype_!=
None)
341 if (scaletype_==UserProvided)
342 oss <<
" (User Scale)";
345 oss << ((scalenormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
346 if (scaletype_==NormOfInitRes)
348 else if (scaletype_==NormOfPrecInitRes)
435 template <
class StorageType,
class MV,
class OP>
438 : tolerance_(Tolerance),
439 currTolerance_(Tolerance),
441 showMaxResNormOnly_(showMaxResNormOnly),
442 resnormtype_(TwoNorm),
443 scaletype_(NormOfInitRes),
444 scalenormtype_(TwoNorm),
445 scalevalue_(
Teuchos::ScalarTraits<MagnitudeType>::one ()),
451 firstcallCheckStatus_(
true),
452 firstcallDefineResForm_(
true),
453 firstcallDefineScaleForm_(
true),
454 lossDetected_(false),
455 ensemble_iterations(StorageType::static_size, 0)
461 template <
class StorageType,
class MV,
class OP>
465 template <
class StorageType,
class MV,
class OP>
474 currTolerance_ = tolerance_;
475 firstcallCheckStatus_ =
true;
476 lossDetected_ =
false;
477 curSoln_ = Teuchos::null;
478 ensemble_iterations = std::vector<int>(StorageType::static_size, 0);
481 template <
class StorageType,
class MV,
class OP>
484 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==
false,StatusTestError,
485 "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
486 firstcallDefineResForm_ =
false;
488 resnormtype_ = TypeOfNorm;
493 template <
class StorageType,
class MV,
class OP>
495 MagnitudeType ScaleValue )
497 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==
false,StatusTestError,
498 "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
499 firstcallDefineScaleForm_ =
false;
501 scaletype_ = TypeOfScaling;
502 scalenormtype_ = TypeOfNorm;
503 scalevalue_ = ScaleValue;
508 template <
class StorageType,
class MV,
class OP>
510 checkStatus (Iteration<ScalarType,MV,OP>* iSolver)
516 const MagnitudeType zero = STM::zero ();
517 const MagnitudeType one = STM::one ();
518 const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem ();
521 if (firstcallCheckStatus_) {
522 StatusType status = firstCallCheckStatusSetup (iSolver);
523 if (status == Failed) {
532 if (curLSNum_ != lp.getLSNumber ()) {
536 curLSNum_ = lp.getLSNumber();
537 curLSIdx_ = lp.getLSIndex();
538 curBlksz_ = (int)curLSIdx_.size();
540 for (
int i=0; i<curBlksz_; ++i) {
541 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
544 curNumRHS_ = validLS;
545 curSoln_ = Teuchos::null;
550 if (status_ == Passed) {
575 std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
576 RCP<const MV> residMV = iSolver->getNativeResiduals (&tmp_resvector);
577 if (! residMV.is_null ()) {
579 tmp_resvector.resize (MVT::GetNumberVecs (*residMV));
580 MVT::MvNorm (*residMV, tmp_resvector, resnormtype_);
581 typename std::vector<int>::iterator p = curLSIdx_.begin();
582 for (
int i=0; p<curLSIdx_.end(); ++p, ++i) {
585 resvector_[*p] = tmp_resvector[i];
589 typename std::vector<int>::iterator p = curLSIdx_.begin();
590 for (
int i=0; p<curLSIdx_.end(); ++p, ++i) {
593 resvector_[*p] = tmp_resvector[i];
600 if (scalevector_.size () > 0) {
602 typename std::vector<int>::iterator p = curLSIdx_.begin();
603 for (; p<curLSIdx_.end(); ++p) {
607 testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
608 mask_assign(scalevector_[ *p ]!= zero, testvector_[ *p ]) /= scalevector_[ *p ];
613 typename std::vector<int>::iterator p = curLSIdx_.begin();
614 for (; p<curLSIdx_.end(); ++p) {
617 testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
630 ind_.resize( curLSIdx_.size() );
631 std::vector<int> lclInd( curLSIdx_.size() );
632 typename std::vector<int>::iterator p = curLSIdx_.begin();
633 for (
int i=0; p<curLSIdx_.end(); ++p, ++i) {
637 const int ensemble_size = StorageType::static_size;
638 bool all_converged =
true;
639 for (
int ii=0; ii<ensemble_size; ++ii) {
640 if (testvector_[ *p ].coeff(ii) > currTolerance_.coeff(ii)) {
641 ++ensemble_iterations[ii];
642 all_converged =
false;
644 else if (!(testvector_[ *p ].coeff(ii) <= currTolerance_.coeff(ii))) {
652 TEUCHOS_TEST_FOR_EXCEPTION(
true, StatusTestError,
"Belos::"
653 "StatusTestImpResNorm::checkStatus(): One or more of the current "
654 "implicit residual norms is NaN.");
672 RCP<MV> cur_update = iSolver->getCurrentUpdate ();
673 curSoln_ = lp.updateSolution (cur_update);
674 RCP<MV> cur_res = MVT::Clone (*curSoln_, MVT::GetNumberVecs (*curSoln_));
675 lp.computeCurrResVec (&*cur_res, &*curSoln_);
676 tmp_resvector.resize (MVT::GetNumberVecs (*cur_res));
677 std::vector<MagnitudeType> tmp_testvector (have);
679 MVT::MvNorm (*cur_res, tmp_resvector, resnormtype_);
682 if ( scalevector_.size() > 0 ) {
683 for (
int i=0; i<have; ++i) {
685 tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
686 mask_assign(scalevector_[ ind_[i] ]!=zero,tmp_testvector[ i ]) /= scalevector_[ ind_[i] ];
690 for (
int i=0; i<have; ++i) {
691 tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
702 for (
int i = 0; i < have; ++i) {
716 const MagnitudeType diff = STM::magnitude (testvector_[ind_[i]] - tmp_testvector[i]);
719 const int ensemble_size = StorageType::static_size;
721 bool all_converged =
true;
723 for (
int ii=0; ii<ensemble_size; ++ii) {
724 if (!(tmp_testvector[i].coeff(ii) <= tolerance_.coeff(ii))) {
725 if (diff.coeff(ii) > currTolerance_.coeff(ii)) {
726 lossDetected_ =
true;
729 const value_type onePointFive = as<value_type>(3) / as<value_type> (2);
730 const value_type oneTenth = as<value_type>(1) / as<value_type> (10);
732 currTolerance_.fastAccessCoeff(ii) = currTolerance_.coeff(ii) - onePointFive * diff.coeff(ii);
733 while (currTolerance_.coeff(ii) < as<value_type>(0)) {
734 currTolerance_.fastAccessCoeff(ii) += oneTenth * diff.coeff(ii);
736 all_converged =
false;
742 ind_[have2] = ind_[i];
752 int need = (quorum_ == -1) ? curNumRHS_: quorum_;
753 status_ = (have >= need) ? Passed : Failed;
761 template <
class StorageType,
class MV,
class OP>
764 for (
int j = 0;
j < indent;
j ++)
766 printStatus(os, status_);
768 if (status_==Undefined)
769 os <<
", tol = " << tolerance_ << std::endl;
772 if(showMaxResNormOnly_ && curBlksz_ > 1) {
773 const MagnitudeType maxRelRes = *std::max_element(
774 testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
776 for (
int j = 0;
j < indent + 13;
j ++)
778 os <<
"max{residual["<<curLSIdx_[0]<<
"..."<<curLSIdx_[curBlksz_-1]<<
"]} = " << maxRelRes
779 << ( maxRelRes <= tolerance_ ?
" <= " :
" > " ) << tolerance_ << std::endl;
782 for (
int i=0; i<numrhs_; i++ ) {
783 for (
int j = 0;
j < indent + 13;
j ++)
785 os <<
"residual [ " << i <<
" ] = " << testvector_[ i ];
786 os << ((testvector_[i]<tolerance_) ?
" < " : (testvector_[i]==tolerance_) ?
" == " : (testvector_[i]>tolerance_) ?
" > " :
" " ) << tolerance_ << std::endl;
793 template <
class StorageType,
class MV,
class OP>
796 os << std::left << std::setw(13) << std::setfill(
'.');
803 os <<
"Unconverged (LoA)";
812 os << std::left << std::setfill(
' ');
816 template <
class StorageType,
class MV,
class OP>
818 firstCallCheckStatusSetup (Iteration<ScalarType,MV,OP>* iSolver)
821 const MagnitudeType zero = STM::zero ();
822 const MagnitudeType one = STM::one ();
823 const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
825 if (firstcallCheckStatus_) {
829 firstcallCheckStatus_ =
false;
831 if (scaletype_== NormOfRHS) {
832 Teuchos::RCP<const MV> rhs = lp.getRHS();
833 numrhs_ = MVT::GetNumberVecs( *rhs );
834 scalevector_.resize( numrhs_ );
835 MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
837 else if (scaletype_==NormOfInitRes) {
838 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
839 numrhs_ = MVT::GetNumberVecs( *init_res );
840 scalevector_.resize( numrhs_ );
841 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
843 else if (scaletype_==NormOfPrecInitRes) {
844 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
845 numrhs_ = MVT::GetNumberVecs( *init_res );
846 scalevector_.resize( numrhs_ );
847 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
850 numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
853 resvector_.resize( numrhs_ );
854 testvector_.resize( numrhs_ );
856 curLSNum_ = lp.getLSNumber();
857 curLSIdx_ = lp.getLSIndex();
858 curBlksz_ = (int)curLSIdx_.size();
860 for (i=0; i<curBlksz_; ++i) {
861 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
864 curNumRHS_ = validLS;
867 for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
870 if (scalevalue_ == zero) {
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
NormType scalenormtype_
Type of norm to use on the scaling (OneNorm, TwoNorm, or InfNorm)
void reset()
Resets the internal configuration to the initial state.
bool lossDetected_
Has a loss in accuracy been detected?
int setShowMaxResNormOnly(bool showMaxResNormOnly)
Set whether the only maximum residual norm is displayed when the print() method is called.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
The type of the magnitude (absolute value) of a ScalarType.
NormType resnormtype_
Type of norm to use on residual (OneNorm, TwoNorm, or InfNorm).
bool getShowMaxResNormOnly()
Returns whether the only maximum residual norm is displayed when the print() method is called.
std::vector< int > ind_
Vector containing the indices for the vectors that passed the test.
MagnitudeType getTolerance() const
"Original" convergence tolerance as set by user.
Teuchos::RCP< MV > curSoln_
Current solution vector.
int quorum_
Number of residuals that must pass the convergence test before Passed is returned.
std::vector< int > ensemble_iterations
The number of iterations at which point each ensemble component converges.
int curNumRHS_
The current number of right-hand sides being solved for.
int defineResForm(NormType TypeOfNorm)
Define form of the residual, its norm and optional weighting vector.
void printStatus(std::ostream &os, StatusType type) const
Print message for each status specific to this stopping test.
Teuchos::ScalarTraits< MagnitudeType > STM
bool firstcallCheckStatus_
Is this the first time CheckStatus is called?
StatusType firstCallCheckStatusSetup(Iteration< ScalarType, MV, OP > *iSolver)
Call to setup initial scaling vector.
StatusType checkStatus(Iteration< ScalarType, MV, OP > *iSolver)
Check convergence status: Passed, Failed, or Undefined.
std::string description() const
Method to return description of the maximum iteration status test
std::vector< MagnitudeType > resvector_
Residual norm vector.
std::vector< MagnitudeType > scalevector_
Scaling vector.
std::vector< MagnitudeType > testvector_
Test vector = resvector_ / scalevector_.
const std::vector< MagnitudeType > * getResNormValue() const
Returns the residual norm value, , computed in most recent call to CheckStatus.
Teuchos::ScalarTraits< ScalarType > STS
StatusType status_
Status.
std::vector< int > convIndices()
Returns the vector containing the indices of the residuals that passed the test.
MultiVecTraits< ScalarType, MV > MVT
ScaleType scaletype_
Type of scaling to use (Norm of RHS, Norm of Initial Residual, None or User provided)
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
bool getLOADetected() const
Returns a boolean indicating a loss of accuracy has been detected in computing the residual.
void print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
std::string resFormStr() const
Description of current residual form.
StatusType getStatus() const
Return the result of the most recent CheckStatus call.
bool firstcallDefineScaleForm_
Is this the first time DefineScaleForm is called?
const std::vector< MagnitudeType > * getScaledNormValue() const
Returns the scaled norm value, .
MagnitudeType currTolerance_
MagnitudeType scalevalue_
Scaling value.
virtual ~StatusTestImpResNorm()
Destructor (virtual for memory safety).
bool showMaxResNormOnly_
Determines if the entries for all of the residuals are shown or just the max.
StatusTestImpResNorm(MagnitudeType Tolerance, int quorum=-1, bool showMaxResNormOnly=false)
Constructor.
int numrhs_
The total number of right-hand sides being solved for.
const std::vector< int > getEnsembleIterations() const
Returns number of ensemble iterations.
int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue=Teuchos::ScalarTraits< MagnitudeType >::one())
Define form of the scaling, its norm, its optional weighting vector, or, alternatively,...
MagnitudeType getCurrTolerance() const
Current convergence tolerance; may be changed to prevent loss of accuracy.
int setQuorum(int quorum)
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
int curBlksz_
The current blocksize of the linear system being solved.
bool firstcallDefineResForm_
Is this the first time DefineResForm is called?
Sacado::MP::Vector< Storage > ScalarType
std::vector< int > curLSIdx_
The indices of the current number of right-hand sides being solved for.
int curLSNum_
The current number of linear systems that have been loaded into the linear problem.
Teuchos::RCP< MV > getSolution()
Returns the current solution estimate that was computed for the most recent residual test.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...