47 #ifndef BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP
48 #define BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP
56 #include <BelosConfigDefs.hpp>
57 #include <BelosTypes.hpp>
58 #include <BelosOperatorT.hpp>
59 #include <BelosXpetraAdapterOperator.hpp>
60 #include <BelosStatusTestGenResSubNorm.hpp>
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 class StatusTestGenResSubNorm<
Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >
69 :
public StatusTestResNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > {
75 typedef Belos::OperatorT<MV>
OP;
79 typedef MultiVecTraits<Scalar, MV>
MVT;
80 typedef OperatorTraits<Scalar, MV, OP>
OT;
97 StatusTestGenResSubNorm(
MagnitudeType Tolerance,
size_t subIdx,
int quorum = -1,
bool showMaxResNormOnly =
false)
98 : tolerance_(Tolerance)
101 , showMaxResNormOnly_(showMaxResNormOnly)
102 , resnormtype_(TwoNorm)
103 , scaletype_(NormOfInitRes)
104 , scalenormtype_(TwoNorm)
111 , firstcallCheckStatus_(true)
112 , firstcallDefineResForm_(true)
113 , firstcallDefineScaleForm_(true)
114 , mapExtractor_(Teuchos::null) {}
117 virtual ~StatusTestGenResSubNorm(){};
130 int defineResForm(NormType TypeOfNorm) {
132 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
133 firstcallDefineResForm_ =
false;
135 resnormtype_ = TypeOfNorm;
163 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
164 firstcallDefineScaleForm_ =
false;
166 scaletype_ = TypeOfScaling;
167 scalenormtype_ = TypeOfNorm;
168 scalevalue_ = ScaleValue;
177 int setTolerance(MagnitudeType tolerance) {
178 tolerance_ = tolerance;
185 int setSubIdx(
size_t subIdx) {
192 int setQuorum(
int quorum) {
198 int setShowMaxResNormOnly(
bool showMaxResNormOnly) {
199 showMaxResNormOnly_ = showMaxResNormOnly;
214 StatusType checkStatus(Iteration<Scalar, MV, OP>* iSolver) {
216 const LinearProblem<Scalar, MV, OP>& lp = iSolver->getProblem();
218 if (firstcallCheckStatus_) {
219 StatusType status = firstCallCheckStatusSetup(iSolver);
220 if (status == Failed) {
229 if (curLSNum_ != lp.getLSNumber()) {
233 curLSNum_ = lp.getLSNumber();
234 curLSIdx_ = lp.getLSIndex();
235 curBlksz_ = (int)curLSIdx_.size();
237 for (
int i = 0; i < curBlksz_; ++i) {
238 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
241 curNumRHS_ = validLS;
242 curSoln_ = Teuchos::null;
248 if (status_ == Passed) {
257 curSoln_ = lp.updateSolution(cur_update);
258 Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
259 lp.computeCurrResVec(&*cur_res, &*curSoln_);
260 std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
261 MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
263 typename std::vector<int>::iterator p = curLSIdx_.begin();
264 for (
int i = 0; p < curLSIdx_.end(); ++p, ++i) {
267 resvector_[*p] = tmp_resvector[i];
274 if (scalevector_.size() > 0) {
275 typename std::vector<int>::iterator pp = curLSIdx_.begin();
276 for (; pp < curLSIdx_.end(); ++pp) {
280 if (scalevector_[*pp] != zero) {
282 testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
284 testvector_[*pp] = resvector_[*pp] / scalevalue_;
289 typename std::vector<int>::iterator pp = curLSIdx_.begin();
290 for (; pp < curLSIdx_.end(); ++pp) {
293 testvector_[*pp] = resvector_[*pp] / scalevalue_;
298 ind_.resize(curLSIdx_.size());
299 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
300 for (; p2 < curLSIdx_.end(); ++p2) {
304 if (testvector_[*p2] > tolerance_) {
308 }
else if (testvector_[*p2] <= tolerance_) {
319 int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
320 status_ = (have >= need) ? Passed : Failed;
326 StatusType getStatus()
const {
return (status_); };
340 firstcallCheckStatus_ =
true;
341 curSoln_ = Teuchos::null;
350 void print(std::ostream& os,
int indent = 0)
const {
351 os.setf(std::ios_base::scientific);
352 for (
int j = 0; j < indent; j++)
354 printStatus(os, status_);
356 if (status_ == Undefined)
357 os <<
", tol = " << tolerance_ << std::endl;
360 if (showMaxResNormOnly_ && curBlksz_ > 1) {
361 const MagnitudeType maxRelRes = *std::max_element(
362 testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
363 for (
int j = 0; j < indent + 13; j++)
365 os <<
"max{residual[" << curLSIdx_[0] <<
"..." << curLSIdx_[curBlksz_ - 1] <<
"]} = " << maxRelRes
366 << (maxRelRes <= tolerance_ ?
" <= " :
" > ") << tolerance_ << std::endl;
368 for (
int i = 0; i < numrhs_; i++) {
369 for (
int j = 0; j < indent + 13; j++)
371 os <<
"residual [ " << i <<
" ] = " << testvector_[i];
372 os << ((testvector_[i] < tolerance_) ?
" < " : (testvector_[i] == tolerance_) ?
" == "
373 : (testvector_[i] > tolerance_) ?
" > "
375 << tolerance_ << std::endl;
383 void printStatus(std::ostream& os, StatusType type)
const {
384 os << std::left << std::setw(13) << std::setfill(
'.');
397 os << std::left << std::setfill(
' ');
410 int getQuorum()
const {
return quorum_; }
413 size_t getSubIdx()
const {
return subIdx_; }
416 bool getShowMaxResNormOnly() {
return showMaxResNormOnly_; }
419 std::vector<int> convIndices() {
return ind_; }
422 MagnitudeType getTolerance()
const {
return (tolerance_); };
425 const std::vector<MagnitudeType>* getTestValue()
const {
return (&testvector_); };
428 const std::vector<MagnitudeType>* getResNormValue()
const {
return (&resvector_); };
431 const std::vector<MagnitudeType>* getScaledNormValue()
const {
return (&scalevector_); };
435 bool getLOADetected()
const {
return false; }
447 StatusType firstCallCheckStatusSetup(Iteration<Scalar, MV, OP>* iSolver) {
451 const LinearProblem<Scalar, MV, OP>& lp = iSolver->getProblem();
453 if (firstcallCheckStatus_) {
457 firstcallCheckStatus_ =
false;
462 Teuchos::rcp_dynamic_cast<
const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Op);
472 mapExtractor_ = bMat->getRangeMapExtractor();
477 if (scaletype_ == NormOfRHS) {
479 numrhs_ = MVT::GetNumberVecs(*rhs);
480 scalevector_.resize(numrhs_);
481 MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
482 }
else if (scaletype_ == NormOfInitRes) {
484 numrhs_ = MVT::GetNumberVecs(*init_res);
485 scalevector_.resize(numrhs_);
486 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
487 }
else if (scaletype_ == NormOfPrecInitRes) {
489 numrhs_ = MVT::GetNumberVecs(*init_res);
490 scalevector_.resize(numrhs_);
491 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
492 }
else if (scaletype_ == NormOfFullInitRes) {
494 numrhs_ = MVT::GetNumberVecs(*init_res);
495 scalevector_.resize(numrhs_);
496 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
498 }
else if (scaletype_ == NormOfFullPrecInitRes) {
500 numrhs_ = MVT::GetNumberVecs(*init_res);
501 scalevector_.resize(numrhs_);
502 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
504 }
else if (scaletype_ == NormOfFullScaledInitRes) {
506 numrhs_ = MVT::GetNumberVecs(*init_res);
507 scalevector_.resize(numrhs_);
508 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
509 MvScalingRatio(*init_res, subIdx_, scalevalue_);
510 }
else if (scaletype_ == NormOfFullScaledPrecInitRes) {
512 numrhs_ = MVT::GetNumberVecs(*init_res);
513 scalevector_.resize(numrhs_);
514 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
515 MvScalingRatio(*init_res, subIdx_, scalevalue_);
517 numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
520 resvector_.resize(numrhs_);
521 testvector_.resize(numrhs_);
523 curLSNum_ = lp.getLSNumber();
524 curLSIdx_ = lp.getLSIndex();
525 curBlksz_ = (int)curLSIdx_.size();
527 for (i = 0; i < curBlksz_; ++i) {
528 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
531 curNumRHS_ = validLS;
534 for (i = 0; i < numrhs_; i++) {
535 testvector_[i] = one;
539 if (scalevalue_ == zero) {
551 std::string description()
const {
552 std::ostringstream oss;
553 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
554 oss <<
", tol = " << tolerance_;
564 std::string resFormStr()
const {
565 std::ostringstream oss;
567 oss << ((resnormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm"
570 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
573 if (scaletype_ !=
None) {
578 if (scaletype_ == UserProvided)
579 oss <<
" (User Scale)";
582 oss << ((scalenormtype_ == OneNorm) ?
"1-Norm" : (resnormtype_ == TwoNorm) ?
"2-Norm"
584 if (scaletype_ == NormOfInitRes)
585 oss <<
" Res0 [" << subIdx_ <<
"]";
586 else if (scaletype_ == NormOfPrecInitRes)
587 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
588 else if (scaletype_ == NormOfFullInitRes)
589 oss <<
" Full Res0 [" << subIdx_ <<
"]";
590 else if (scaletype_ == NormOfFullPrecInitRes)
591 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
592 else if (scaletype_ == NormOfFullScaledInitRes)
593 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
594 else if (scaletype_ == NormOfFullScaledPrecInitRes)
595 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
597 oss <<
" RHS [" << subIdx_ <<
"]";
617 MVT::MvNorm(*SubVec, normVec, type);
621 void MvScalingRatio(
const MV& mv,
size_t block, MagnitudeType& lengthRatio) {
626 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
634 MagnitudeType tolerance_;
643 bool showMaxResNormOnly_;
646 NormType resnormtype_;
649 ScaleType scaletype_;
652 NormType scalenormtype_;
655 MagnitudeType scalevalue_;
658 std::vector<MagnitudeType> scalevector_;
661 std::vector<MagnitudeType> resvector_;
664 std::vector<MagnitudeType> testvector_;
667 std::vector<int> ind_;
682 std::vector<int> curLSIdx_;
691 bool firstcallCheckStatus_;
694 bool firstcallDefineResForm_;
697 bool firstcallDefineScaleForm_;
Exception indicating invalid cast attempted.
SCT::magnitudeType MagnitudeType
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MultiVecTraits< Scalar, MV > MVT
MueLu::DefaultScalar Scalar
OperatorTraits< Scalar, MV, OP > OT
Teuchos::ScalarTraits< Scalar > SCT
Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > BCRS
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > ME
Exception throws to report errors in the internal logical of the program.
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MV
Belos::OperatorT< MV > OP