BelosStatusTestImpResNorm.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_IMPRESNORM_H
00030 #define BELOS_STATUS_TEST_IMPRESNORM_H
00031 
00037 #include "BelosStatusTestResNorm.hpp"
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosMultiVecTraits.hpp"
00040 
00062 namespace Belos {
00063 
00064 template <class ScalarType, class MV, class OP>
00065 class StatusTestImpResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00066 
00067  public:
00068 
00069   // Convenience typedefs
00070   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00071   typedef typename SCT::magnitudeType MagnitudeType;
00072   typedef MultiVecTraits<ScalarType,MV>  MVT;
00073 
00075 
00076 
00077 
00090   StatusTestImpResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
00091 
00093   virtual ~StatusTestImpResNorm();
00095 
00097 
00098 
00100 
00106   int defineResForm( NormType TypeOfNorm);
00107   
00109 
00129   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00130 
00132 
00135   int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
00136 
00139   int setQuorum(int quorum) {quorum_ = quorum; return(0);}
00140 
00142   int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
00143 
00145 
00147 
00148 
00149 
00155   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00156 
00158   StatusType getStatus() const {return(status_);};
00160 
00162 
00163  
00165   void reset();
00166 
00168 
00170 
00171 
00173   void print(std::ostream& os, int indent = 0) const;
00174 
00176   void printStatus(std::ostream& os, StatusType type) const; 
00178 
00180 
00181 
00183   Teuchos::RCP<MV> getSolution() { return curSoln_; }
00184 
00187   int getQuorum() const { return quorum_; }
00188 
00190   bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
00191 
00193   std::vector<int> convIndices() { return ind_; }
00194   
00196   MagnitudeType getTolerance() const {return(tolerance_);};
00197   
00199   MagnitudeType getCurrTolerance() const {return(currTolerance_);};
00200   
00202   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00203 
00205   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00206 
00208   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00209 
00211   bool getLOADetected() const { return lossDetected_; }
00213 
00214 
00217 
00223   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00225 
00228 
00230   std::string description() const
00231   {
00232     std::ostringstream oss;
00233     oss << "Belos::StatusTestImpResNorm<>: " << resFormStr();
00234     oss << ", tol = " << tolerance_;
00235     return oss.str();
00236   }
00238 
00239  protected:
00240 
00241  private:
00242  
00244 
00245 
00246   std::string resFormStr() const
00247   {
00248     std::ostringstream oss;
00249     oss << "(";
00250     oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00251     oss << " Res Vec) ";
00252     if (scaletype_!=None)
00253       oss << "/ ";
00254     if (scaletype_==UserProvided)
00255       oss << " (User Scale)";
00256     else {
00257       oss << "(";
00258       oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00259       if (scaletype_==NormOfInitRes)
00260         oss << " Res0";
00261       else if (scaletype_==NormOfPrecInitRes)
00262         oss << " Prec Res0";
00263       else
00264         oss << " RHS ";
00265       oss << ")";
00266     }
00267     return oss.str();
00268   }
00269 
00271 
00272   
00274   MagnitudeType tolerance_, currTolerance_;
00275 
00277   int quorum_;
00278 
00280   bool showMaxResNormOnly_;
00281  
00283   NormType resnormtype_;
00284   
00286   ScaleType scaletype_;
00287   
00289   NormType scalenormtype_;
00290 
00292   MagnitudeType scalevalue_;
00293 
00295   std::vector<MagnitudeType> scalevector_;
00296   
00298   std::vector<MagnitudeType> resvector_;
00299 
00301   std::vector<MagnitudeType> testvector_;
00302 
00304   Teuchos::RCP<MV> curSoln_; 
00305 
00307   std::vector<int> ind_;
00308   
00310   StatusType status_;
00311   
00313   int curBlksz_;
00314 
00316   int curNumRHS_;
00317 
00319   std::vector<int> curLSIdx_;
00320 
00322   int curLSNum_;
00323 
00325   int numrhs_;
00326 
00328   bool firstcallCheckStatus_;
00329 
00331   bool firstcallDefineResForm_;
00332 
00334   bool firstcallDefineScaleForm_;
00335 
00337   bool lossDetected_;
00339 
00340 };
00341 
00342 template <class ScalarType, class MV, class OP>
00343 StatusTestImpResNorm<ScalarType,MV,OP>::StatusTestImpResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00344   : tolerance_(Tolerance),
00345     currTolerance_(Tolerance),
00346     quorum_(quorum),
00347     showMaxResNormOnly_(showMaxResNormOnly),
00348     resnormtype_(TwoNorm),  
00349     scaletype_(NormOfInitRes),
00350     scalenormtype_(TwoNorm),
00351     scalevalue_(1.0),
00352     status_(Undefined),
00353     curBlksz_(0),
00354     curLSNum_(0),
00355     numrhs_(0),
00356     firstcallCheckStatus_(true),
00357     firstcallDefineResForm_(true),
00358     firstcallDefineScaleForm_(true),
00359     lossDetected_(false)
00360 {
00361   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00362   // the implicit residual vector.
00363 }
00364 
00365 template <class ScalarType, class MV, class OP>
00366 StatusTestImpResNorm<ScalarType,MV,OP>::~StatusTestImpResNorm() 
00367 {}
00368 
00369 template <class ScalarType, class MV, class OP>
00370 void StatusTestImpResNorm<ScalarType,MV,OP>::reset() 
00371 {
00372   status_ = Undefined;
00373   curBlksz_ = 0;
00374   curLSNum_ = 0;
00375   curLSIdx_.resize(0);
00376   numrhs_ = 0;
00377   ind_.resize(0);
00378   currTolerance_ = tolerance_;
00379   firstcallCheckStatus_ = true;
00380   lossDetected_ = false;
00381   curSoln_ = Teuchos::null;
00382 }
00383 
00384 template <class ScalarType, class MV, class OP>
00385 int StatusTestImpResNorm<ScalarType,MV,OP>::defineResForm( NormType TypeOfNorm )
00386 {    
00387   TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00388   "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
00389   firstcallDefineResForm_ = false;
00390     
00391   resnormtype_ = TypeOfNorm;
00392     
00393   return(0);
00394 }
00395 
00396 template <class ScalarType, class MV, class OP> 
00397 int StatusTestImpResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00398                                                          MagnitudeType ScaleValue )
00399 {
00400   TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00401   "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
00402   firstcallDefineScaleForm_ = false;
00403     
00404   scaletype_ = TypeOfScaling;
00405   scalenormtype_ = TypeOfNorm;
00406   scalevalue_ = ScaleValue;
00407     
00408   return(0);
00409 }
00410 
00411 template <class ScalarType, class MV, class OP>
00412 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00413 {
00414   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00415   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00416   // Compute scaling term (done once for each block that's being solved)
00417   if (firstcallCheckStatus_) {
00418     StatusType status = firstCallCheckStatusSetup(iSolver);
00419     if(status==Failed) {
00420       status_ = Failed;
00421       return(status_);
00422     }
00423   }
00424   //
00425   // This section computes the norm of the residual vector
00426   //
00427   if ( curLSNum_ != lp.getLSNumber() ) {
00428     //
00429     // We have moved on to the next rhs block
00430     //
00431     curLSNum_ = lp.getLSNumber();
00432     curLSIdx_ = lp.getLSIndex();
00433     curBlksz_ = (int)curLSIdx_.size();
00434     int validLS = 0;
00435     for (int i=0; i<curBlksz_; ++i) {
00436       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00437   validLS++;
00438     }
00439     curNumRHS_ = validLS; 
00440     curSoln_ = Teuchos::null;
00441     //
00442   } else {
00443     //
00444     // We are in the same rhs block, return if we are converged
00445     //
00446     if (status_==Passed) { return status_; }
00447   }
00448   //
00449   // get the native residual norms from the solver for this block of right-hand sides.
00450   // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00451   // Otherwise the native residual is assumed to be stored in the resvector_.
00452   //
00453   std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00454   RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00455   if ( residMV != Teuchos::null ) { 
00456     tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00457     MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00458     typename std::vector<int>::iterator p = curLSIdx_.begin();
00459     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00460       // Check if this index is valid
00461       if (*p != -1)  
00462   resvector_[*p] = tmp_resvector[i]; 
00463     }
00464   } else {
00465     typename std::vector<int>::iterator p = curLSIdx_.begin();
00466     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00467       // Check if this index is valid
00468       if (*p != -1)
00469   resvector_[*p] = tmp_resvector[i];
00470     }
00471   }
00472   //
00473   // Compute the new linear system residuals for testing.
00474   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00475   //
00476   if ( scalevector_.size() > 0 ) {
00477     typename std::vector<int>::iterator p = curLSIdx_.begin();
00478     for (; p<curLSIdx_.end(); ++p) {
00479       // Check if this index is valid
00480       if (*p != -1) {     
00481         // Scale the vector accordingly
00482         if ( scalevector_[ *p ] != zero ) {
00483           // Don't intentionally divide by zero.
00484           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00485         } else {
00486           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00487         }
00488       }
00489     }
00490   }
00491   else {
00492     typename std::vector<int>::iterator p = curLSIdx_.begin();
00493     for (; p<curLSIdx_.end(); ++p) {
00494       // Check if this index is valid
00495       if (*p != -1)     
00496         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00497     }
00498   } 
00499 
00500   // Check status of new linear system residuals and see if we have the quorum.
00501   int have = 0;
00502   ind_.resize( curLSIdx_.size() );
00503   std::vector<int> lclInd( curLSIdx_.size() );
00504   typename std::vector<int>::iterator p = curLSIdx_.begin();
00505   for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00506     // Check if this index is valid
00507     if (*p != -1) {     
00508       // Check if any of the residuals are larger than the tolerance.
00509       if (testvector_[ *p ] > currTolerance_) {
00510         // do nothing.
00511       } else if (testvector_[ *p ] <= currTolerance_) { 
00512         ind_[have] = *p;
00513         lclInd[have] = i;
00514         have++;
00515       } else {
00516         // Throw an std::exception if a NaN is found.
00517         status_ = Failed;
00518         TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestResNorm::checkStatus(): NaN has been detected.");
00519       }
00520     }
00521   } 
00522   ind_.resize(have);
00523   lclInd.resize(have);
00524   
00525   // Now check the exact residuals
00526   if (have) {
00527     RCP<MV> cur_update = iSolver->getCurrentUpdate();
00528     curSoln_ = lp.updateSolution( cur_update );
00529     RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_) );
00530     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00531     tmp_resvector.resize( MVT::GetNumberVecs( *cur_res ) );
00532     std::vector<MagnitudeType> tmp_testvector( have );
00533     MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
00534     
00535     if ( scalevector_.size() > 0 ) {
00536       for (int i=0; i<have; ++i) {
00537   // Scale the vector accordingly
00538   if ( scalevector_[ ind_[i] ] != zero ) {
00539     // Don't intentionally divide by zero.
00540     tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevector_[ ind_[i] ] / scalevalue_;
00541   } else {
00542     tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00543   }
00544       }
00545     }
00546     else {
00547       for (int i=0; i<have; ++i) {
00548   tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00549       }
00550     } 
00551 
00552     // Check if we want to keep the linear system and try to reduce the residual more.
00553     int have2 = 0;
00554     for (int i=0; i<have; ++i) {
00555       MagnitudeType diff = Teuchos::ScalarTraits<MagnitudeType>::magnitude( testvector_[ ind_[i] ]-tmp_testvector[ i ] );
00556       if (tmp_testvector[ i ] <= currTolerance_) {
00557         ind_[have2] = ind_[i];
00558         have2++;
00559       }
00560       else if (diff > currTolerance_) {
00561         lossDetected_ = true;
00562         ind_[have2] = ind_[i];
00563         have2++;
00564       } 
00565       else {
00566         currTolerance_ = currTolerance_ - 1.5*diff;
00567         while (currTolerance_ < 0.0) currTolerance_ += 0.1*diff;
00568       }  
00569 
00570     }
00571     have = have2;
00572     ind_.resize(have);
00573   }
00574 
00575   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00576   status_ = (have >= need) ? Passed : Failed;
00577   
00578   // Return the current status
00579   return status_;
00580 }
00581 
00582 template <class ScalarType, class MV, class OP>
00583 void StatusTestImpResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00584 {
00585   for (int j = 0; j < indent; j ++)
00586     os << ' ';
00587   printStatus(os, status_);
00588   os << resFormStr();
00589   if (status_==Undefined)
00590     os << ", tol = " << tolerance_ << std::endl;
00591   else {
00592     os << std::endl;
00593     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00594       const MagnitudeType maxRelRes = *std::max_element(
00595         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00596         );
00597       for (int j = 0; j < indent + 13; j ++)
00598         os << ' ';
00599       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00600          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00601     }
00602     else {
00603       for ( int i=0; i<numrhs_; i++ ) {
00604         for (int j = 0; j < indent + 13; j ++)
00605           os << ' ';
00606         os << "residual [ " << i << " ] = " << testvector_[ i ];
00607         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00608       }
00609     }
00610   }
00611   os << std::endl;
00612 }
00613 
00614 template <class ScalarType, class MV, class OP>
00615 void StatusTestImpResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00616 {
00617   os << std::left << std::setw(13) << std::setfill('.');
00618   switch (type) {
00619   case  Passed:
00620     os << "Converged";
00621     break;
00622   case  Failed:
00623     if (lossDetected_)
00624       os << "Unconverged (LoA)";
00625     else
00626       os << "Unconverged";
00627     break;
00628   case  Undefined:
00629   default:
00630     os << "**";
00631     break;
00632   }
00633   os << std::left << std::setfill(' ');
00634     return;
00635 }
00636 
00637 template <class ScalarType, class MV, class OP>
00638 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00639 {
00640   int i;
00641   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00642   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00643   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00644   // Compute scaling term (done once for each block that's being solved)
00645   if (firstcallCheckStatus_) {
00646     //
00647     // Get some current solver information.
00648     //
00649     firstcallCheckStatus_ = false;
00650 
00651     if (scaletype_== NormOfRHS) {
00652       RCP<const MV> rhs = lp.getRHS();
00653       numrhs_ = MVT::GetNumberVecs( *rhs );
00654       scalevector_.resize( numrhs_ );
00655       resvector_.resize( numrhs_ ); 
00656       testvector_.resize( numrhs_ );
00657       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00658     }
00659     else if (scaletype_==NormOfInitRes) {
00660       RCP<const MV> init_res = lp.getInitResVec();
00661       numrhs_ = MVT::GetNumberVecs( *init_res );
00662       scalevector_.resize( numrhs_ );
00663       resvector_.resize( numrhs_ ); 
00664       testvector_.resize( numrhs_ );
00665       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00666     }
00667     else if (scaletype_==NormOfPrecInitRes) {
00668       RCP<const MV> init_res = lp.getInitPrecResVec();
00669       numrhs_ = MVT::GetNumberVecs( *init_res );
00670       scalevector_.resize( numrhs_ );
00671       resvector_.resize( numrhs_ ); 
00672       testvector_.resize( numrhs_ );
00673       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00674     }
00675 
00676     curLSNum_ = lp.getLSNumber();
00677     curLSIdx_ = lp.getLSIndex();
00678     curBlksz_ = (int)curLSIdx_.size();
00679     int validLS = 0;
00680     for (i=0; i<curBlksz_; ++i) {
00681       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00682         validLS++;
00683     }
00684     curNumRHS_ = validLS;
00685     //
00686     // Initialize the testvector.
00687     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00688 
00689     // Return an error if the scaling is zero.
00690     if (scalevalue_ == zero) {
00691       return Failed;
00692     }
00693   }
00694   return Undefined;
00695 }
00696 
00697 } // end namespace Belos
00698 
00699 #endif /* BELOS_STATUS_TEST_IMPRESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:15 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3