BelosStatusTestGenResNorm.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_STATUS_TEST_GEN_RESNORM_H
00030 #define BELOS_STATUS_TEST_GEN_RESNORM_H
00031 
00037 #include "BelosStatusTestResNorm.hpp"
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosMultiVecTraits.hpp"
00040 
00063 namespace Belos {
00064 
00065 template <class ScalarType, class MV, class OP>
00066 class StatusTestGenResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00067 
00068  public:
00069 
00070   // Convenience typedefs
00071   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00072   typedef typename SCT::magnitudeType MagnitudeType;
00073   typedef MultiVecTraits<ScalarType,MV>  MVT;
00074 
00076 
00077 
00081   enum ResType {Implicit, 
00082     Explicit  
00084   };
00085 
00087 
00089 
00090 
00091 
00104   StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
00105 
00107   virtual ~StatusTestGenResNorm();
00109 
00111 
00112 
00114 
00121   int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
00122   
00124 
00144   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00145 
00147 
00150   int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
00151 
00154   int setQuorum(int quorum) {quorum_ = quorum; return(0);}
00155 
00157   int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
00158 
00160 
00162 
00163 
00164 
00170   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00171 
00173   StatusType getStatus() const {return(status_);};
00175 
00177 
00178  
00180   void reset();
00181 
00183 
00185 
00186 
00188   void print(std::ostream& os, int indent = 0) const;
00189 
00191   void printStatus(std::ostream& os, StatusType type) const; 
00193 
00195 
00196 
00200   Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
00201 
00204   int getQuorum() const { return quorum_; }
00205 
00207   std::vector<int> convIndices() { return ind_; }
00208 
00210   MagnitudeType getTolerance() const {return(tolerance_);};
00211   
00213   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00214 
00216   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00217 
00219   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00220 
00223   bool getLOADetected() const { return false; }
00224 
00226 
00227 
00230 
00236   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00238 
00241 
00243   std::string description() const
00244   {
00245     std::ostringstream oss;
00246     oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
00247     oss << ", tol = " << tolerance_;
00248     return oss.str();
00249   }
00251 
00252  protected:
00253 
00254  private:
00255 
00257 
00258 
00259   std::string resFormStr() const
00260   {
00261     std::ostringstream oss;
00262     oss << "(";
00263     oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00264     oss << ((restype_==Explicit) ? " Exp" : " Imp");
00265     oss << " Res Vec) ";
00266     if (scaletype_!=None)
00267       oss << "/ ";
00268     if (scaletype_==UserProvided)
00269       oss << " (User Scale)";
00270     else {
00271       oss << "(";
00272       oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00273       if (scaletype_==NormOfInitRes)
00274         oss << " Res0";
00275       else if (scaletype_==NormOfPrecInitRes)
00276         oss << " Prec Res0";
00277       else
00278         oss << " RHS ";
00279       oss << ")";
00280     }
00281     return oss.str();
00282   }
00283 
00285 
00287 
00288   
00290   MagnitudeType tolerance_;
00291 
00293   int quorum_;
00294 
00296   bool showMaxResNormOnly_;
00297  
00299   ResType restype_;
00300   
00302   NormType resnormtype_;
00303   
00305   ScaleType scaletype_;
00306   
00308   NormType scalenormtype_;
00309 
00311   MagnitudeType scalevalue_;
00312 
00314   std::vector<MagnitudeType> scalevector_;
00315   
00317   std::vector<MagnitudeType> resvector_;
00318 
00320   std::vector<MagnitudeType> testvector_;
00321 
00323   std::vector<int> ind_;
00324 
00326   Teuchos::RCP<MV> curSoln_;
00327   
00329   StatusType status_;
00330   
00332   int curBlksz_;
00333 
00335   int curNumRHS_;
00336 
00338   std::vector<int> curLSIdx_;
00339 
00341   int curLSNum_;
00342 
00344   int numrhs_;
00345 
00347   bool firstcallCheckStatus_;
00348 
00350   bool firstcallDefineResForm_;
00351 
00353   bool firstcallDefineScaleForm_;
00354 
00356 
00357 };
00358 
00359 template <class ScalarType, class MV, class OP>
00360 StatusTestGenResNorm<ScalarType,MV,OP>::StatusTestGenResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00361   : tolerance_(Tolerance),
00362     quorum_(quorum),
00363     showMaxResNormOnly_(showMaxResNormOnly),
00364     restype_(Implicit),
00365     resnormtype_(TwoNorm),  
00366     scaletype_(NormOfInitRes),
00367     scalenormtype_(TwoNorm),
00368     scalevalue_(1.0),
00369     status_(Undefined),
00370     curBlksz_(0),
00371     curLSNum_(0),
00372     numrhs_(0),
00373     firstcallCheckStatus_(true),
00374     firstcallDefineResForm_(true),
00375     firstcallDefineScaleForm_(true)
00376 {
00377   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00378   // the implicit residual std::vector.
00379 }
00380 
00381 template <class ScalarType, class MV, class OP>
00382 StatusTestGenResNorm<ScalarType,MV,OP>::~StatusTestGenResNorm() 
00383 {}
00384 
00385 template <class ScalarType, class MV, class OP>
00386 void StatusTestGenResNorm<ScalarType,MV,OP>::reset() 
00387 {
00388   status_ = Undefined;
00389   curBlksz_ = 0;
00390   curLSNum_ = 0;
00391   curLSIdx_.resize(0);
00392   numrhs_ = 0;
00393   ind_.resize(0);
00394   firstcallCheckStatus_ = true;
00395   curSoln_ = Teuchos::null;
00396 }
00397 
00398 template <class ScalarType, class MV, class OP>
00399 int StatusTestGenResNorm<ScalarType,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
00400 {    
00401   TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00402   "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
00403   firstcallDefineResForm_ = false;
00404     
00405   restype_ = TypeOfResidual;
00406   resnormtype_ = TypeOfNorm;
00407     
00408   return(0);
00409 }
00410 
00411 template <class ScalarType, class MV, class OP> 
00412 int StatusTestGenResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00413                                                          MagnitudeType ScaleValue )
00414 {
00415   TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00416   "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
00417   firstcallDefineScaleForm_ = false;
00418     
00419   scaletype_ = TypeOfScaling;
00420   scalenormtype_ = TypeOfNorm;
00421   scalevalue_ = ScaleValue;
00422     
00423   return(0);
00424 }
00425 
00426 template <class ScalarType, class MV, class OP>
00427 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00428 {
00429   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00430   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00431   // Compute scaling term (done once for each block that's being solved)
00432   if (firstcallCheckStatus_) {
00433     StatusType status = firstCallCheckStatusSetup(iSolver);
00434     if(status==Failed) {
00435       status_ = Failed;
00436       return(status_);
00437     }
00438   }
00439   //
00440   // This section computes the norm of the residual std::vector
00441   //
00442   if ( curLSNum_ != lp.getLSNumber() ) {
00443     //
00444     // We have moved on to the next rhs block
00445     //
00446     curLSNum_ = lp.getLSNumber();
00447     curLSIdx_ = lp.getLSIndex();
00448     curBlksz_ = (int)curLSIdx_.size();
00449     int validLS = 0;
00450     for (int i=0; i<curBlksz_; ++i) {
00451       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00452   validLS++;
00453     }
00454     curNumRHS_ = validLS; 
00455     curSoln_ = Teuchos::null;
00456     //
00457   } else {
00458     //
00459     // We are in the same rhs block, return if we are converged
00460     //
00461     if (status_==Passed) { return status_; }
00462   }
00463   if (restype_==Implicit) {
00464     //
00465     // get the native residual norms from the solver for this block of right-hand sides.
00466     // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00467     // Otherwise the native residual is assumed to be stored in the resvector_.
00468     //
00469     std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00470     RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00471     if ( residMV != Teuchos::null ) { 
00472       tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00473       MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00474       typename std::vector<int>::iterator p = curLSIdx_.begin();
00475       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00476         // Check if this index is valid
00477         if (*p != -1)  
00478           resvector_[*p] = tmp_resvector[i]; 
00479       }
00480     } else {
00481       typename std::vector<int>::iterator p = curLSIdx_.begin();
00482       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00483         // Check if this index is valid
00484         if (*p != -1)
00485           resvector_[*p] = tmp_resvector[i];
00486       }
00487     }
00488   }
00489   else if (restype_==Explicit) {
00490     //
00491     // Request the true residual for this block of right-hand sides.
00492     //
00493     RCP<MV> cur_update = iSolver->getCurrentUpdate();
00494     curSoln_ = lp.updateSolution( cur_update );
00495     RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
00496     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00497     std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
00498     MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
00499     typename std::vector<int>::iterator p = curLSIdx_.begin();
00500     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00501       // Check if this index is valid
00502       if (*p != -1)
00503         resvector_[*p] = tmp_resvector[i];
00504     }
00505   }
00506   //
00507   // Compute the new linear system residuals for testing.
00508   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00509   //
00510   if ( scalevector_.size() > 0 ) {
00511     typename std::vector<int>::iterator p = curLSIdx_.begin();
00512     for (; p<curLSIdx_.end(); ++p) {
00513       // Check if this index is valid
00514       if (*p != -1) {     
00515         // Scale the std::vector accordingly
00516         if ( scalevector_[ *p ] != zero ) {
00517           // Don't intentionally divide by zero.
00518           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00519         } else {
00520           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00521         }
00522       }
00523     }
00524   }
00525   else {
00526     typename std::vector<int>::iterator p = curLSIdx_.begin();
00527     for (; p<curLSIdx_.end(); ++p) {
00528       // Check if this index is valid
00529       if (*p != -1)     
00530         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00531     }
00532   } 
00533 
00534   // Check status of new linear system residuals and see if we have the quorum.
00535   int have = 0;
00536   ind_.resize( curLSIdx_.size() );
00537   typename std::vector<int>::iterator p = curLSIdx_.begin();
00538   for (; p<curLSIdx_.end(); ++p) {
00539     // Check if this index is valid
00540     if (*p != -1) {     
00541       // Check if any of the residuals are larger than the tolerance.
00542       if (testvector_[ *p ] > tolerance_) {
00543         // do nothing.
00544       } else if (testvector_[ *p ] <= tolerance_) { 
00545         ind_[have] = *p;
00546         have++;
00547       } else {
00548         // Throw an std::exception if a NaN is found.
00549         status_ = Failed;
00550         TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
00551       }
00552     }
00553   } 
00554   ind_.resize(have);
00555   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00556   status_ = (have >= need) ? Passed : Failed;
00557   
00558   // Return the current status
00559   return status_;
00560 }
00561 
00562 template <class ScalarType, class MV, class OP>
00563 void StatusTestGenResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00564 {
00565   for (int j = 0; j < indent; j ++)
00566     os << ' ';
00567   printStatus(os, status_);
00568   os << resFormStr();
00569   if (status_==Undefined)
00570     os << ", tol = " << tolerance_ << std::endl;
00571   else {
00572     os << std::endl;
00573     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00574       const MagnitudeType maxRelRes = *std::max_element(
00575         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00576         );
00577       for (int j = 0; j < indent + 13; j ++)
00578         os << ' ';
00579       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00580          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00581     }
00582     else {
00583       for ( int i=0; i<numrhs_; i++ ) {
00584         for (int j = 0; j < indent + 13; j ++)
00585           os << ' ';
00586         os << "residual [ " << i << " ] = " << testvector_[ i ];
00587         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00588       }
00589     }
00590   }
00591   os << std::endl;
00592 }
00593 
00594 template <class ScalarType, class MV, class OP>
00595 void StatusTestGenResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00596 {
00597   os << std::left << std::setw(13) << std::setfill('.');
00598   switch (type) {
00599   case  Passed:
00600     os << "Converged";
00601     break;
00602   case  Failed:
00603     os << "Unconverged";
00604     break;
00605   case  Undefined:
00606   default:
00607     os << "**";
00608     break;
00609   }
00610   os << std::left << std::setfill(' ');
00611     return;
00612 }
00613 
00614 template <class ScalarType, class MV, class OP>
00615 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00616 {
00617   int i;
00618   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00619   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00620   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00621   // Compute scaling term (done once for each block that's being solved)
00622   if (firstcallCheckStatus_) {
00623     //
00624     // Get some current solver information.
00625     //
00626     firstcallCheckStatus_ = false;
00627 
00628     if (scaletype_== NormOfRHS) {
00629       RCP<const MV> rhs = lp.getRHS();
00630       numrhs_ = MVT::GetNumberVecs( *rhs );
00631       scalevector_.resize( numrhs_ );
00632       resvector_.resize( numrhs_ ); 
00633       testvector_.resize( numrhs_ );
00634       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00635     }
00636     else if (scaletype_==NormOfInitRes) {
00637       RCP<const MV> init_res = lp.getInitResVec();
00638       numrhs_ = MVT::GetNumberVecs( *init_res );
00639       scalevector_.resize( numrhs_ );
00640       resvector_.resize( numrhs_ ); 
00641       testvector_.resize( numrhs_ );
00642       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00643     }
00644     else if (scaletype_==NormOfPrecInitRes) {
00645       RCP<const MV> init_res = lp.getInitPrecResVec();
00646       numrhs_ = MVT::GetNumberVecs( *init_res );
00647       scalevector_.resize( numrhs_ );
00648       resvector_.resize( numrhs_ ); 
00649       testvector_.resize( numrhs_ );
00650       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00651     }
00652 
00653     curLSNum_ = lp.getLSNumber();
00654     curLSIdx_ = lp.getLSIndex();
00655     curBlksz_ = (int)curLSIdx_.size();
00656     int validLS = 0;
00657     for (i=0; i<curBlksz_; ++i) {
00658       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00659         validLS++;
00660     }
00661     curNumRHS_ = validLS;
00662     //
00663     // Initialize the testvector.
00664     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00665 
00666     // Return an error if the scaling is zero.
00667     if (scalevalue_ == zero) {
00668       return Failed;
00669     }
00670   }
00671   return Undefined;
00672 }
00673 
00674 } // end namespace Belos
00675 
00676 #endif /* BELOS_STATUS_TEST_RESNORM_H */

Generated on Wed May 12 21:30:09 2010 for Belos by  doxygen 1.4.7