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

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