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> 60 #include <BelosStatusTestGenResSubNorm.hpp> 68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 class StatusTestGenResSubNorm<Scalar,
Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
Belos::OperatorT<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > >
70 :
public StatusTestResNorm<Scalar,Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>,Belos::OperatorT<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > > {
81 typedef MultiVecTraits<Scalar,MV>
MVT;
82 typedef OperatorTraits<Scalar,MV,OP>
OT;
99 StatusTestGenResSubNorm(
MagnitudeType Tolerance,
size_t subIdx,
int quorum = -1,
bool showMaxResNormOnly =
false )
100 : tolerance_(Tolerance),
103 showMaxResNormOnly_(showMaxResNormOnly),
104 resnormtype_(TwoNorm),
105 scaletype_(NormOfInitRes),
106 scalenormtype_(TwoNorm),
113 firstcallCheckStatus_(true),
114 firstcallDefineResForm_(true),
115 firstcallDefineScaleForm_(true),
116 mapExtractor_(
Teuchos::null) { }
119 virtual ~StatusTestGenResSubNorm() { };
132 int defineResForm(NormType TypeOfNorm) {
134 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
135 firstcallDefineResForm_ =
false;
137 resnormtype_ = TypeOfNorm;
165 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
166 firstcallDefineScaleForm_ =
false;
168 scaletype_ = TypeOfScaling;
169 scalenormtype_ = TypeOfNorm;
170 scalevalue_ = ScaleValue;
179 int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance;
return(0);}
184 int setSubIdx (
size_t subIdx ) { subIdx_ = subIdx;
return(0);}
188 int setQuorum(
int quorum) {quorum_ = quorum;
return(0);}
191 int setShowMaxResNormOnly(
bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly;
return(0);}
204 StatusType checkStatus(Iteration<Scalar,MV,OP>* iSolver) {
206 const LinearProblem<Scalar,MV,OP>& lp = iSolver->getProblem();
208 if (firstcallCheckStatus_) {
209 StatusType status = firstCallCheckStatusSetup(iSolver);
219 if ( curLSNum_ != lp.getLSNumber() ) {
223 curLSNum_ = lp.getLSNumber();
224 curLSIdx_ = lp.getLSIndex();
225 curBlksz_ = (int)curLSIdx_.size();
227 for (
int i=0; i<curBlksz_; ++i) {
228 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
231 curNumRHS_ = validLS;
232 curSoln_ = Teuchos::null;
238 if (status_==Passed) {
return status_; }
245 curSoln_ = lp.updateSolution( cur_update );
246 Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
247 lp.computeCurrResVec( &*cur_res, &*curSoln_ );
248 std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
249 MvSubNorm( *cur_res, subIdx_, tmp_resvector, resnormtype_ );
251 typename std::vector<int>::iterator p = curLSIdx_.begin();
252 for (
int i=0; p<curLSIdx_.end(); ++p, ++i) {
255 resvector_[*p] = tmp_resvector[i];
262 if ( scalevector_.size() > 0 ) {
263 typename std::vector<int>::iterator pp = curLSIdx_.begin();
264 for (; pp<curLSIdx_.end(); ++pp) {
268 if ( scalevector_[ *pp ] != zero ) {
270 testvector_[ *pp ] = resvector_[ *pp ] / scalevector_[ *pp ] / scalevalue_;
272 testvector_[ *pp ] = resvector_[ *pp ] / scalevalue_;
278 typename std::vector<int>::iterator pp = curLSIdx_.begin();
279 for (; pp<curLSIdx_.end(); ++pp) {
282 testvector_[ *pp ] = resvector_[ *pp ] / scalevalue_;
287 ind_.resize( curLSIdx_.size() );
288 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
289 for (; p2<curLSIdx_.end(); ++p2) {
293 if (testvector_[ *p2 ] > tolerance_) {
297 }
else if (testvector_[ *p2 ] <= tolerance_) {
308 int need = (quorum_ == -1) ? curNumRHS_: quorum_;
309 status_ = (have >= need) ? Passed : Failed;
315 StatusType getStatus()
const {
return(status_);};
329 firstcallCheckStatus_ =
true;
330 curSoln_ = Teuchos::null;
339 void print(std::ostream& os,
int indent = 0)
const {
340 os.setf(std::ios_base::scientific);
341 for (
int j = 0; j < indent; j ++)
343 printStatus(os, status_);
346 os <<
", tol = " << tolerance_ << std::endl;
349 if(showMaxResNormOnly_ && curBlksz_ > 1) {
350 const MagnitudeType maxRelRes = *std::max_element(
351 testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
353 for (
int j = 0; j < indent + 13; j ++)
355 os <<
"max{residual["<<curLSIdx_[0]<<
"..."<<curLSIdx_[curBlksz_-1]<<
"]} = " << maxRelRes
356 << ( maxRelRes <= tolerance_ ?
" <= " :
" > " ) << tolerance_ << std::endl;
359 for (
int i=0; i<numrhs_; i++ ) {
360 for (
int j = 0; j < indent + 13; j ++)
362 os <<
"residual [ " << i <<
" ] = " << testvector_[ i ];
363 os << ((testvector_[i]<tolerance_) ?
" < " : (testvector_[i]==tolerance_) ?
" == " : (testvector_[i]>tolerance_) ?
" > " :
" " ) << tolerance_ << std::endl;
371 void printStatus(std::ostream& os, StatusType type)
const {
372 os << std::left << std::setw(13) << std::setfill(
'.');
385 os << std::left << std::setfill(
' ');
398 int getQuorum()
const {
return quorum_; }
401 size_t getSubIdx()
const {
return subIdx_; }
404 bool getShowMaxResNormOnly() {
return showMaxResNormOnly_; }
407 std::vector<int> convIndices() {
return ind_; }
410 MagnitudeType getTolerance()
const {
return(tolerance_);};
413 const std::vector<MagnitudeType>* getTestValue()
const {
return(&testvector_);};
416 const std::vector<MagnitudeType>* getResNormValue()
const {
return(&resvector_);};
419 const std::vector<MagnitudeType>* getScaledNormValue()
const {
return(&scalevector_);};
423 bool getLOADetected()
const {
return false; }
436 StatusType firstCallCheckStatusSetup(Iteration<Scalar,MV,OP>* iSolver) {
440 const LinearProblem<Scalar,MV,OP>& lp = iSolver->getProblem();
442 if (firstcallCheckStatus_) {
446 firstcallCheckStatus_ =
false;
458 mapExtractor_ = bMat->getRangeMapExtractor();
463 if (scaletype_== NormOfRHS) {
465 numrhs_ = MVT::GetNumberVecs( *rhs );
466 scalevector_.resize( numrhs_ );
467 MvSubNorm( *rhs, subIdx_, scalevector_, scalenormtype_ );
469 else if (scaletype_==NormOfInitRes) {
471 numrhs_ = MVT::GetNumberVecs( *init_res );
472 scalevector_.resize( numrhs_ );
473 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
475 else if (scaletype_==NormOfPrecInitRes) {
477 numrhs_ = MVT::GetNumberVecs( *init_res );
478 scalevector_.resize( numrhs_ );
479 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
481 else if (scaletype_==NormOfFullInitRes) {
483 numrhs_ = MVT::GetNumberVecs( *init_res );
484 scalevector_.resize( numrhs_ );
485 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
488 else if (scaletype_==NormOfFullPrecInitRes) {
490 numrhs_ = MVT::GetNumberVecs( *init_res );
491 scalevector_.resize( numrhs_ );
492 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
495 else if (scaletype_==NormOfFullScaledInitRes) {
497 numrhs_ = MVT::GetNumberVecs( *init_res );
498 scalevector_.resize( numrhs_ );
499 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
500 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
502 else if (scaletype_==NormOfFullScaledPrecInitRes) {
504 numrhs_ = MVT::GetNumberVecs( *init_res );
505 scalevector_.resize( numrhs_ );
506 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
507 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
510 numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
513 resvector_.resize( numrhs_ );
514 testvector_.resize( numrhs_ );
516 curLSNum_ = lp.getLSNumber();
517 curLSIdx_ = lp.getLSIndex();
518 curBlksz_ = (int)curLSIdx_.size();
520 for (i=0; i<curBlksz_; ++i) {
521 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
524 curNumRHS_ = validLS;
527 for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
530 if (scalevalue_ == zero) {
542 std::string description()
const 544 std::ostringstream oss;
545 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
546 oss <<
", tol = " << tolerance_;
558 std::string resFormStr()
const 560 std::ostringstream oss;
562 oss << ((resnormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
564 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
567 if (scaletype_!=
None)
573 if (scaletype_==UserProvided)
574 oss <<
" (User Scale)";
577 oss << ((scalenormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
578 if (scaletype_==NormOfInitRes)
579 oss <<
" Res0 [" << subIdx_ <<
"]";
580 else if (scaletype_==NormOfPrecInitRes)
581 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
582 else if (scaletype_==NormOfFullInitRes)
583 oss <<
" Full Res0 [" << subIdx_ <<
"]";
584 else if (scaletype_==NormOfFullPrecInitRes)
585 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
586 else if (scaletype_==NormOfFullScaledInitRes)
587 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
588 else if (scaletype_==NormOfFullScaledPrecInitRes)
589 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
591 oss <<
" RHS [" << subIdx_ <<
"]";
612 typedef MultiVecTraits<Scalar, MV> MVT;
613 MVT::MvNorm(*SubVec,normVec,type);
617 void MvScalingRatio(
const MV& mv,
size_t block, MagnitudeType& lengthRatio) {
622 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
630 MagnitudeType tolerance_;
639 bool showMaxResNormOnly_;
642 NormType resnormtype_;
645 ScaleType scaletype_;
648 NormType scalenormtype_;
651 MagnitudeType scalevalue_;
654 std::vector<MagnitudeType> scalevector_;
657 std::vector<MagnitudeType> resvector_;
660 std::vector<MagnitudeType> testvector_;
663 std::vector<int> ind_;
678 std::vector<int> curLSIdx_;
687 bool firstcallCheckStatus_;
690 bool firstcallDefineResForm_;
693 bool firstcallDefineScaleForm_;
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > ME
Exception indicating invalid cast attempted.
SCT::magnitudeType MagnitudeType
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MV
OperatorTraits< Scalar, MV, OP > OT
Implementation of the Belos::XpetraOp. It derives from the Belos::OperatorT templated on the Xpetra::...
Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > BCRS
Teuchos::ScalarTraits< Scalar > SCT
MultiVecTraits< Scalar, MV > MVT
Exception throws to report errors in the internal logical of the program.
Belos::OperatorT< MV > OP