47 #ifndef BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_MP_VECTORHPP
48 #define BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_MP_VECTORHPP
50 #include "Stokhos_config.h"
52 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
54 #include "Xpetra_ConfigDefs.hpp"
56 #include "Xpetra_BlockedCrsMatrix.hpp"
58 #include "MueLu_Exceptions.hpp"
60 #include <BelosConfigDefs.hpp>
61 #include <BelosTypes.hpp>
62 #include <BelosOperatorT.hpp>
63 #include <BelosXpetraAdapterOperator.hpp>
64 #include <BelosStatusTestGenResSubNorm.hpp>
65 #include <BelosXpetraStatusTestGenResSubNorm.hpp>
74 template <
class Storage,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 class StatusTestGenResSubNorm<typename Teuchos::ScalarTraits<Sacado::MP::Vector<Storage>>::magnitudeType, Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>>>
76 :
public StatusTestResNorm<typename Teuchos::ScalarTraits<Sacado::MP::Vector<Storage>>::magnitudeType, Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>>>
81 typedef Sacado::MP::Vector<Storage> ScalarType;
83 typedef Xpetra::MultiVector<ScalarType, LocalOrdinal, GlobalOrdinal, Node> MV;
84 typedef Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node> BCRS;
85 typedef Xpetra::MapExtractor<ScalarType, LocalOrdinal, GlobalOrdinal, Node> ME;
86 typedef Belos::OperatorT<MV> OP;
88 typedef Teuchos::ScalarTraits<ScalarType> SCT;
89 typedef typename SCT::magnitudeType MagnitudeType;
90 typedef MultiVecTraits<ScalarType, MV> MVT;
91 typedef OperatorTraits<ScalarType, MV, OP> OT;
108 StatusTestGenResSubNorm(MagnitudeType Tolerance,
size_t subIdx,
int quorum = -1,
bool showMaxResNormOnly =
false)
109 : tolerance_(Tolerance),
112 showMaxResNormOnly_(showMaxResNormOnly),
113 resnormtype_(TwoNorm),
114 scaletype_(NormOfInitRes),
115 scalenormtype_(TwoNorm),
116 scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one()),
122 firstcallCheckStatus_(true),
123 firstcallDefineResForm_(true),
124 firstcallDefineScaleForm_(true),
125 mapExtractor_(Teuchos::null),
129 virtual ~StatusTestGenResSubNorm(){};
142 int defineResForm(NormType TypeOfNorm)
144 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ ==
false, StatusTestError,
145 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
146 firstcallDefineResForm_ =
false;
148 resnormtype_ = TypeOfNorm;
174 int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one())
176 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ ==
false, StatusTestError,
177 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
178 firstcallDefineScaleForm_ =
false;
180 scaletype_ = TypeOfScaling;
181 scalenormtype_ = TypeOfNorm;
182 scalevalue_ = ScaleValue;
191 int setTolerance(MagnitudeType tolerance)
193 tolerance_ = tolerance;
200 int setSubIdx(
size_t subIdx)
208 int setQuorum(
int quorum)
215 int setShowMaxResNormOnly(
bool showMaxResNormOnly)
217 showMaxResNormOnly_ = showMaxResNormOnly;
232 StatusType checkStatus(Iteration<MagnitudeType, MV, OP> *iSolver)
235 const int ensemble_size = ET::size;
236 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
237 const LinearProblem<MagnitudeType, MV, OP> &lp = iSolver->getProblem();
239 if (firstcallCheckStatus_)
241 StatusType status = firstCallCheckStatusSetup(iSolver);
242 if (status == Failed)
252 if (curLSNum_ != lp.getLSNumber())
257 curLSNum_ = lp.getLSNumber();
258 curLSIdx_ = lp.getLSIndex();
259 curBlksz_ = (int)curLSIdx_.size();
261 for (
int i = 0; i < curBlksz_; ++i)
263 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
266 curNumRHS_ = validLS;
267 curSoln_ = Teuchos::null;
275 if (status_ == Passed)
284 Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
285 curSoln_ = lp.updateSolution(cur_update);
286 Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
287 lp.computeCurrResVec(&*cur_res, &*curSoln_);
288 std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
289 MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
291 typename std::vector<int>::iterator p = curLSIdx_.begin();
292 for (
int i = 0; p < curLSIdx_.end(); ++p, ++i)
296 resvector_[*p] = tmp_resvector[i];
303 if (scalevector_.size() > 0)
305 typename std::vector<int>::iterator pp = curLSIdx_.begin();
306 for (; pp < curLSIdx_.end(); ++pp)
312 if (scalevector_[*pp] != zero)
315 testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
319 testvector_[*pp] = resvector_[*pp] / scalevalue_;
326 typename std::vector<int>::iterator pp = curLSIdx_.begin();
327 for (; pp < curLSIdx_.end(); ++pp)
331 testvector_[*pp] = resvector_[*pp] / scalevalue_;
336 ind_.resize(curLSIdx_.size());
337 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
338 for (; p2 < curLSIdx_.end(); ++p2)
344 bool all_converged =
true;
345 for (
int ii = 0; ii < ensemble_size; ++ii)
347 if (ET::coeff(testvector_[*p2], ii) > ET::coeff(tolerance_, ii))
349 ++ensemble_iterations[ii];
350 all_converged =
false;
352 else if (!(ET::coeff(testvector_[*p2], ii) <= ET::coeff(tolerance_, ii)))
361 TEUCHOS_TEST_FOR_EXCEPTION(
true, StatusTestError,
"StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
372 int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
373 status_ = (have >= need) ? Passed : Failed;
379 StatusType getStatus()
const {
return (status_); };
394 firstcallCheckStatus_ =
true;
395 curSoln_ = Teuchos::null;
404 void print(std::ostream &os,
int indent = 0)
const
406 os.setf(std::ios_base::scientific);
407 for (
int j = 0; j < indent; j++)
409 printStatus(os, status_);
411 if (status_ == Undefined)
412 os <<
", tol = " << tolerance_ << std::endl;
416 if (showMaxResNormOnly_ && curBlksz_ > 1)
418 const MagnitudeType maxRelRes = *std::max_element(
419 testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
420 for (
int j = 0; j < indent + 13; j++)
422 os <<
"max{residual[" << curLSIdx_[0] <<
"..." << curLSIdx_[curBlksz_ - 1] <<
"]} = " << maxRelRes
423 << (maxRelRes <= tolerance_ ?
" <= " :
" > ") << tolerance_ << std::endl;
427 for (
int i = 0; i < numrhs_; i++)
429 for (
int j = 0; j < indent + 13; j++)
431 os <<
"residual [ " << i <<
" ] = " << testvector_[i];
432 os << ((testvector_[i] < tolerance_) ?
" < " : (testvector_[i] == tolerance_) ?
" == " : (testvector_[i] > tolerance_) ?
" > " :
" ") << tolerance_ << std::endl;
440 void printStatus(std::ostream &os, StatusType type)
const
442 os << std::left << std::setw(13) << std::setfill(
'.');
456 os << std::left << std::setfill(
' ');
465 Teuchos::RCP<MV> getSolution() {
return curSoln_; }
469 int getQuorum()
const {
return quorum_; }
472 size_t getSubIdx()
const {
return subIdx_; }
475 bool getShowMaxResNormOnly() {
return showMaxResNormOnly_; }
478 std::vector<int> convIndices() {
return ind_; }
481 MagnitudeType getTolerance()
const {
return (tolerance_); };
484 const std::vector<MagnitudeType> *getTestValue()
const {
return (&testvector_); };
487 const std::vector<MagnitudeType> *getResNormValue()
const {
return (&resvector_); };
490 const std::vector<MagnitudeType> *getScaledNormValue()
const {
return (&scalevector_); };
494 bool getLOADetected()
const {
return false; }
496 const std::vector<int> getEnsembleIterations()
const {
return ensemble_iterations; }
508 StatusType firstCallCheckStatusSetup(Iteration<MagnitudeType, MV, OP> *iSolver)
511 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
512 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
513 const LinearProblem<MagnitudeType, MV, OP> &lp = iSolver->getProblem();
515 if (firstcallCheckStatus_)
520 firstcallCheckStatus_ =
false;
523 Teuchos::RCP<const OP> Op = lp.getOperator();
524 Teuchos::RCP<const Belos::XpetraOp<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xOp =
525 Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(Op);
526 TEUCHOS_TEST_FOR_EXCEPTION(xOp.is_null(), MueLu::Exceptions::BadCast,
"Bad cast from \'const Belos::OperatorT\' to \'const Belos::XpetraOp\'. The origin type is " <<
typeid(
const OP).name() <<
".");
527 Teuchos::RCP<const Xpetra::Operator<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xIntOp =
529 TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast,
"Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
530 Teuchos::RCP<const Xpetra::Matrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xMat =
531 Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(xIntOp);
532 TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError,
"Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
533 Teuchos::RCP<const Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(xMat);
534 TEUCHOS_TEST_FOR_EXCEPTION(bMat.is_null(), MueLu::Exceptions::BadCast,
"Bad cast from \'const Xpetra::Matrix\' to \'const Xpetra::BlockedCrsMatrix\'. The origin type is " <<
typeid(
const Xpetra::Matrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>).name() <<
". Note: you need a BlockedCrsMatrix object for the StatusTestGenResSubNorm to work!");
535 mapExtractor_ = bMat->getRangeMapExtractor();
536 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError,
"Could not extract map extractor from BlockedCrsMatrix. Error.");
537 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_->NumMaps() <= subIdx_, MueLu::Exceptions::RuntimeError,
"The multivector is only split into " << mapExtractor_->NumMaps() <<
" sub parts. Cannot access sub-block " << subIdx_ <<
".");
540 if (scaletype_ == NormOfRHS)
542 Teuchos::RCP<const MV> rhs = lp.getRHS();
543 numrhs_ = MVT::GetNumberVecs(*rhs);
544 scalevector_.resize(numrhs_);
545 MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
547 else if (scaletype_ == NormOfInitRes)
549 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
550 numrhs_ = MVT::GetNumberVecs(*init_res);
551 scalevector_.resize(numrhs_);
552 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
554 else if (scaletype_ == NormOfPrecInitRes)
556 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
557 numrhs_ = MVT::GetNumberVecs(*init_res);
558 scalevector_.resize(numrhs_);
559 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
561 else if (scaletype_ == NormOfFullInitRes)
563 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
564 numrhs_ = MVT::GetNumberVecs(*init_res);
565 scalevector_.resize(numrhs_);
566 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
569 else if (scaletype_ == NormOfFullPrecInitRes)
571 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
572 numrhs_ = MVT::GetNumberVecs(*init_res);
573 scalevector_.resize(numrhs_);
574 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
577 else if (scaletype_ == NormOfFullScaledInitRes)
579 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
580 numrhs_ = MVT::GetNumberVecs(*init_res);
581 scalevector_.resize(numrhs_);
582 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
583 MvScalingRatio(*init_res, subIdx_, scalevalue_);
585 else if (scaletype_ == NormOfFullScaledPrecInitRes)
587 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
588 numrhs_ = MVT::GetNumberVecs(*init_res);
589 scalevector_.resize(numrhs_);
590 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
591 MvScalingRatio(*init_res, subIdx_, scalevalue_);
595 numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
598 resvector_.resize(numrhs_);
599 testvector_.resize(numrhs_);
601 curLSNum_ = lp.getLSNumber();
602 curLSIdx_ = lp.getLSIndex();
603 curBlksz_ = (int)curLSIdx_.size();
605 for (i = 0; i < curBlksz_; ++i)
607 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
610 curNumRHS_ = validLS;
613 for (i = 0; i < numrhs_; i++)
615 testvector_[i] = one;
619 if (scalevalue_ == zero)
632 std::string description()
const
634 std::ostringstream oss;
635 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
636 oss <<
", tol = " << tolerance_;
646 std::string resFormStr()
const
648 std::ostringstream oss;
650 oss << ((resnormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm" :
"Inf-Norm");
652 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
655 if (scaletype_ != None)
661 if (scaletype_ == UserProvided)
662 oss <<
" (User Scale)";
666 oss << ((scalenormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm" :
"Inf-Norm");
667 if (scaletype_ == NormOfInitRes)
668 oss <<
" Res0 [" << subIdx_ <<
"]";
669 else if (scaletype_ == NormOfPrecInitRes)
670 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
671 else if (scaletype_ == NormOfFullInitRes)
672 oss <<
" Full Res0 [" << subIdx_ <<
"]";
673 else if (scaletype_ == NormOfFullPrecInitRes)
674 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
675 else if (scaletype_ == NormOfFullScaledInitRes)
676 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
677 else if (scaletype_ == NormOfFullScaledPrecInitRes)
678 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
680 oss <<
" RHS [" << subIdx_ <<
"]";
696 void MvSubNorm(
const MV &mv,
size_t block, std::vector<MagnitudeType> &normVec, NormType type = TwoNorm)
699 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
701 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
702 typedef MultiVecTraits<ScalarType, MV> MVT;
703 MVT::MvNorm(*SubVec, normVec, type);
707 void MvScalingRatio(
const MV &mv,
size_t block, MagnitudeType &lengthRatio)
709 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
711 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
713 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
721 MagnitudeType tolerance_;
730 bool showMaxResNormOnly_;
733 NormType resnormtype_;
736 ScaleType scaletype_;
739 NormType scalenormtype_;
742 MagnitudeType scalevalue_;
745 std::vector<MagnitudeType> scalevector_;
748 std::vector<MagnitudeType> resvector_;
751 std::vector<MagnitudeType> testvector_;
754 std::vector<int> ind_;
757 Teuchos::RCP<MV> curSoln_;
769 std::vector<int> curLSIdx_;
778 bool firstcallCheckStatus_;
781 bool firstcallDefineResForm_;
784 bool firstcallDefineScaleForm_;
787 Teuchos::RCP<const ME> mapExtractor_;
790 std::vector<int> ensemble_iterations;
799 #include "Xpetra_ConfigDefs.hpp"
801 #include "Xpetra_BlockedCrsMatrix.hpp"
803 #include "MueLu_Exceptions.hpp"
805 #include <BelosConfigDefs.hpp>
806 #include <BelosTypes.hpp>
807 #include <BelosOperatorT.hpp>
808 #include <BelosXpetraAdapterOperator.hpp>
809 #include <BelosStatusTestGenResSubNorm.hpp>
818 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
819 class StatusTestGenResSubNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>>
820 :
public StatusTestResNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>>
825 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>
MV;
826 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>
BCRS;
827 typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>
ME;
828 typedef Belos::OperatorT<MV>
OP;
830 typedef Teuchos::ScalarTraits<Scalar>
SCT;
832 typedef MultiVecTraits<Scalar, MV>
MVT;
833 typedef OperatorTraits<Scalar, MV, OP>
OT;
848 StatusTestGenResSubNorm(
MagnitudeType Tolerance,
size_t subIdx,
int quorum = -1,
bool showMaxResNormOnly =
false)
849 : tolerance_(Tolerance),
852 showMaxResNormOnly_(showMaxResNormOnly),
853 resnormtype_(TwoNorm),
854 scaletype_(NormOfInitRes),
855 scalenormtype_(TwoNorm),
862 firstcallCheckStatus_(true),
863 firstcallDefineResForm_(true),
864 firstcallDefineScaleForm_(true),
865 mapExtractor_(Teuchos::null),
869 virtual ~StatusTestGenResSubNorm(){};
882 int defineResForm(NormType TypeOfNorm)
884 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ ==
false, StatusTestError,
885 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
886 firstcallDefineResForm_ =
false;
888 resnormtype_ = TypeOfNorm;
913 int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one())
915 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ ==
false, StatusTestError,
916 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
917 firstcallDefineScaleForm_ =
false;
919 scaletype_ = TypeOfScaling;
920 scalenormtype_ = TypeOfNorm;
921 scalevalue_ = ScaleValue;
930 int setTolerance(MagnitudeType tolerance)
932 tolerance_ = tolerance;
939 int setSubIdx(
size_t subIdx)
947 int setQuorum(
int quorum)
954 int setShowMaxResNormOnly(
bool showMaxResNormOnly)
956 showMaxResNormOnly_ = showMaxResNormOnly;
970 StatusType checkStatus(Iteration<Scalar, MV, OP> *iSolver)
973 const int ensemble_size = ET::size;
974 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
975 const LinearProblem<Scalar, MV, OP> &lp = iSolver->getProblem();
977 if (firstcallCheckStatus_)
979 StatusType status = firstCallCheckStatusSetup(iSolver);
980 if (status == Failed)
990 if (curLSNum_ != lp.getLSNumber())
995 curLSNum_ = lp.getLSNumber();
996 curLSIdx_ = lp.getLSIndex();
997 curBlksz_ = (int)curLSIdx_.size();
999 for (
int i = 0; i < curBlksz_; ++i)
1001 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
1004 curNumRHS_ = validLS;
1005 curSoln_ = Teuchos::null;
1013 if (status_ == Passed)
1022 Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
1023 curSoln_ = lp.updateSolution(cur_update);
1024 Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
1025 lp.computeCurrResVec(&*cur_res, &*curSoln_);
1026 std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
1027 MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
1029 typename std::vector<int>::iterator p = curLSIdx_.begin();
1030 for (
int i = 0; p < curLSIdx_.end(); ++p, ++i)
1034 resvector_[*p] = tmp_resvector[i];
1041 if (scalevector_.size() > 0)
1043 typename std::vector<int>::iterator pp = curLSIdx_.begin();
1044 for (; pp < curLSIdx_.end(); ++pp)
1050 if (scalevector_[*pp] != zero)
1053 testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
1057 testvector_[*pp] = resvector_[*pp] / scalevalue_;
1064 typename std::vector<int>::iterator pp = curLSIdx_.begin();
1065 for (; pp < curLSIdx_.end(); ++pp)
1069 testvector_[*pp] = resvector_[*pp] / scalevalue_;
1074 ind_.resize(curLSIdx_.size());
1075 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
1076 for (; p2 < curLSIdx_.end(); ++p2)
1082 bool all_converged =
true;
1083 for (
int ii = 0; ii < ensemble_size; ++ii)
1085 if (ET::coeff(testvector_[*p2], ii) > ET::coeff(tolerance_, ii))
1087 ++ensemble_iterations[ii];
1088 all_converged =
false;
1090 else if (!(ET::coeff(testvector_[*p2], ii) <= ET::coeff(tolerance_, ii)))
1099 TEUCHOS_TEST_FOR_EXCEPTION(
true, StatusTestError,
"StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
1110 int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
1111 status_ = (have >= need) ? Passed : Failed;
1117 StatusType getStatus()
const {
return (status_); };
1126 status_ = Undefined;
1129 curLSIdx_.resize(0);
1132 firstcallCheckStatus_ =
true;
1133 curSoln_ = Teuchos::null;
1142 void print(std::ostream &os,
int indent = 0)
const
1144 os.setf(std::ios_base::scientific);
1145 for (
int j = 0; j < indent; j++)
1147 printStatus(os, status_);
1149 if (status_ == Undefined)
1150 os <<
", tol = " << tolerance_ << std::endl;
1154 if (showMaxResNormOnly_ && curBlksz_ > 1)
1156 const MagnitudeType maxRelRes = *std::max_element(
1157 testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
1158 for (
int j = 0; j < indent + 13; j++)
1160 os <<
"max{residual[" << curLSIdx_[0] <<
"..." << curLSIdx_[curBlksz_ - 1] <<
"]} = " << maxRelRes
1161 << (maxRelRes <= tolerance_ ?
" <= " :
" > ") << tolerance_ << std::endl;
1165 for (
int i = 0; i < numrhs_; i++)
1167 for (
int j = 0; j < indent + 13; j++)
1169 os <<
"residual [ " << i <<
" ] = " << testvector_[i];
1170 os << ((testvector_[i] < tolerance_) ?
" < " : (testvector_[i] == tolerance_) ?
" == " : (testvector_[i] > tolerance_) ?
" > " :
" ") << tolerance_ << std::endl;
1178 void printStatus(std::ostream &os, StatusType type)
const
1180 os << std::left << std::setw(13) << std::setfill(
'.');
1187 os <<
"Unconverged";
1194 os << std::left << std::setfill(
' ');
1203 Teuchos::RCP<MV> getSolution() {
return curSoln_; }
1207 int getQuorum()
const {
return quorum_; }
1210 size_t getSubIdx()
const {
return subIdx_; }
1213 bool getShowMaxResNormOnly() {
return showMaxResNormOnly_; }
1216 std::vector<int> convIndices() {
return ind_; }
1219 MagnitudeType getTolerance()
const {
return (tolerance_); };
1222 const std::vector<MagnitudeType> *getTestValue()
const {
return (&testvector_); };
1225 const std::vector<MagnitudeType> *getResNormValue()
const {
return (&resvector_); };
1228 const std::vector<MagnitudeType> *getScaledNormValue()
const {
return (&scalevector_); };
1232 bool getLOADetected()
const {
return false; }
1234 const std::vector<int> getEnsembleIterations()
const {
return ensemble_iterations; }
1246 StatusType firstCallCheckStatusSetup(Iteration<Scalar, MV, OP> *iSolver)
1249 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
1250 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1251 const LinearProblem<Scalar, MV, OP> &lp = iSolver->getProblem();
1253 if (firstcallCheckStatus_)
1258 firstcallCheckStatus_ =
false;
1261 Teuchos::RCP<const OP> Op = lp.getOperator();
1262 Teuchos::RCP<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xOp =
1263 Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
1264 TEUCHOS_TEST_FOR_EXCEPTION(xOp.is_null(), MueLu::Exceptions::BadCast,
"Bad cast from \'const Belos::OperatorT\' to \'const Belos::XpetraOp\'. The origin type is " <<
typeid(
const OP).name() <<
".");
1265 Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xIntOp =
1267 TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast,
"Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
1268 Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xMat =
1269 Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xIntOp);
1270 TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError,
"Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
1271 Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xMat);
1272 TEUCHOS_TEST_FOR_EXCEPTION(bMat.is_null(), MueLu::Exceptions::BadCast,
"Bad cast from \'const Xpetra::Matrix\' to \'const Xpetra::BlockedCrsMatrix\'. The origin type is " <<
typeid(
const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>).name() <<
". Note: you need a BlockedCrsMatrix object for the StatusTestGenResSubNorm to work!");
1273 mapExtractor_ = bMat->getRangeMapExtractor();
1274 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError,
"Could not extract map extractor from BlockedCrsMatrix. Error.");
1275 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_->NumMaps() <= subIdx_, MueLu::Exceptions::RuntimeError,
"The multivector is only split into " << mapExtractor_->NumMaps() <<
" sub parts. Cannot access sub-block " << subIdx_ <<
".");
1278 if (scaletype_ == NormOfRHS)
1280 Teuchos::RCP<const MV> rhs = lp.getRHS();
1281 numrhs_ = MVT::GetNumberVecs(*rhs);
1282 scalevector_.resize(numrhs_);
1283 MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
1285 else if (scaletype_ == NormOfInitRes)
1287 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1288 numrhs_ = MVT::GetNumberVecs(*init_res);
1289 scalevector_.resize(numrhs_);
1290 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
1292 else if (scaletype_ == NormOfPrecInitRes)
1294 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1295 numrhs_ = MVT::GetNumberVecs(*init_res);
1296 scalevector_.resize(numrhs_);
1297 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
1299 else if (scaletype_ == NormOfFullInitRes)
1301 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1302 numrhs_ = MVT::GetNumberVecs(*init_res);
1303 scalevector_.resize(numrhs_);
1304 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1307 else if (scaletype_ == NormOfFullPrecInitRes)
1309 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1310 numrhs_ = MVT::GetNumberVecs(*init_res);
1311 scalevector_.resize(numrhs_);
1312 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1315 else if (scaletype_ == NormOfFullScaledInitRes)
1317 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1318 numrhs_ = MVT::GetNumberVecs(*init_res);
1319 scalevector_.resize(numrhs_);
1320 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1321 MvScalingRatio(*init_res, subIdx_, scalevalue_);
1323 else if (scaletype_ == NormOfFullScaledPrecInitRes)
1325 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1326 numrhs_ = MVT::GetNumberVecs(*init_res);
1327 scalevector_.resize(numrhs_);
1328 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1329 MvScalingRatio(*init_res, subIdx_, scalevalue_);
1333 numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
1336 resvector_.resize(numrhs_);
1337 testvector_.resize(numrhs_);
1339 curLSNum_ = lp.getLSNumber();
1340 curLSIdx_ = lp.getLSIndex();
1341 curBlksz_ = (int)curLSIdx_.size();
1343 for (i = 0; i < curBlksz_; ++i)
1345 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
1348 curNumRHS_ = validLS;
1351 for (i = 0; i < numrhs_; i++)
1353 testvector_[i] = one;
1357 if (scalevalue_ == zero)
1370 std::string description()
const
1372 std::ostringstream oss;
1373 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
1374 oss <<
", tol = " << tolerance_;
1384 std::string resFormStr()
const
1386 std::ostringstream oss;
1388 oss << ((resnormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm" :
"Inf-Norm");
1390 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
1393 if (scaletype_ != None)
1399 if (scaletype_ == UserProvided)
1400 oss <<
" (User Scale)";
1404 oss << ((scalenormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm" :
"Inf-Norm");
1405 if (scaletype_ == NormOfInitRes)
1406 oss <<
" Res0 [" << subIdx_ <<
"]";
1407 else if (scaletype_ == NormOfPrecInitRes)
1408 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
1409 else if (scaletype_ == NormOfFullInitRes)
1410 oss <<
" Full Res0 [" << subIdx_ <<
"]";
1411 else if (scaletype_ == NormOfFullPrecInitRes)
1412 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
1413 else if (scaletype_ == NormOfFullScaledInitRes)
1414 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
1415 else if (scaletype_ == NormOfFullScaledPrecInitRes)
1416 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
1418 oss <<
" RHS [" << subIdx_ <<
"]";
1434 void MvSubNorm(
const MV &mv,
size_t block, std::vector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normVec, NormType type = TwoNorm)
1437 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
1439 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
1440 typedef MultiVecTraits<Scalar, MV> MVT;
1441 MVT::MvNorm(*SubVec, normVec, type);
1445 void MvScalingRatio(
const MV &mv,
size_t block, MagnitudeType &lengthRatio)
1447 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
1449 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
1451 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
1459 MagnitudeType tolerance_;
1468 bool showMaxResNormOnly_;
1471 NormType resnormtype_;
1474 ScaleType scaletype_;
1477 NormType scalenormtype_;
1480 MagnitudeType scalevalue_;
1483 std::vector<MagnitudeType> scalevector_;
1486 std::vector<MagnitudeType> resvector_;
1489 std::vector<MagnitudeType> testvector_;
1492 std::vector<int> ind_;
1495 Teuchos::RCP<MV> curSoln_;
1507 std::vector<int> curLSIdx_;
1516 bool firstcallCheckStatus_;
1519 bool firstcallDefineResForm_;
1522 bool firstcallDefineScaleForm_;
1525 Teuchos::RCP<const ME> mapExtractor_;
1528 std::vector<int> ensemble_iterations;