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   bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
00208 
00210   std::vector<int> convIndices() { return ind_; }
00211 
00213   MagnitudeType getTolerance() const {return(tolerance_);};
00214   
00216   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00217 
00219   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00220 
00222   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00223 
00226   bool getLOADetected() const { return false; }
00227 
00229 
00230 
00233 
00239   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00241 
00244 
00246   std::string description() const
00247   {
00248     std::ostringstream oss;
00249     oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
00250     oss << ", tol = " << tolerance_;
00251     return oss.str();
00252   }
00254 
00255  protected:
00256 
00257  private:
00258 
00260 
00261 
00262   std::string resFormStr() const
00263   {
00264     std::ostringstream oss;
00265     oss << "(";
00266     oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00267     oss << ((restype_==Explicit) ? " Exp" : " Imp");
00268     oss << " Res Vec) ";
00269     if (scaletype_!=None)
00270       oss << "/ ";
00271     if (scaletype_==UserProvided)
00272       oss << " (User Scale)";
00273     else {
00274       oss << "(";
00275       oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00276       if (scaletype_==NormOfInitRes)
00277         oss << " Res0";
00278       else if (scaletype_==NormOfPrecInitRes)
00279         oss << " Prec Res0";
00280       else
00281         oss << " RHS ";
00282       oss << ")";
00283     }
00284     return oss.str();
00285   }
00286 
00288 
00290 
00291   
00293   MagnitudeType tolerance_;
00294 
00296   int quorum_;
00297 
00299   bool showMaxResNormOnly_;
00300  
00302   ResType restype_;
00303   
00305   NormType resnormtype_;
00306   
00308   ScaleType scaletype_;
00309   
00311   NormType scalenormtype_;
00312 
00314   MagnitudeType scalevalue_;
00315 
00317   std::vector<MagnitudeType> scalevector_;
00318   
00320   std::vector<MagnitudeType> resvector_;
00321 
00323   std::vector<MagnitudeType> testvector_;
00324 
00326   std::vector<int> ind_;
00327 
00329   Teuchos::RCP<MV> curSoln_;
00330   
00332   StatusType status_;
00333   
00335   int curBlksz_;
00336 
00338   int curNumRHS_;
00339 
00341   std::vector<int> curLSIdx_;
00342 
00344   int curLSNum_;
00345 
00347   int numrhs_;
00348 
00350   bool firstcallCheckStatus_;
00351 
00353   bool firstcallDefineResForm_;
00354 
00356   bool firstcallDefineScaleForm_;
00357 
00359 
00360 };
00361 
00362 template <class ScalarType, class MV, class OP>
00363 StatusTestGenResNorm<ScalarType,MV,OP>::StatusTestGenResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00364   : tolerance_(Tolerance),
00365     quorum_(quorum),
00366     showMaxResNormOnly_(showMaxResNormOnly),
00367     restype_(Implicit),
00368     resnormtype_(TwoNorm),  
00369     scaletype_(NormOfInitRes),
00370     scalenormtype_(TwoNorm),
00371     scalevalue_(1.0),
00372     status_(Undefined),
00373     curBlksz_(0),
00374     curLSNum_(0),
00375     numrhs_(0),
00376     firstcallCheckStatus_(true),
00377     firstcallDefineResForm_(true),
00378     firstcallDefineScaleForm_(true)
00379 {
00380   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00381   // the implicit residual std::vector.
00382 }
00383 
00384 template <class ScalarType, class MV, class OP>
00385 StatusTestGenResNorm<ScalarType,MV,OP>::~StatusTestGenResNorm() 
00386 {}
00387 
00388 template <class ScalarType, class MV, class OP>
00389 void StatusTestGenResNorm<ScalarType,MV,OP>::reset() 
00390 {
00391   status_ = Undefined;
00392   curBlksz_ = 0;
00393   curLSNum_ = 0;
00394   curLSIdx_.resize(0);
00395   numrhs_ = 0;
00396   ind_.resize(0);
00397   firstcallCheckStatus_ = true;
00398   curSoln_ = Teuchos::null;
00399 }
00400 
00401 template <class ScalarType, class MV, class OP>
00402 int StatusTestGenResNorm<ScalarType,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
00403 {    
00404   TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00405   "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
00406   firstcallDefineResForm_ = false;
00407     
00408   restype_ = TypeOfResidual;
00409   resnormtype_ = TypeOfNorm;
00410     
00411   return(0);
00412 }
00413 
00414 template <class ScalarType, class MV, class OP> 
00415 int StatusTestGenResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00416                                                          MagnitudeType ScaleValue )
00417 {
00418   TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00419   "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
00420   firstcallDefineScaleForm_ = false;
00421     
00422   scaletype_ = TypeOfScaling;
00423   scalenormtype_ = TypeOfNorm;
00424   scalevalue_ = ScaleValue;
00425     
00426   return(0);
00427 }
00428 
00429 template <class ScalarType, class MV, class OP>
00430 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00431 {
00432   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00433   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00434   // Compute scaling term (done once for each block that's being solved)
00435   if (firstcallCheckStatus_) {
00436     StatusType status = firstCallCheckStatusSetup(iSolver);
00437     if(status==Failed) {
00438       status_ = Failed;
00439       return(status_);
00440     }
00441   }
00442   //
00443   // This section computes the norm of the residual std::vector
00444   //
00445   if ( curLSNum_ != lp.getLSNumber() ) {
00446     //
00447     // We have moved on to the next rhs block
00448     //
00449     curLSNum_ = lp.getLSNumber();
00450     curLSIdx_ = lp.getLSIndex();
00451     curBlksz_ = (int)curLSIdx_.size();
00452     int validLS = 0;
00453     for (int i=0; i<curBlksz_; ++i) {
00454       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00455   validLS++;
00456     }
00457     curNumRHS_ = validLS; 
00458     curSoln_ = Teuchos::null;
00459     //
00460   } else {
00461     //
00462     // We are in the same rhs block, return if we are converged
00463     //
00464     if (status_==Passed) { return status_; }
00465   }
00466   if (restype_==Implicit) {
00467     //
00468     // get the native residual norms from the solver for this block of right-hand sides.
00469     // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00470     // Otherwise the native residual is assumed to be stored in the resvector_.
00471     //
00472     std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00473     RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00474     if ( residMV != Teuchos::null ) { 
00475       tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00476       MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00477       typename std::vector<int>::iterator p = curLSIdx_.begin();
00478       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00479         // Check if this index is valid
00480         if (*p != -1)  
00481           resvector_[*p] = tmp_resvector[i]; 
00482       }
00483     } else {
00484       typename std::vector<int>::iterator p = curLSIdx_.begin();
00485       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00486         // Check if this index is valid
00487         if (*p != -1)
00488           resvector_[*p] = tmp_resvector[i];
00489       }
00490     }
00491   }
00492   else if (restype_==Explicit) {
00493     //
00494     // Request the true residual for this block of right-hand sides.
00495     //
00496     RCP<MV> cur_update = iSolver->getCurrentUpdate();
00497     curSoln_ = lp.updateSolution( cur_update );
00498     RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
00499     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00500     std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
00501     MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
00502     typename std::vector<int>::iterator p = curLSIdx_.begin();
00503     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00504       // Check if this index is valid
00505       if (*p != -1)
00506         resvector_[*p] = tmp_resvector[i];
00507     }
00508   }
00509   //
00510   // Compute the new linear system residuals for testing.
00511   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00512   //
00513   if ( scalevector_.size() > 0 ) {
00514     typename std::vector<int>::iterator p = curLSIdx_.begin();
00515     for (; p<curLSIdx_.end(); ++p) {
00516       // Check if this index is valid
00517       if (*p != -1) {     
00518         // Scale the std::vector accordingly
00519         if ( scalevector_[ *p ] != zero ) {
00520           // Don't intentionally divide by zero.
00521           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00522         } else {
00523           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00524         }
00525       }
00526     }
00527   }
00528   else {
00529     typename std::vector<int>::iterator p = curLSIdx_.begin();
00530     for (; p<curLSIdx_.end(); ++p) {
00531       // Check if this index is valid
00532       if (*p != -1)     
00533         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00534     }
00535   } 
00536 
00537   // Check status of new linear system residuals and see if we have the quorum.
00538   int have = 0;
00539   ind_.resize( curLSIdx_.size() );
00540   typename std::vector<int>::iterator p = curLSIdx_.begin();
00541   for (; p<curLSIdx_.end(); ++p) {
00542     // Check if this index is valid
00543     if (*p != -1) {     
00544       // Check if any of the residuals are larger than the tolerance.
00545       if (testvector_[ *p ] > tolerance_) {
00546         // do nothing.
00547       } else if (testvector_[ *p ] <= tolerance_) { 
00548         ind_[have] = *p;
00549         have++;
00550       } else {
00551         // Throw an std::exception if a NaN is found.
00552         status_ = Failed;
00553         TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
00554       }
00555     }
00556   } 
00557   ind_.resize(have);
00558   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00559   status_ = (have >= need) ? Passed : Failed;
00560   
00561   // Return the current status
00562   return status_;
00563 }
00564 
00565 template <class ScalarType, class MV, class OP>
00566 void StatusTestGenResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00567 {
00568   for (int j = 0; j < indent; j ++)
00569     os << ' ';
00570   printStatus(os, status_);
00571   os << resFormStr();
00572   if (status_==Undefined)
00573     os << ", tol = " << tolerance_ << std::endl;
00574   else {
00575     os << std::endl;
00576     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00577       const MagnitudeType maxRelRes = *std::max_element(
00578         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00579         );
00580       for (int j = 0; j < indent + 13; j ++)
00581         os << ' ';
00582       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00583          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00584     }
00585     else {
00586       for ( int i=0; i<numrhs_; i++ ) {
00587         for (int j = 0; j < indent + 13; j ++)
00588           os << ' ';
00589         os << "residual [ " << i << " ] = " << testvector_[ i ];
00590         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00591       }
00592     }
00593   }
00594   os << std::endl;
00595 }
00596 
00597 template <class ScalarType, class MV, class OP>
00598 void StatusTestGenResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00599 {
00600   os << std::left << std::setw(13) << std::setfill('.');
00601   switch (type) {
00602   case  Passed:
00603     os << "Converged";
00604     break;
00605   case  Failed:
00606     os << "Unconverged";
00607     break;
00608   case  Undefined:
00609   default:
00610     os << "**";
00611     break;
00612   }
00613   os << std::left << std::setfill(' ');
00614     return;
00615 }
00616 
00617 template <class ScalarType, class MV, class OP>
00618 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00619 {
00620   int i;
00621   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00622   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00623   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00624   // Compute scaling term (done once for each block that's being solved)
00625   if (firstcallCheckStatus_) {
00626     //
00627     // Get some current solver information.
00628     //
00629     firstcallCheckStatus_ = false;
00630 
00631     if (scaletype_== NormOfRHS) {
00632       RCP<const MV> rhs = lp.getRHS();
00633       numrhs_ = MVT::GetNumberVecs( *rhs );
00634       scalevector_.resize( numrhs_ );
00635       resvector_.resize( numrhs_ ); 
00636       testvector_.resize( numrhs_ );
00637       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00638     }
00639     else if (scaletype_==NormOfInitRes) {
00640       RCP<const MV> init_res = lp.getInitResVec();
00641       numrhs_ = MVT::GetNumberVecs( *init_res );
00642       scalevector_.resize( numrhs_ );
00643       resvector_.resize( numrhs_ ); 
00644       testvector_.resize( numrhs_ );
00645       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00646     }
00647     else if (scaletype_==NormOfPrecInitRes) {
00648       RCP<const MV> init_res = lp.getInitPrecResVec();
00649       numrhs_ = MVT::GetNumberVecs( *init_res );
00650       scalevector_.resize( numrhs_ );
00651       resvector_.resize( numrhs_ ); 
00652       testvector_.resize( numrhs_ );
00653       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00654     }
00655 
00656     curLSNum_ = lp.getLSNumber();
00657     curLSIdx_ = lp.getLSIndex();
00658     curBlksz_ = (int)curLSIdx_.size();
00659     int validLS = 0;
00660     for (i=0; i<curBlksz_; ++i) {
00661       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00662         validLS++;
00663     }
00664     curNumRHS_ = validLS;
00665     //
00666     // Initialize the testvector.
00667     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00668 
00669     // Return an error if the scaling is zero.
00670     if (scalevalue_ == zero) {
00671       return Failed;
00672     }
00673   }
00674   return Undefined;
00675 }
00676 
00677 } // end namespace Belos
00678 
00679 #endif /* BELOS_STATUS_TEST_RESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:07 2011 for Belos by  doxygen 1.6.3