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);
159 int defineResForm( NormType TypeOfNorm);
182 int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm,
MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
189 tolerance_ = tolerance;
202 showMaxResNormOnly_ = showMaxResNormOnly;
217 StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
235 void print(std::ostream& os,
int indent = 0)
const;
238 void printStatus(std::ostream& os, StatusType type)
const;
276 return currTolerance_;
280 const std::vector<MagnitudeType>*
getTestValue()
const {
return(&testvector_);}
304 StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
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>
462 StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::~StatusTestImpResNorm()
465 template <
class StorageType,
class MV,
class OP>
466 void StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::reset()
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>
482 int StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineResForm( NormType TypeOfNorm )
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>
494 int StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
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>
509 StatusType StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,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>
762 void StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::print(std::ostream& os,
int indent)
const 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>
794 void StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::printStatus(std::ostream& os, StatusType type)
const 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>
817 StatusType StatusTestImpResNorm<Sacado::MP::Vector<StorageType>,MV,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) {
NormType scalenormtype_
Type of norm to use on the scaling (OneNorm, TwoNorm, or InfNorm)
MagnitudeType getTolerance() const
"Original" convergence tolerance as set by user.
MultiVecTraits< ScalarType, MV > MVT
std::vector< MagnitudeType > testvector_
Test vector = resvector_ / scalevector_.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
MagnitudeType tolerance_
Current tolerance used to determine convergence and the default tolerance set by user.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
int curBlksz_
The current blocksize of the linear system being solved.
Sacado::MP::Vector< Storage > ScalarType
const std::vector< MagnitudeType > * getScaledNormValue() const
Returns the scaled norm value, .
int numrhs_
The total number of right-hand sides being solved for.
int curNumRHS_
The current number of right-hand sides being solved for.
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
std::vector< MagnitudeType > scalevector_
Scaling vector.
std::vector< int > convIndices()
Returns the vector containing the indices of the residuals that passed the test.
bool firstcallDefineScaleForm_
Is this the first time DefineScaleForm is called?
int setQuorum(int quorum)
std::vector< int > curLSIdx_
The indices of the current number of right-hand sides being solved for.
bool getLOADetected() const
Returns a boolean indicating a loss of accuracy has been detected in computing the residual...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
The type of the magnitude (absolute value) of a ScalarType.
MagnitudeType scalevalue_
Scaling value.
Teuchos::ScalarTraits< MagnitudeType > STM
std::string resFormStr() const
Description of current residual form.
int curLSNum_
The current number of linear systems that have been loaded into the linear problem.
Teuchos::RCP< MV > curSoln_
Current solution vector.
const std::vector< MagnitudeType > * getResNormValue() const
Returns the residual norm value, , computed in most recent call to CheckStatus.
int setShowMaxResNormOnly(bool showMaxResNormOnly)
Set whether the only maximum residual norm is displayed when the print() method is called...
MagnitudeType getCurrTolerance() const
Current convergence tolerance; may be changed to prevent loss of accuracy.
bool showMaxResNormOnly_
Determines if the entries for all of the residuals are shown or just the max.
bool firstcallCheckStatus_
Is this the first time CheckStatus is called?
Teuchos::ScalarTraits< ScalarType > STS
bool lossDetected_
Has a loss in accuracy been detected?
StatusType status_
Status.
StatusType getStatus() const
Return the result of the most recent CheckStatus call.
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.
ScaleType scaletype_
Type of scaling to use (Norm of RHS, Norm of Initial Residual, None or User provided) ...
bool firstcallDefineResForm_
Is this the first time DefineResForm is called?
Teuchos::RCP< MV > getSolution()
Returns the current solution estimate that was computed for the most recent residual test...
const std::vector< int > getEnsembleIterations() const
Returns number of ensemble iterations.
std::vector< int > ind_
Vector containing the indices for the vectors that passed the test.
NormType resnormtype_
Type of norm to use on residual (OneNorm, TwoNorm, or InfNorm).
std::string description() const
Method to return description of the maximum iteration status test.
std::vector< MagnitudeType > resvector_
Residual norm vector.
bool getShowMaxResNormOnly()
Returns whether the only maximum residual norm is displayed when the print() method is called...