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 
167  NormType getResNormType() {return resnormtype_;}
168 
170 
173  int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
174 
177  int setQuorum(int quorum) {quorum_ = quorum; return(0);}
178 
180  int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
181 
183 
185 
186 
193  StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
194 
196  StatusType getStatus() const {return(status_);}
198 
200 
201 
203  void reset();
204 
206 
208 
209 
211  void print(std::ostream& os, int indent = 0) const;
212 
214  void printStatus(std::ostream& os, StatusType type) const;
216 
218 
219 
223  Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
224 
227  int getQuorum() const { return quorum_; }
228 
230  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
231 
233  std::vector<int> convIndices() { return ind_; }
234 
236  MagnitudeType getTolerance() const {return(tolerance_);}
237 
239  const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);}
240 
242  const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);}
243 
245  const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);}
246 
249  bool getLOADetected() const { return false; }
250 
252  const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
253 
255 
256 
259 
265  StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
267 
270 
272  std::string description() const
273  {
274  std::ostringstream oss;
275  oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
276  oss << ", tol = " << tolerance_;
277  return oss.str();
278  }
280 
281  protected:
282 
283  private:
284 
286 
287 
288  std::string resFormStr() const
289  {
290  std::ostringstream oss;
291  oss << "(";
292  oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
293  oss << ((restype_==Explicit) ? " Exp" : " Imp");
294  oss << " Res Vec) ";
295 
296  // If there is no residual scaling, return current string.
297  if (scaletype_!=None)
298  {
299  // Insert division sign.
300  oss << "/ ";
301 
302  // Determine output string for scaling, if there is any.
303  if (scaletype_==UserProvided)
304  oss << " (User Scale)";
305  else {
306  oss << "(";
307  oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
308  if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes)
309  oss << " Res0";
310  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes)
311  oss << " Prec Res0";
312  else
313  oss << " RHS ";
314  oss << ")";
315  }
316  }
317 
318  return oss.str();
319  }
320 
322 
324 
325 
328 
330  int quorum_;
331 
334 
337 
339  NormType resnormtype_;
340 
342  ScaleType scaletype_;
343 
345  NormType scalenormtype_;
346 
349 
351  std::vector<MagnitudeType> scalevector_;
352 
354  std::vector<MagnitudeType> resvector_;
355 
357  std::vector<MagnitudeType> testvector_;
358 
360  std::vector<int> ind_;
361 
363  Teuchos::RCP<MV> curSoln_;
364 
366  StatusType status_;
367 
370 
373 
375  std::vector<int> curLSIdx_;
376 
379 
381  int numrhs_;
382 
385 
388 
391 
393  std::vector<int> ensemble_converged;
394 
396  std::vector<int> ensemble_iterations;
397 
399 
400 };
401 
402 template <class StorageType, class MV, class OP>
404 StatusTestGenResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
405  : tolerance_(Tolerance),
406  quorum_(quorum),
407  showMaxResNormOnly_(showMaxResNormOnly),
408  restype_(Implicit),
409  resnormtype_(TwoNorm),
410  scaletype_(NormOfInitRes),
411  scalenormtype_(TwoNorm),
412  scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
413  status_(Undefined),
414  curBlksz_(0),
415  curNumRHS_(0),
416  curLSNum_(0),
417  numrhs_(0),
418  firstcallCheckStatus_(true),
419  firstcallDefineResForm_(true),
420  firstcallDefineScaleForm_(true),
421  ensemble_converged(StorageType::static_size, 0),
422  ensemble_iterations(StorageType::static_size, 0)
423 {
424  // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
425  // the implicit residual std::vector.
426 }
427 
428 template <class StorageType, class MV, class OP>
429 StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::~StatusTestGenResNorm()
430 {}
431 
432 template <class StorageType, class MV, class OP>
433 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::reset()
434 {
435  status_ = Undefined;
436  curBlksz_ = 0;
437  curLSNum_ = 0;
438  curLSIdx_.resize(0);
439  numrhs_ = 0;
440  ind_.resize(0);
441  firstcallCheckStatus_ = true;
442  curSoln_ = Teuchos::null;
443  ensemble_converged = std::vector<int>(StorageType::static_size, 0);
444  ensemble_iterations = std::vector<int>(StorageType::static_size, 0);
445 }
446 
447 template <class StorageType, class MV, class OP>
448 int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
449 {
450  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
451  "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
452  firstcallDefineResForm_ = false;
453 
454  restype_ = TypeOfResidual;
455  resnormtype_ = TypeOfNorm;
456 
457  return(0);
458 }
459 
460 template <class StorageType, class MV, class OP>
461 int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
462  MagnitudeType ScaleValue )
463 {
464  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
465  "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
466  firstcallDefineScaleForm_ = false;
467 
468  scaletype_ = TypeOfScaling;
469  scalenormtype_ = TypeOfNorm;
470  scalevalue_ = ScaleValue;
471 
472  return(0);
473 }
474 
475 template <class StorageType, class MV, class OP>
476 StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
477 {
478  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
479  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
480  // Compute scaling term (done once for each block that's being solved)
481  if (firstcallCheckStatus_) {
482  StatusType status = firstCallCheckStatusSetup(iSolver);
483  if(status==Failed) {
484  status_ = Failed;
485  return(status_);
486  }
487  }
488  //
489  // This section computes the norm of the residual std::vector
490  //
491  if ( curLSNum_ != lp.getLSNumber() ) {
492  //
493  // We have moved on to the next rhs block
494  //
495  curLSNum_ = lp.getLSNumber();
496  curLSIdx_ = lp.getLSIndex();
497  curBlksz_ = (int)curLSIdx_.size();
498  int validLS = 0;
499  for (int i=0; i<curBlksz_; ++i) {
500  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
501  validLS++;
502  }
503  curNumRHS_ = validLS;
504  curSoln_ = Teuchos::null;
505  //
506  } else {
507  //
508  // We are in the same rhs block, return if we are converged
509  //
510  if (status_==Passed) { return status_; }
511  }
512  if (restype_==Implicit) {
513  //
514  // get the native residual norms from the solver for this block of right-hand sides.
515  // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
516  // Otherwise the native residual is assumed to be stored in the resvector_.
517  //
518  std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
519  Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );
520  if ( residMV != Teuchos::null ) {
521  tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
522  MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );
523  typename std::vector<int>::iterator p = curLSIdx_.begin();
524  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
525  // Check if this index is valid
526  if (*p != -1)
527  resvector_[*p] = tmp_resvector[i];
528  }
529  } else {
530  typename std::vector<int>::iterator p = curLSIdx_.begin();
531  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
532  // Check if this index is valid
533  if (*p != -1)
534  resvector_[*p] = tmp_resvector[i];
535  }
536  }
537  }
538  else if (restype_==Explicit) {
539  //
540  // Request the true residual for this block of right-hand sides.
541  //
542  Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
543  curSoln_ = lp.updateSolution( cur_update );
544  Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
545  lp.computeCurrResVec( &*cur_res, &*curSoln_ );
546  std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
547  MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
548  typename std::vector<int>::iterator p = curLSIdx_.begin();
549  for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
550  // Check if this index is valid
551  if (*p != -1)
552  resvector_[*p] = tmp_resvector[i];
553  }
554  }
555  //
556  // Compute the new linear system residuals for testing.
557  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
558  //
559  if ( scalevector_.size() > 0 ) {
560  typename std::vector<int>::iterator p = curLSIdx_.begin();
561  for (; p<curLSIdx_.end(); ++p) {
562  // Check if this index is valid
563  if (*p != -1) {
564  // Scale the std::vector accordingly
565  if ( scalevector_[ *p ] != zero ) {
566  // Don't intentionally divide by zero.
567  testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
568  } else {
569  testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
570  }
571  }
572  }
573  }
574  else {
575  typename std::vector<int>::iterator p = curLSIdx_.begin();
576  for (; p<curLSIdx_.end(); ++p) {
577  // Check if this index is valid
578  if (*p != -1)
579  testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
580  }
581  }
582 
583  // Check status of new linear system residuals and see if we have the quorum.
584  int have = 0;
585  ind_.resize( curLSIdx_.size() );
586  typename std::vector<int>::iterator p = curLSIdx_.begin();
587  for (; p<curLSIdx_.end(); ++p) {
588  // Check if this index is valid
589  if (*p != -1) {
590  // Check if any of the residuals are larger than the tolerance.
591  const int ensemble_size = StorageType::static_size;
592  bool all_converged = true;
593  for (int i=0; i<ensemble_size; ++i) {
594  if (!ensemble_converged[i]) {
595  if (testvector_[ *p ].coeff(i) > tolerance_) {
596  ++ensemble_iterations[i];
597  all_converged = false;
598  }
599  else if (testvector_[ *p ].coeff(i) <= tolerance_) {
600  ensemble_converged[i] = 1;
601  }
602  else {
603  // Throw an std::exception if a NaN is found.
604  status_ = Failed;
605  TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
606  }
607  }
608  }
609  if (all_converged) {
610  ind_[have] = *p;
611  have++;
612  }
613  }
614  }
615  ind_.resize(have);
616  int need = (quorum_ == -1) ? curNumRHS_: quorum_;
617  status_ = (have >= need) ? Passed : Failed;
618 
619  // Return the current status
620  return status_;
621 }
622 
623 template <class StorageType, class MV, class OP>
624 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::print(std::ostream& os, int indent) const
625 {
626  for (int j = 0; j < indent; j ++)
627  os << ' ';
628  printStatus(os, status_);
629  os << resFormStr();
630  if (status_==Undefined)
631  os << ", tol = " << tolerance_ << std::endl;
632  else {
633  os << std::endl;
634  if(showMaxResNormOnly_ && curBlksz_ > 1) {
635  const MagnitudeType maxRelRes = *std::max_element(
636  testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
637  );
638  for (int j = 0; j < indent + 13; j ++)
639  os << ' ';
640  os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
641  << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
642  }
643  else {
644  for ( int i=0; i<numrhs_; i++ ) {
645  for (int j = 0; j < indent + 13; j ++)
646  os << ' ';
647  os << "residual [ " << i << " ] = " << testvector_[ i ];
648  os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " " ) << tolerance_ << std::endl;
649  }
650  }
651  }
652  os << std::endl;
653 }
654 
655 template <class StorageType, class MV, class OP>
656 void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::printStatus(std::ostream& os, StatusType type) const
657 {
658  os << std::left << std::setw(13) << std::setfill('.');
659  switch (type) {
660  case Passed:
661  os << "Converged";
662  break;
663  case Failed:
664  os << "Unconverged";
665  break;
666  case Undefined:
667  default:
668  os << "**";
669  break;
670  }
671  os << std::left << std::setfill(' ');
672  return;
673 }
674 
675 template <class StorageType, class MV, class OP>
676 StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
677 {
678  int i;
679  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
680  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
681  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
682  // Compute scaling term (done once for each block that's being solved)
683  if (firstcallCheckStatus_) {
684  //
685  // Get some current solver information.
686  //
687  firstcallCheckStatus_ = false;
688 
689  if (scaletype_== NormOfRHS) {
690  Teuchos::RCP<const MV> rhs = lp.getRHS();
691  numrhs_ = MVT::GetNumberVecs( *rhs );
692  scalevector_.resize( numrhs_ );
693  MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
694  }
695  else if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes) {
696  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
697  numrhs_ = MVT::GetNumberVecs( *init_res );
698  scalevector_.resize( numrhs_ );
699  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
700  }
701  else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes) {
702  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
703  numrhs_ = MVT::GetNumberVecs( *init_res );
704  scalevector_.resize( numrhs_ );
705  MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
706  }
707  else {
708  numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
709  }
710 
711  resvector_.resize( numrhs_ );
712  testvector_.resize( numrhs_ );
713 
714  curLSNum_ = lp.getLSNumber();
715  curLSIdx_ = lp.getLSIndex();
716  curBlksz_ = (int)curLSIdx_.size();
717  int validLS = 0;
718  for (i=0; i<curBlksz_; ++i) {
719  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
720  validLS++;
721  }
722  curNumRHS_ = validLS;
723  //
724  // Initialize the testvector.
725  for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
726 
727  // Return an error if the scaling is zero.
728  if (scalevalue_ == zero) {
729  return Failed;
730  }
731  }
732  return Undefined;
733 }
734 
735 } // end namespace Belos
736 
737 #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.