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

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7