42 #ifndef BELOS_STATUS_TEST_WEIGHTED_GEN_RESNORM_MP_VECTORH
43 #define BELOS_STATUS_TEST_WEIGHTED_GEN_RESNORM_MP_VECTORH
50 #include "BelosStatusTestResNorm.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosMultiVecTraits.hpp"
83 template <
class ScalarType,
class MV,
class OP>
89 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
91 typedef MultiVecTraits<ScalarType, MV>
MVT;
164 int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
204 StatusType
checkStatus(Iteration<ScalarType, MV, OP> *iSolver);
222 void print(std::ostream &os,
int indent = 0)
const;
225 void printStatus(std::ostream &os, StatusType type)
const;
283 std::ostringstream oss;
284 oss <<
"Belos::StatusTestWeightedGenResNorm<>: " <<
resFormStr();
297 std::ostringstream oss;
311 oss <<
" (User Scale)";
335 void MvWNorm(
const MV &mv,
const MV &w, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normVec, NormType type = TwoNorm)
338 Teuchos::RCP<MV> wmv = MVT::CloneCopy(mv);
340 Teuchos::ArrayRCP<typename MV::scalar_type> wmv_entries;
341 Teuchos::ArrayRCP<const typename MV::scalar_type> w_entries;
343 for (
size_t j = 0; j < MVT::GetNumberVecs(mv); ++j)
345 wmv_entries = wmv->getDataNonConst(j);
346 w_entries = w.getData(j);
347 for (
size_t i = 0; i < wmv_entries.size(); ++i)
348 wmv_entries[i] *= w_entries[i];
350 MVT::MvNorm(*wmv, normVec, type);
429 template <
class ScalarType,
class MV,
class OP>
432 : tolerance_(Tolerance),
434 showMaxResNormOnly_(showMaxResNormOnly),
435 resnormtype_(TwoNorm),
436 scaletype_(NormOfInitRes),
437 scalenormtype_(TwoNorm),
444 firstcallCheckStatus_(true),
445 firstcallDefineResForm_(true),
446 firstcallDefineScaleForm_(true),
454 template <
class ScalarType,
class MV,
class OP>
459 template <
class ScalarType,
class MV,
class OP>
468 firstcallCheckStatus_ =
true;
469 curSoln_ = Teuchos::null;
472 template <
class ScalarType,
class MV,
class OP>
475 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ ==
false, StatusTestError,
476 "StatusTestWeightedGenResNorm::defineResForm(): The residual form has already been defined.");
477 firstcallDefineResForm_ =
false;
479 resnormtype_ = TypeOfNorm;
484 template <
class ScalarType,
class MV,
class OP>
488 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ ==
false, StatusTestError,
489 "StatusTestWeightedGenResNorm::defineScaleForm(): The scaling type has already been defined.");
490 firstcallDefineScaleForm_ =
false;
492 scaletype_ = TypeOfScaling;
493 scalenormtype_ = TypeOfNorm;
494 scalevalue_ = ScaleValue;
499 template <
class ScalarType,
class MV,
class OP>
503 const int ensemble_size = ET::size;
504 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
505 const LinearProblem<ScalarType, MV, OP> &lp = iSolver->getProblem();
507 if (firstcallCheckStatus_)
509 StatusType status = firstCallCheckStatusSetup(iSolver);
510 if (status == Failed)
519 if (curLSNum_ != lp.getLSNumber())
524 curLSNum_ = lp.getLSNumber();
525 curLSIdx_ = lp.getLSIndex();
526 curBlksz_ = (int)curLSIdx_.size();
528 for (
int i = 0; i < curBlksz_; ++i)
530 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
533 curNumRHS_ = validLS;
534 curSoln_ = Teuchos::null;
542 if (status_ == Passed)
551 Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
552 curSoln_ = lp.updateSolution(cur_update);
553 Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
554 lp.computeCurrResVec(&*cur_res, &*curSoln_);
555 std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
556 MvWNorm(*cur_res, *weights_, tmp_resvector, resnormtype_);
557 typename std::vector<int>::iterator p = curLSIdx_.begin();
558 for (
int i = 0; p < curLSIdx_.end(); ++p, ++i)
562 resvector_[*p] = tmp_resvector[i];
569 if (scalevector_.size() > 0)
571 typename std::vector<int>::iterator p = curLSIdx_.begin();
572 for (; p < curLSIdx_.end(); ++p)
578 if (scalevector_[*p] != zero)
581 testvector_[*p] = resvector_[*p] / scalevector_[*p] / scalevalue_;
585 testvector_[*p] = resvector_[*p] / scalevalue_;
592 typename std::vector<int>::iterator p = curLSIdx_.begin();
593 for (; p < curLSIdx_.end(); ++p)
597 testvector_[*p] = resvector_[*p] / scalevalue_;
603 ind_.resize(curLSIdx_.size());
604 typename std::vector<int>::iterator p = curLSIdx_.begin();
605 for (; p < curLSIdx_.end(); ++p)
611 bool all_converged =
true;
612 for (
int ii = 0; ii < ensemble_size; ++ii)
614 if (ET::coeff(testvector_[*p], ii) > ET::coeff(tolerance_, ii))
616 ++ensemble_iterations[ii];
617 all_converged =
false;
619 else if (!(ET::coeff(testvector_[*p], ii) <= ET::coeff(tolerance_, ii)))
623 TEUCHOS_TEST_FOR_EXCEPTION(
true, StatusTestError,
"StatusTestWeightedGenResNorm::checkStatus(): NaN has been detected.");
634 int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
635 status_ = (have >= need) ? Passed : Failed;
641 template <
class ScalarType,
class MV,
class OP>
644 for (
int j = 0; j < indent; j++)
646 printStatus(os, status_);
648 if (status_ == Undefined)
649 os <<
", tol = " << tolerance_ << std::endl;
653 if (showMaxResNormOnly_ && curBlksz_ > 1)
656 testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
657 for (
int j = 0; j < indent + 13; j++)
659 os <<
"max{residual[" << curLSIdx_[0] <<
"..." << curLSIdx_[curBlksz_ - 1] <<
"]} = " << maxRelRes
660 << (maxRelRes <= tolerance_ ?
" <= " :
" > ") << tolerance_ << std::endl;
664 for (
int i = 0; i < numrhs_; i++)
666 for (
int j = 0; j < indent + 13; j++)
668 os <<
"residual [ " << i <<
" ] = " << testvector_[i];
669 os << ((testvector_[i] < tolerance_) ?
" < " : (testvector_[i] == tolerance_) ?
" == " : (testvector_[i] > tolerance_) ?
" > " :
" ") << tolerance_ << std::endl;
676 template <
class ScalarType,
class MV,
class OP>
679 os << std::left << std::setw(13) << std::setfill(
'.');
693 os << std::left << std::setfill(
' ');
697 template <
class ScalarType,
class MV,
class OP>
701 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
702 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
703 const LinearProblem<ScalarType, MV, OP> &lp = iSolver->getProblem();
705 if (firstcallCheckStatus_)
710 firstcallCheckStatus_ =
false;
712 if (scaletype_ == NormOfRHS)
714 Teuchos::RCP<const MV> rhs = lp.getRHS();
715 numrhs_ = MVT::GetNumberVecs(*rhs);
716 scalevector_.resize(numrhs_);
717 MvWNorm(*rhs, *weights_, scalevector_, scalenormtype_);
719 else if (scaletype_ == NormOfInitRes || scaletype_ == NormOfFullInitRes || scaletype_ == NormOfFullScaledInitRes)
721 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
722 numrhs_ = MVT::GetNumberVecs(*init_res);
723 scalevector_.resize(numrhs_);
724 MvWNorm(*init_res, *weights_, scalevector_, scalenormtype_);
726 else if (scaletype_ == NormOfPrecInitRes || scaletype_ == NormOfFullPrecInitRes || scaletype_ == NormOfFullScaledPrecInitRes)
728 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
729 numrhs_ = MVT::GetNumberVecs(*init_res);
730 scalevector_.resize(numrhs_);
731 MvWNorm(*init_res, *weights_, scalevector_, scalenormtype_);
735 numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
738 resvector_.resize(numrhs_);
739 testvector_.resize(numrhs_);
741 curLSNum_ = lp.getLSNumber();
742 curLSIdx_ = lp.getLSIndex();
743 curBlksz_ = (int)curLSIdx_.size();
745 for (i = 0; i < curBlksz_; ++i)
747 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
750 curNumRHS_ = validLS;
753 for (i = 0; i < numrhs_; i++)
755 testvector_[i] = one;
759 if (scalevalue_ == zero)