42#ifndef BELOS_STATUS_TEST_GEN_IMPRESNORM_MP_VECTOR_HPP
43#define BELOS_STATUS_TEST_GEN_IMPRESNORM_MP_VECTOR_HPP
54#include "BelosStatusTestImpResNorm.hpp"
109template <
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)
435template <
class StorageType,
class MV,
class OP>
437StatusTestImpResNorm (MagnitudeType Tolerance,
int quorum,
bool showMaxResNormOnly)
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)
461template <
class StorageType,
class MV,
class OP>
465template <
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);
481template <
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;
493template <
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;
508template <
class StorageType,
class MV,
class OP>
510checkStatus (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;
720 typedef typename StorageType::value_type
value_type;
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;
761template <
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;
793template <
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(
' ');
816template <
class StorageType,
class MV,
class OP>
818firstCallCheckStatusSetup (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)
const std::vector< MagnitudeType > * getResNormValue() const
Returns the residual norm value, , computed in most recent call to CheckStatus.
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?
Teuchos::RCP< MV > getSolution()
Returns the current solution estimate that was computed for the most recent residual test.
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_.
Teuchos::ScalarTraits< ScalarType > STS
StatusType status_
Status.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
MultiVecTraits< ScalarType, MV > MVT
ScaleType scaletype_
Type of scaling to use (Norm of RHS, Norm of Initial Residual, None or User provided)
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< int > getEnsembleIterations() const
Returns number of ensemble iterations.
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.
std::vector< int > convIndices()
Returns the vector containing the indices of the residuals that passed the test.
int numrhs_
The total number of right-hand sides being solved for.
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.
const std::vector< MagnitudeType > * getScaledNormValue() const
Returns the scaled norm value, .
int curLSNum_
The current number of linear systems that have been loaded into the linear problem.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...