Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosStatusTestGenResNorm.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_STATUS_TEST_GEN_RESNORM_H
11 #define BELOS_STATUS_TEST_GEN_RESNORM_H
12 
19 #include "BelosLinearProblem.hpp"
20 #include "BelosMultiVecTraits.hpp"
21 
44 namespace Belos {
45 
46 template <class ScalarType, class MV, class OP>
47 class StatusTestGenResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
48 
49  public:
50 
51  // Convenience typedefs
55 
57 
58 
62  enum ResType {Implicit,
65  };
66 
68 
70 
71 
85  StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
86 
88  virtual ~StatusTestGenResNorm();
90 
92 
93 
95 
102  int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
103 
105 
125  int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
126 
127  NormType getResNormType() {return resnormtype_;}
128 
130 
133  int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
134 
137  int setQuorum(int quorum) {quorum_ = quorum; return(0);}
138 
140  int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
141 
143 
145 
146 
154 
156  StatusType getStatus() const {return(status_);};
158 
160 
161 
163  void reset();
164 
166 
168 
169 
171  void print(std::ostream& os, int indent = 0) const;
172 
174  void printStatus(std::ostream& os, StatusType type) const;
176 
178 
179 
183  Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
184 
187  int getQuorum() const { return quorum_; }
188 
190  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
191 
193  std::vector<int> convIndices() { return ind_; }
194 
196  MagnitudeType getTolerance() const {return(tolerance_);};
197 
199  const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
200 
202  const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
203 
205  const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
206 
209  bool getLOADetected() const { return false; }
210 
212 
213 
216 
224 
227 
229  std::string description() const
230  {
231  std::ostringstream oss;
232  oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
233  oss << ", tol = " << tolerance_;
234  return oss.str();
235  }
237 
238  protected:
239 
240  private:
241 
243 
244 
245  std::string resFormStr() const
246  {
247  std::ostringstream oss;
248  oss << "(";
249  oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
250  oss << ((restype_==Explicit) ? " Exp" : " Imp");
251  oss << " Res Vec) ";
252 
253  // If there is no residual scaling, return current string.
254  if (scaletype_!=None)
255  {
256  // Insert division sign.
257  oss << "/ ";
258 
259  // Determine output string for scaling, if there is any.
260  if (scaletype_==UserProvided)
261  oss << " (User Scale)";
262  else {
263  oss << "(";
264  oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
265  if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes)
266  oss << " Res0";
267  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes)
268  oss << " Prec Res0";
269  else
270  oss << " RHS ";
271  oss << ")";
272  }
273  }
274 
275  return oss.str();
276  }
277 
279 
281 
282 
284  MagnitudeType tolerance_;
285 
287  int quorum_;
288 
290  bool showMaxResNormOnly_;
291 
293  ResType restype_;
294 
296  NormType resnormtype_;
297 
299  ScaleType scaletype_;
300 
302  NormType scalenormtype_;
303 
305  MagnitudeType scalevalue_;
306 
308  std::vector<MagnitudeType> scalevector_;
309 
311  std::vector<MagnitudeType> resvector_;
312 
314  std::vector<MagnitudeType> testvector_;
315 
317  std::vector<int> ind_;
318 
320  Teuchos::RCP<MV> curSoln_;
321 
323  StatusType status_;
324 
326  int curBlksz_;
327 
329  int curNumRHS_;
330 
332  std::vector<int> curLSIdx_;
333 
335  int curLSNum_;
336 
338  int numrhs_;
339 
341  bool firstcallCheckStatus_;
342 
344  bool firstcallDefineResForm_;
345 
347  bool firstcallDefineScaleForm_;
348 
350 
351 };
352 
353 template <class ScalarType, class MV, class OP>
355 StatusTestGenResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
356  : tolerance_(Tolerance),
357  quorum_(quorum),
358  showMaxResNormOnly_(showMaxResNormOnly),
359  restype_(Implicit),
360  resnormtype_(TwoNorm),
361  scaletype_(NormOfInitRes),
362  scalenormtype_(TwoNorm),
363  scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
364  status_(Undefined),
365  curBlksz_(0),
366  curNumRHS_(0),
367  curLSNum_(0),
368  numrhs_(0),
369  firstcallCheckStatus_(true),
370  firstcallDefineResForm_(true),
371  firstcallDefineScaleForm_(true)
372 {
373  // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
374  // the implicit residual std::vector.
375 }
376 
377 template <class ScalarType, class MV, class OP>
379 {}
380 
381 template <class ScalarType, class MV, class OP>
383 {
384  status_ = Undefined;
385  curBlksz_ = 0;
386  curLSNum_ = 0;
387  curLSIdx_.resize(0);
388  numrhs_ = 0;
389  ind_.resize(0);
390  firstcallCheckStatus_ = true;
391  curSoln_ = Teuchos::null;
392 }
393 
394 template <class ScalarType, class MV, class OP>
396 {
397  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
398  "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
399  firstcallDefineResForm_ = false;
400 
401  restype_ = TypeOfResidual;
402  resnormtype_ = TypeOfNorm;
403 
404  return(0);
405 }
406 
407 template <class ScalarType, class MV, class OP>
409  MagnitudeType ScaleValue )
410 {
411  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
412  "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
413  firstcallDefineScaleForm_ = false;
414 
415  scaletype_ = TypeOfScaling;
416  scalenormtype_ = TypeOfNorm;
417  scalevalue_ = ScaleValue;
418 
419  return(0);
420 }
421 
422 template <class ScalarType, class MV, class OP>
424 {
426  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
427  // Compute scaling term (done once for each block that's being solved)
428  if (firstcallCheckStatus_) {
429  StatusType status = firstCallCheckStatusSetup(iSolver);
430  if(status==Failed) {
431  status_ = Failed;
432  return(status_);
433  }
434  }
435  //
436  // This section computes the norm of the residual std::vector
437  //
438  if ( curLSNum_ != lp.getLSNumber() ) {
439  //
440  // We have moved on to the next rhs block
441  //
442  curLSNum_ = lp.getLSNumber();
443  curLSIdx_ = lp.getLSIndex();
444  curBlksz_ = (int)curLSIdx_.size();
445  int validLS = 0;
446  for (int i=0; i<curBlksz_; ++i) {
447  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
448  validLS++;
449  }
450  curNumRHS_ = validLS;
451  curSoln_ = Teuchos::null;
452  //
453  } else {
454  //
455  // We are in the same rhs block, return if we are converged
456  //
457  if (status_==Passed) { return status_; }
458  }
459  if (restype_==Implicit) {
460  //
461  // get the native residual norms from the solver for this block of right-hand sides.
462  // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
463  // Otherwise the native residual is assumed to be stored in the resvector_.
464  //
465  std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
466  Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );
467  if ( residMV != Teuchos::null ) {
468  tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
469  MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );
470  typename std::vector<int>::iterator p = curLSIdx_.begin();
471  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
472  // Check if this index is valid
473  if (*p != -1)
474  resvector_[*p] = tmp_resvector[i];
475  }
476  } else {
477  typename std::vector<int>::iterator p = curLSIdx_.begin();
478  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
479  // Check if this index is valid
480  if (*p != -1)
481  resvector_[*p] = tmp_resvector[i];
482  }
483  }
484  }
485  else if (restype_==Explicit) {
486  //
487  // Request the true residual for this block of right-hand sides.
488  //
489  Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
490  curSoln_ = lp.updateSolution( cur_update );
491  Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
492  lp.computeCurrResVec( &*cur_res, &*curSoln_ );
493  std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
494  MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
495  typename std::vector<int>::iterator p = curLSIdx_.begin();
496  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
497  // Check if this index is valid
498  if (*p != -1)
499  resvector_[*p] = tmp_resvector[i];
500  }
501  }
502  //
503  // Compute the new linear system residuals for testing.
504  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
505  //
506  if ( scalevector_.size() > 0 ) {
507  typename std::vector<int>::iterator p = curLSIdx_.begin();
508  for (; p<curLSIdx_.end(); ++p) {
509  // Check if this index is valid
510  if (*p != -1) {
511  // Scale the std::vector accordingly
512  testvector_[ *p ] =
513  scalevector_[ *p ] != zero? resvector_[ *p ] / (scalevector_[ *p ] * scalevalue_) : resvector_[ *p ] / scalevalue_;
514  }
515  }
516  }
517  else {
518  typename std::vector<int>::iterator p = curLSIdx_.begin();
519  for (; p<curLSIdx_.end(); ++p) {
520  // Check if this index is valid
521  if (*p != -1)
522  testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
523  }
524  }
525 
526  // Check status of new linear system residuals and see if we have the quorum.
527  int have = 0;
528  ind_.resize( curLSIdx_.size() );
529  typename std::vector<int>::iterator p = curLSIdx_.begin();
530  for (; p<curLSIdx_.end(); ++p) {
531  // Check if this index is valid
532  if (*p != -1) {
533  // Check if any of the residuals are larger than the tolerance.
534  if (testvector_[ *p ] > tolerance_) {
535  // do nothing.
536  } else if (testvector_[ *p ] <= tolerance_) {
537  ind_[have] = *p;
538  have++;
539  } else {
540  // Throw an std::exception if a NaN is found.
541  status_ = Failed;
542  TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestNaNError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
543  }
544  }
545  }
546  ind_.resize(have);
547  int need = (quorum_ == -1) ? curNumRHS_: quorum_;
548  status_ = (have >= need) ? Passed : Failed;
549 
550  // Return the current status
551  return status_;
552 }
553 
554 template <class ScalarType, class MV, class OP>
555 void StatusTestGenResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
556 {
557  for (int j = 0; j < indent; j ++)
558  os << ' ';
559  printStatus(os, status_);
560  os << resFormStr();
561  if (status_==Undefined)
562  os << ", tol = " << tolerance_ << std::endl;
563  else {
564  os << std::endl;
565  if(showMaxResNormOnly_ && curBlksz_ > 1) {
566  const MagnitudeType maxRelRes = *std::max_element(
567  testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
568  );
569  for (int j = 0; j < indent + 13; j ++)
570  os << ' ';
571  os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
572  << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
573  }
574  else {
575  for ( int i=0; i<numrhs_; i++ ) {
576  for (int j = 0; j < indent + 13; j ++)
577  os << ' ';
578  os << "residual [ " << i << " ] = " << testvector_[ i ];
579  os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " " ) << tolerance_ << std::endl;
580  }
581  }
582  }
583  os << std::endl;
584 }
585 
586 template <class ScalarType, class MV, class OP>
588 {
589  os << std::left << std::setw(13) << std::setfill('.');
590  switch (type) {
591  case Passed:
592  os << "Converged";
593  break;
594  case Failed:
595  os << "Unconverged";
596  break;
597  case Undefined:
598  default:
599  os << "**";
600  break;
601  }
602  os << std::left << std::setfill(' ');
603  return;
604 }
605 
606 template <class ScalarType, class MV, class OP>
608 {
609  int i;
612  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
613  // Compute scaling term (done once for each block that's being solved)
614  if (firstcallCheckStatus_) {
615  //
616  // Get some current solver information.
617  //
618  firstcallCheckStatus_ = false;
619 
620  if (scaletype_== NormOfRHS) {
621  Teuchos::RCP<const MV> rhs = lp.getRHS();
622  numrhs_ = MVT::GetNumberVecs( *rhs );
623  scalevector_.resize( numrhs_ );
624  MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
625  }
626  else if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes) {
627  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
628  numrhs_ = MVT::GetNumberVecs( *init_res );
629  scalevector_.resize( numrhs_ );
630  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
631  }
632  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes) {
634  numrhs_ = MVT::GetNumberVecs( *init_res );
635  scalevector_.resize( numrhs_ );
636  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
637  }
638  else {
639  numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
640  }
641 
642  resvector_.resize( numrhs_ );
643  testvector_.resize( numrhs_ );
644 
645  curLSNum_ = lp.getLSNumber();
646  curLSIdx_ = lp.getLSIndex();
647  curBlksz_ = (int)curLSIdx_.size();
648  int validLS = 0;
649  for (i=0; i<curBlksz_; ++i) {
650  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
651  validLS++;
652  }
653  curNumRHS_ = validLS;
654  //
655  // Initialize the testvector.
656  for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
657 
658  // Return an error if the scaling is zero.
659  if (scalevalue_ == zero) {
660  return Failed;
661  }
662  }
663  return Undefined;
664 }
665 
666 } // end namespace Belos
667 
668 #endif /* BELOS_STATUS_TEST_RESNORM_H */
virtual Teuchos::RCP< const MV > getNativeResiduals(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *norms) const =0
Get the residuals native to the solver.
int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue=Teuchos::ScalarTraits< MagnitudeType >::one())
Define form of the scaling, its norm, its optional weighting std::vector, or, alternatively, define an explicit value.
MagnitudeType getTolerance() const
Returns the value of the tolerance, , set in the constructor.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:88
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
int setShowMaxResNormOnly(bool showMaxResNormOnly)
Set whether the only maximum residual norm is displayed when the print() method is called...
StatusType firstCallCheckStatusSetup(Iteration< ScalarType, MV, OP > *iSolver)
Call to setup initial scaling std::vector.
virtual Teuchos::RCP< MV > updateSolution(const Teuchos::RCP< MV > &update=Teuchos::null, bool updateLP=false, ScalarType scale=Teuchos::ScalarTraits< ScalarType >::one())
Compute the new solution to the linear system using the given update vector.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual Teuchos::RCP< MV > getCurrentUpdate() const =0
Get the current update to the linear system.
Declaration of basic traits for the multivector type.
std::string description() const
Method to return description of the maximum iteration status test.
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
An abstract class of StatusTest for stopping criteria using residual norms.
void printStatus(std::ostream &os, StatusType type) const
Print message for each status specific to this stopping test.
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
StatusType
Whether the StatusTest wants iteration to stop.
Definition: BelosTypes.hpp:157
const std::vector< MagnitudeType > * getResNormValue() const
Returns the residual norm value, , computed in most recent call to CheckStatus.
Traits class which defines basic operations on multivectors.
const std::vector< int > & getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
MultiVecTraits< ScalarType, MV > MVT
bool getShowMaxResNormOnly()
Returns whether the only maximum residual norm is displayed when the print() method is called...
int defineResForm(ResType TypeOfResidual, NormType TypeOfNorm)
Define form of the residual, its norm and optional weighting std::vector.
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
int getLSNumber() const
The number of linear systems that have been set.
int setQuorum(int quorum)
Sets the number of residuals that must pass the convergence test before Passed is returned...
std::vector< int > convIndices()
Returns the std::vector containing the indices of the residuals that passed the test.
virtual const LinearProblem< ScalarType, MV, OP > & getProblem() const =0
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
Belos::StatusTest abstract class for specifying a residual norm stopping criteria.
StatusType getStatus() const
Return the result of the most recent CheckStatus call.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:65
StatusType checkStatus(Iteration< ScalarType, MV, OP > *iSolver)
Check convergence status: Passed, Failed, or Undefined.
virtual void computeCurrResVec(MV *R, const MV *X=0, const MV *B=0) const
Compute a residual R for this operator given a solution X, and right-hand side B. ...
void print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
ResType
Select how the residual std::vector is produced.
void reset()
Resets the internal configuration to the initial state.
int getQuorum() const
Returns the number of residuals that must pass the convergence test before Passed is returned...
StatusTestGenResNorm(MagnitudeType Tolerance, int quorum=-1, bool showMaxResNormOnly=false)
Constructor.
const std::vector< MagnitudeType > * getScaledNormValue() const
Returns the scaled norm value, .
Teuchos::RCP< MV > getSolution()
Returns the current solution estimate that was computed for the most recent residual test...
bool getLOADetected() const
Returns a boolean indicating a loss of accuracy has been detected in computing the residual...

Generated on Fri Nov 22 2024 09:23:06 for Belos by doxygen 1.8.5