Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Belos_StatusTest_GenResNorm_MP_Vector.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
43 #define BELOS_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
44 
54 #include "BelosStatusTestGenResNorm.hpp"
55 
82 namespace Belos {
83 
84 template <class Storage, class MV, class OP>
85 class StatusTestGenResNorm<Sacado::MP::Vector<Storage>, MV, OP> :
86  public StatusTestResNorm<Sacado::MP::Vector<Storage>,MV,OP> {
87 
88  public:
89 
90  // Convenience typedefs
92  typedef Teuchos::ScalarTraits<ScalarType> SCT;
93  typedef typename SCT::magnitudeType MagnitudeType;
94  typedef MultiVecTraits<ScalarType,MV> MVT;
95 
97 
98 
102  enum ResType {Implicit,
103  Explicit
105  };
106 
108 
110 
111 
125  StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
126 
128  virtual ~StatusTestGenResNorm();
130 
132 
133 
135 
142  int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
143 
145 
165  int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
166 
168 
171  int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
172 
175  int setQuorum(int quorum) {quorum_ = quorum; return(0);}
176 
178  int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
179 
181 
183 
184 
191  StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
192 
194  StatusType getStatus() const {return(status_);}
196 
198 
199 
201  void reset();
202 
204 
206 
207 
209  void print(std::ostream& os, int indent = 0) const;
210 
212  void printStatus(std::ostream& os, StatusType type) const;
214 
216 
217 
221  Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
222 
225  int getQuorum() const { return quorum_; }
226 
228  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
229 
231  std::vector<int> convIndices() { return ind_; }
232 
234  MagnitudeType getTolerance() const {return(tolerance_);}
235 
237  const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);}
238 
240  const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);}
241 
243  const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);}
244 
247  bool getLOADetected() const { return false; }
248 
250  const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
251 
253 
254 
257 
263  StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
265 
268 
270  std::string description() const
271  {
272  std::ostringstream oss;
273  oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
274  oss << ", tol = " << tolerance_;
275  return oss.str();
276  }
278 
279  protected:
280 
281  private:
282 
284 
285 
286  std::string resFormStr() const
287  {
288  std::ostringstream oss;
289  oss << "(";
290  oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
291  oss << ((restype_==Explicit) ? " Exp" : " Imp");
292  oss << " Res Vec) ";
293 
294  // If there is no residual scaling, return current string.
295  if (scaletype_!=None)
296  {
297  // Insert division sign.
298  oss << "/ ";
299 
300  // Determine output string for scaling, if there is any.
301  if (scaletype_==UserProvided)
302  oss << " (User Scale)";
303  else {
304  oss << "(";
305  oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
306  if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes)
307  oss << " Res0";
308  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes)
309  oss << " Prec Res0";
310  else
311  oss << " RHS ";
312  oss << ")";
313  }
314  }
315 
316  return oss.str();
317  }
318 
320 
322 
323 
326 
328  int quorum_;
329 
332 
335 
337  NormType resnormtype_;
338 
340  ScaleType scaletype_;
341 
343  NormType scalenormtype_;
344 
347 
349  std::vector<MagnitudeType> scalevector_;
350 
352  std::vector<MagnitudeType> resvector_;
353 
355  std::vector<MagnitudeType> testvector_;
356 
358  std::vector<int> ind_;
359 
361  Teuchos::RCP<MV> curSoln_;
362 
364  StatusType status_;
365 
368 
371 
373  std::vector<int> curLSIdx_;
374 
377 
379  int numrhs_;
380 
383 
386 
389 
391  std::vector<int> ensemble_converged;
392 
394  std::vector<int> ensemble_iterations;
395 
397 
398 };
399 
400 template <class StorageType, class MV, class OP>
402 StatusTestGenResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
403  : tolerance_(Tolerance),
404  quorum_(quorum),
405  showMaxResNormOnly_(showMaxResNormOnly),
406  restype_(Implicit),
407  resnormtype_(TwoNorm),
408  scaletype_(NormOfInitRes),
409  scalenormtype_(TwoNorm),
410  scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
411  status_(Undefined),
412  curBlksz_(0),
413  curNumRHS_(0),
414  curLSNum_(0),
415  numrhs_(0),
416  firstcallCheckStatus_(true),
417  firstcallDefineResForm_(true),
418  firstcallDefineScaleForm_(true),
419  ensemble_converged(StorageType::static_size, 0),
420  ensemble_iterations(StorageType::static_size, 0)
421 {
422  // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
423  // the implicit residual std::vector.
424 }
425 
426 template <class StorageType, class MV, class OP>
427 StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::~StatusTestGenResNorm()
428 {}
429 
430 template <class StorageType, class MV, class OP>
431 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::reset()
432 {
433  status_ = Undefined;
434  curBlksz_ = 0;
435  curLSNum_ = 0;
436  curLSIdx_.resize(0);
437  numrhs_ = 0;
438  ind_.resize(0);
439  firstcallCheckStatus_ = true;
440  curSoln_ = Teuchos::null;
441  ensemble_converged = std::vector<int>(StorageType::static_size, 0);
442  ensemble_iterations = std::vector<int>(StorageType::static_size, 0);
443 }
444 
445 template <class StorageType, class MV, class OP>
446 int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
447 {
448  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
449  "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
450  firstcallDefineResForm_ = false;
451 
452  restype_ = TypeOfResidual;
453  resnormtype_ = TypeOfNorm;
454 
455  return(0);
456 }
457 
458 template <class StorageType, class MV, class OP>
459 int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
460  MagnitudeType ScaleValue )
461 {
462  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
463  "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
464  firstcallDefineScaleForm_ = false;
465 
466  scaletype_ = TypeOfScaling;
467  scalenormtype_ = TypeOfNorm;
468  scalevalue_ = ScaleValue;
469 
470  return(0);
471 }
472 
473 template <class StorageType, class MV, class OP>
474 StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
475 {
476  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
477  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
478  // Compute scaling term (done once for each block that's being solved)
479  if (firstcallCheckStatus_) {
480  StatusType status = firstCallCheckStatusSetup(iSolver);
481  if(status==Failed) {
482  status_ = Failed;
483  return(status_);
484  }
485  }
486  //
487  // This section computes the norm of the residual std::vector
488  //
489  if ( curLSNum_ != lp.getLSNumber() ) {
490  //
491  // We have moved on to the next rhs block
492  //
493  curLSNum_ = lp.getLSNumber();
494  curLSIdx_ = lp.getLSIndex();
495  curBlksz_ = (int)curLSIdx_.size();
496  int validLS = 0;
497  for (int i=0; i<curBlksz_; ++i) {
498  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
499  validLS++;
500  }
501  curNumRHS_ = validLS;
502  curSoln_ = Teuchos::null;
503  //
504  } else {
505  //
506  // We are in the same rhs block, return if we are converged
507  //
508  if (status_==Passed) { return status_; }
509  }
510  if (restype_==Implicit) {
511  //
512  // get the native residual norms from the solver for this block of right-hand sides.
513  // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
514  // Otherwise the native residual is assumed to be stored in the resvector_.
515  //
516  std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
517  Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );
518  if ( residMV != Teuchos::null ) {
519  tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
520  MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );
521  typename std::vector<int>::iterator p = curLSIdx_.begin();
522  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
523  // Check if this index is valid
524  if (*p != -1)
525  resvector_[*p] = tmp_resvector[i];
526  }
527  } else {
528  typename std::vector<int>::iterator p = curLSIdx_.begin();
529  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
530  // Check if this index is valid
531  if (*p != -1)
532  resvector_[*p] = tmp_resvector[i];
533  }
534  }
535  }
536  else if (restype_==Explicit) {
537  //
538  // Request the true residual for this block of right-hand sides.
539  //
540  Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
541  curSoln_ = lp.updateSolution( cur_update );
542  Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
543  lp.computeCurrResVec( &*cur_res, &*curSoln_ );
544  std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
545  MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
546  typename std::vector<int>::iterator p = curLSIdx_.begin();
547  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
548  // Check if this index is valid
549  if (*p != -1)
550  resvector_[*p] = tmp_resvector[i];
551  }
552  }
553  //
554  // Compute the new linear system residuals for testing.
555  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
556  //
557  if ( scalevector_.size() > 0 ) {
558  typename std::vector<int>::iterator p = curLSIdx_.begin();
559  for (; p<curLSIdx_.end(); ++p) {
560  // Check if this index is valid
561  if (*p != -1) {
562  // Scale the std::vector accordingly
563  if ( scalevector_[ *p ] != zero ) {
564  // Don't intentionally divide by zero.
565  testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
566  } else {
567  testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
568  }
569  }
570  }
571  }
572  else {
573  typename std::vector<int>::iterator p = curLSIdx_.begin();
574  for (; p<curLSIdx_.end(); ++p) {
575  // Check if this index is valid
576  if (*p != -1)
577  testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
578  }
579  }
580 
581  // Check status of new linear system residuals and see if we have the quorum.
582  int have = 0;
583  ind_.resize( curLSIdx_.size() );
584  typename std::vector<int>::iterator p = curLSIdx_.begin();
585  for (; p<curLSIdx_.end(); ++p) {
586  // Check if this index is valid
587  if (*p != -1) {
588  // Check if any of the residuals are larger than the tolerance.
589  const int ensemble_size = StorageType::static_size;
590  bool all_converged = true;
591  for (int i=0; i<ensemble_size; ++i) {
592  if (!ensemble_converged[i]) {
593  if (testvector_[ *p ].coeff(i) > tolerance_) {
594  ++ensemble_iterations[i];
595  all_converged = false;
596  }
597  else if (testvector_[ *p ].coeff(i) <= tolerance_) {
598  ensemble_converged[i] = 1;
599  }
600  else {
601  // Throw an std::exception if a NaN is found.
602  status_ = Failed;
603  TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
604  }
605  }
606  }
607  if (all_converged) {
608  ind_[have] = *p;
609  have++;
610  }
611  }
612  }
613  ind_.resize(have);
614  int need = (quorum_ == -1) ? curNumRHS_: quorum_;
615  status_ = (have >= need) ? Passed : Failed;
616 
617  // Return the current status
618  return status_;
619 }
620 
621 template <class StorageType, class MV, class OP>
622 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::print(std::ostream& os, int indent) const
623 {
624  for (int j = 0; j < indent; j ++)
625  os << ' ';
626  printStatus(os, status_);
627  os << resFormStr();
628  if (status_==Undefined)
629  os << ", tol = " << tolerance_ << std::endl;
630  else {
631  os << std::endl;
632  if(showMaxResNormOnly_ && curBlksz_ > 1) {
633  const MagnitudeType maxRelRes = *std::max_element(
634  testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
635  );
636  for (int j = 0; j < indent + 13; j ++)
637  os << ' ';
638  os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
639  << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
640  }
641  else {
642  for ( int i=0; i<numrhs_; i++ ) {
643  for (int j = 0; j < indent + 13; j ++)
644  os << ' ';
645  os << "residual [ " << i << " ] = " << testvector_[ i ];
646  os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " " ) << tolerance_ << std::endl;
647  }
648  }
649  }
650  os << std::endl;
651 }
652 
653 template <class StorageType, class MV, class OP>
654 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::printStatus(std::ostream& os, StatusType type) const
655 {
656  os << std::left << std::setw(13) << std::setfill('.');
657  switch (type) {
658  case Passed:
659  os << "Converged";
660  break;
661  case Failed:
662  os << "Unconverged";
663  break;
664  case Undefined:
665  default:
666  os << "**";
667  break;
668  }
669  os << std::left << std::setfill(' ');
670  return;
671 }
672 
673 template <class StorageType, class MV, class OP>
674 StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
675 {
676  int i;
677  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
678  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
679  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
680  // Compute scaling term (done once for each block that's being solved)
681  if (firstcallCheckStatus_) {
682  //
683  // Get some current solver information.
684  //
685  firstcallCheckStatus_ = false;
686 
687  if (scaletype_== NormOfRHS) {
688  Teuchos::RCP<const MV> rhs = lp.getRHS();
689  numrhs_ = MVT::GetNumberVecs( *rhs );
690  scalevector_.resize( numrhs_ );
691  MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
692  }
693  else if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes) {
694  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
695  numrhs_ = MVT::GetNumberVecs( *init_res );
696  scalevector_.resize( numrhs_ );
697  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
698  }
699  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes) {
700  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
701  numrhs_ = MVT::GetNumberVecs( *init_res );
702  scalevector_.resize( numrhs_ );
703  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
704  }
705  else {
706  numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
707  }
708 
709  resvector_.resize( numrhs_ );
710  testvector_.resize( numrhs_ );
711 
712  curLSNum_ = lp.getLSNumber();
713  curLSIdx_ = lp.getLSIndex();
714  curBlksz_ = (int)curLSIdx_.size();
715  int validLS = 0;
716  for (i=0; i<curBlksz_; ++i) {
717  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
718  validLS++;
719  }
720  curNumRHS_ = validLS;
721  //
722  // Initialize the testvector.
723  for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
724 
725  // Return an error if the scaling is zero.
726  if (scalevalue_ == zero) {
727  return Failed;
728  }
729  }
730  return Undefined;
731 }
732 
733 } // end namespace Belos
734 
735 #endif /* BELOS_STATUS_TEST_RESNORM_H */
std::string description() const
Method to return description of the maximum iteration status test.
int curLSNum_
The current number of linear systems that have been loaded into the linear problem.
StatusType getStatus() const
Return the result of the most recent CheckStatus call.
NormType scalenormtype_
Type of norm to use on the scaling (OneNorm, TwoNorm, or InfNorm)
std::vector< int > ind_
Vector containing the indices for the vectors that passed the test.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
const std::vector< MagnitudeType > * getResNormValue() const
Returns the residual norm value, , computed in most recent call to CheckStatus.
bool getShowMaxResNormOnly()
Returns whether the only maximum residual norm is displayed when the print() method is called...
std::vector< int > curLSIdx_
The indices of the current number of right-hand sides being solved for.
std::vector< MagnitudeType > testvector_
Test std::vector = resvector_ / scalevector_.
An implementation of StatusTestResNorm using a family of residual norms.
std::vector< int > ensemble_iterations
The number of iterations at which point each ensemble component converges.
bool showMaxResNormOnly_
Determines if the entries for all of the residuals are shown or just the max.
int setShowMaxResNormOnly(bool showMaxResNormOnly)
Set whether the only maximum residual norm is displayed when the print() method is called...
const std::vector< MagnitudeType > * getScaledNormValue() const
Returns the scaled norm value, .
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
int quorum_
Number of residuals that must pass the convergence test before Passed is returned.
MagnitudeType getTolerance() const
Returns the value of the tolerance, , set in the constructor.
NormType resnormtype_
Type of norm to use on residual (OneNorm, TwoNorm, or InfNorm).
std::vector< int > convIndices()
Returns the std::vector containing the indices of the residuals that passed the test.
ScaleType scaletype_
Type of scaling to use (Norm of RHS, Norm of Initial Residual, None or User provided) ...
const std::vector< int > getEnsembleIterations() const
Returns number of ensemble iterations.
std::vector< int > ensemble_converged
Which ensemble components have converged.
Teuchos::RCP< MV > curSoln_
Most recent solution vector used by this status test.