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>
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> > > {
77 typedef Belos::OperatorT<MV>
OP;
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;
451 Teuchos::rcp_dynamic_cast<
const Belos::XpetraOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
461 mapExtractor_ = bMat->getRangeMapExtractor();
466 if (scaletype_== NormOfRHS) {
468 numrhs_ = MVT::GetNumberVecs( *rhs );
469 scalevector_.resize( numrhs_ );
470 MvSubNorm( *rhs, subIdx_, scalevector_, scalenormtype_ );
472 else if (scaletype_==NormOfInitRes) {
474 numrhs_ = MVT::GetNumberVecs( *init_res );
475 scalevector_.resize( numrhs_ );
476 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
478 else if (scaletype_==NormOfPrecInitRes) {
480 numrhs_ = MVT::GetNumberVecs( *init_res );
481 scalevector_.resize( numrhs_ );
482 MvSubNorm( *init_res, subIdx_, scalevector_, scalenormtype_ );
484 else if (scaletype_==NormOfFullInitRes) {
486 numrhs_ = MVT::GetNumberVecs( *init_res );
487 scalevector_.resize( numrhs_ );
488 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
491 else if (scaletype_==NormOfFullPrecInitRes) {
493 numrhs_ = MVT::GetNumberVecs( *init_res );
494 scalevector_.resize( numrhs_ );
495 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
498 else if (scaletype_==NormOfFullScaledInitRes) {
500 numrhs_ = MVT::GetNumberVecs( *init_res );
501 scalevector_.resize( numrhs_ );
502 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
503 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
505 else if (scaletype_==NormOfFullScaledPrecInitRes) {
507 numrhs_ = MVT::GetNumberVecs( *init_res );
508 scalevector_.resize( numrhs_ );
509 MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
510 MvScalingRatio( *init_res, subIdx_, scalevalue_ );
513 numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
516 resvector_.resize( numrhs_ );
517 testvector_.resize( numrhs_ );
519 curLSNum_ = lp.getLSNumber();
520 curLSIdx_ = lp.getLSIndex();
521 curBlksz_ = (int)curLSIdx_.size();
523 for (i=0; i<curBlksz_; ++i) {
524 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
527 curNumRHS_ = validLS;
530 for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
533 if (scalevalue_ == zero) {
545 std::string description()
const
547 std::ostringstream oss;
548 oss <<
"Belos::StatusTestGenResSubNorm<>: " << resFormStr();
549 oss <<
", tol = " << tolerance_;
561 std::string resFormStr()
const
563 std::ostringstream oss;
565 oss << ((resnormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
567 oss <<
" Res Vec [" << subIdx_ <<
"]) ";
570 if (scaletype_!=
None)
576 if (scaletype_==UserProvided)
577 oss <<
" (User Scale)";
580 oss << ((scalenormtype_==OneNorm) ?
"1-Norm" : (resnormtype_==TwoNorm) ?
"2-Norm" :
"Inf-Norm");
581 if (scaletype_==NormOfInitRes)
582 oss <<
" Res0 [" << subIdx_ <<
"]";
583 else if (scaletype_==NormOfPrecInitRes)
584 oss <<
" Prec Res0 [" << subIdx_ <<
"]";
585 else if (scaletype_==NormOfFullInitRes)
586 oss <<
" Full Res0 [" << subIdx_ <<
"]";
587 else if (scaletype_==NormOfFullPrecInitRes)
588 oss <<
" Full Prec Res0 [" << subIdx_ <<
"]";
589 else if (scaletype_==NormOfFullScaledInitRes)
590 oss <<
" scaled Full Res0 [" << subIdx_ <<
"]";
591 else if (scaletype_==NormOfFullScaledPrecInitRes)
592 oss <<
" scaled Full Prec Res0 [" << subIdx_ <<
"]";
594 oss <<
" RHS [" << subIdx_ <<
"]";
615 typedef MultiVecTraits<Scalar, MV> MVT;
616 MVT::MvNorm(*SubVec,normVec,type);
620 void MvScalingRatio(
const MV& mv,
size_t block, MagnitudeType& lengthRatio) {
625 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
633 MagnitudeType tolerance_;
642 bool showMaxResNormOnly_;
645 NormType resnormtype_;
648 ScaleType scaletype_;
651 NormType scalenormtype_;
654 MagnitudeType scalevalue_;
657 std::vector<MagnitudeType> scalevector_;
660 std::vector<MagnitudeType> resvector_;
663 std::vector<MagnitudeType> testvector_;
666 std::vector<int> ind_;
681 std::vector<int> curLSIdx_;
690 bool firstcallCheckStatus_;
693 bool firstcallDefineResForm_;
696 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
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