Belos Version of the Day
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_STATUS_TEST_IMPRESNORM_H
00043 #define BELOS_STATUS_TEST_IMPRESNORM_H
00044 
00050 #include "BelosStatusTestResNorm.hpp"
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosMultiVecTraits.hpp"
00053 
00075 namespace Belos {
00076 
00077 template <class ScalarType, class MV, class OP>
00078 class StatusTestImpResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00079 
00080  public:
00081 
00082   // Convenience typedefs
00083   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00084   typedef typename SCT::magnitudeType MagnitudeType;
00085   typedef MultiVecTraits<ScalarType,MV>  MVT;
00086 
00088 
00089 
00090 
00103   StatusTestImpResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
00104 
00106   virtual ~StatusTestImpResNorm();
00108 
00110 
00111 
00113 
00119   int defineResForm( NormType TypeOfNorm);
00120   
00122 
00142   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00143 
00145 
00148   int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
00149 
00152   int setQuorum(int quorum) {quorum_ = quorum; return(0);}
00153 
00155   int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
00156 
00158 
00160 
00161 
00162 
00168   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00169 
00171   StatusType getStatus() const {return(status_);};
00173 
00175 
00176  
00178   void reset();
00179 
00181 
00183 
00184 
00186   void print(std::ostream& os, int indent = 0) const;
00187 
00189   void printStatus(std::ostream& os, StatusType type) const; 
00191 
00193 
00194 
00196   Teuchos::RCP<MV> getSolution() { return curSoln_; }
00197 
00200   int getQuorum() const { return quorum_; }
00201 
00203   bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
00204 
00206   std::vector<int> convIndices() { return ind_; }
00207   
00209   MagnitudeType getTolerance() const {return(tolerance_);};
00210   
00212   MagnitudeType getCurrTolerance() const {return(currTolerance_);};
00213   
00215   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00216 
00218   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00219 
00221   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00222 
00224   bool getLOADetected() const { return lossDetected_; }
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::StatusTestImpResNorm<>: " << 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 << " Res Vec) ";
00265 
00266     // If there is no residual scaling, return current string.
00267     if (scaletype_!=None)
00268     {
00269       // Insert division sign.
00270       oss << "/ ";
00271       
00272       // Determine output string for scaling, if there is any.
00273       if (scaletype_==UserProvided)
00274         oss << " (User Scale)";
00275       else {
00276         oss << "(";
00277         oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00278         if (scaletype_==NormOfInitRes)
00279           oss << " Res0";
00280         else if (scaletype_==NormOfPrecInitRes)
00281           oss << " Prec Res0";
00282         else
00283           oss << " RHS ";
00284         oss << ")";
00285       }
00286     } 
00287 
00288     return oss.str();
00289   }
00290 
00292 
00293   
00295   MagnitudeType tolerance_, currTolerance_;
00296 
00298   int quorum_;
00299 
00301   bool showMaxResNormOnly_;
00302  
00304   NormType resnormtype_;
00305   
00307   ScaleType scaletype_;
00308   
00310   NormType scalenormtype_;
00311 
00313   MagnitudeType scalevalue_;
00314 
00316   std::vector<MagnitudeType> scalevector_;
00317   
00319   std::vector<MagnitudeType> resvector_;
00320 
00322   std::vector<MagnitudeType> testvector_;
00323 
00325   Teuchos::RCP<MV> curSoln_; 
00326 
00328   std::vector<int> ind_;
00329   
00331   StatusType status_;
00332   
00334   int curBlksz_;
00335 
00337   int curNumRHS_;
00338 
00340   std::vector<int> curLSIdx_;
00341 
00343   int curLSNum_;
00344 
00346   int numrhs_;
00347 
00349   bool firstcallCheckStatus_;
00350 
00352   bool firstcallDefineResForm_;
00353 
00355   bool firstcallDefineScaleForm_;
00356 
00358   bool lossDetected_;
00360 
00361 };
00362 
00363 template <class ScalarType, class MV, class OP>
00364 StatusTestImpResNorm<ScalarType,MV,OP>::StatusTestImpResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00365   : tolerance_(Tolerance),
00366     currTolerance_(Tolerance),
00367     quorum_(quorum),
00368     showMaxResNormOnly_(showMaxResNormOnly),
00369     resnormtype_(TwoNorm),  
00370     scaletype_(NormOfInitRes),
00371     scalenormtype_(TwoNorm),
00372     scalevalue_(1.0),
00373     status_(Undefined),
00374     curBlksz_(0),
00375     curLSNum_(0),
00376     numrhs_(0),
00377     firstcallCheckStatus_(true),
00378     firstcallDefineResForm_(true),
00379     firstcallDefineScaleForm_(true),
00380     lossDetected_(false)
00381 {
00382   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00383   // the implicit residual vector.
00384 }
00385 
00386 template <class ScalarType, class MV, class OP>
00387 StatusTestImpResNorm<ScalarType,MV,OP>::~StatusTestImpResNorm() 
00388 {}
00389 
00390 template <class ScalarType, class MV, class OP>
00391 void StatusTestImpResNorm<ScalarType,MV,OP>::reset() 
00392 {
00393   status_ = Undefined;
00394   curBlksz_ = 0;
00395   curLSNum_ = 0;
00396   curLSIdx_.resize(0);
00397   numrhs_ = 0;
00398   ind_.resize(0);
00399   currTolerance_ = tolerance_;
00400   firstcallCheckStatus_ = true;
00401   lossDetected_ = false;
00402   curSoln_ = Teuchos::null;
00403 }
00404 
00405 template <class ScalarType, class MV, class OP>
00406 int StatusTestImpResNorm<ScalarType,MV,OP>::defineResForm( NormType TypeOfNorm )
00407 {    
00408   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00409   "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
00410   firstcallDefineResForm_ = false;
00411     
00412   resnormtype_ = TypeOfNorm;
00413     
00414   return(0);
00415 }
00416 
00417 template <class ScalarType, class MV, class OP> 
00418 int StatusTestImpResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00419                                                          MagnitudeType ScaleValue )
00420 {
00421   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00422   "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
00423   firstcallDefineScaleForm_ = false;
00424     
00425   scaletype_ = TypeOfScaling;
00426   scalenormtype_ = TypeOfNorm;
00427   scalevalue_ = ScaleValue;
00428     
00429   return(0);
00430 }
00431 
00432 template <class ScalarType, class MV, class OP>
00433 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00434 {
00435   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00436   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00437   // Compute scaling term (done once for each block that's being solved)
00438   if (firstcallCheckStatus_) {
00439     StatusType status = firstCallCheckStatusSetup(iSolver);
00440     if(status==Failed) {
00441       status_ = Failed;
00442       return(status_);
00443     }
00444   }
00445   //
00446   // This section computes the norm of the residual vector
00447   //
00448   if ( curLSNum_ != lp.getLSNumber() ) {
00449     //
00450     // We have moved on to the next rhs block
00451     //
00452     curLSNum_ = lp.getLSNumber();
00453     curLSIdx_ = lp.getLSIndex();
00454     curBlksz_ = (int)curLSIdx_.size();
00455     int validLS = 0;
00456     for (int i=0; i<curBlksz_; ++i) {
00457       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00458   validLS++;
00459     }
00460     curNumRHS_ = validLS; 
00461     curSoln_ = Teuchos::null;
00462     //
00463   } else {
00464     //
00465     // We are in the same rhs block, return if we are converged
00466     //
00467     if (status_==Passed) { return status_; }
00468   }
00469   //
00470   // get the native residual norms from the solver for this block of right-hand sides.
00471   // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00472   // Otherwise the native residual is assumed to be stored in the resvector_.
00473   //
00474   std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00475   Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00476   if ( residMV != Teuchos::null ) { 
00477     tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00478     MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00479     typename std::vector<int>::iterator p = curLSIdx_.begin();
00480     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00481       // Check if this index is valid
00482       if (*p != -1)  
00483   resvector_[*p] = tmp_resvector[i]; 
00484     }
00485   } else {
00486     typename std::vector<int>::iterator p = curLSIdx_.begin();
00487     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00488       // Check if this index is valid
00489       if (*p != -1)
00490   resvector_[*p] = tmp_resvector[i];
00491     }
00492   }
00493   //
00494   // Compute the new linear system residuals for testing.
00495   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00496   //
00497   if ( scalevector_.size() > 0 ) {
00498     typename std::vector<int>::iterator p = curLSIdx_.begin();
00499     for (; p<curLSIdx_.end(); ++p) {
00500       // Check if this index is valid
00501       if (*p != -1) {     
00502         // Scale the vector accordingly
00503         if ( scalevector_[ *p ] != zero ) {
00504           // Don't intentionally divide by zero.
00505           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00506         } else {
00507           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00508         }
00509       }
00510     }
00511   }
00512   else {
00513     typename std::vector<int>::iterator p = curLSIdx_.begin();
00514     for (; p<curLSIdx_.end(); ++p) {
00515       // Check if this index is valid
00516       if (*p != -1)     
00517         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00518     }
00519   } 
00520 
00521   // Check status of new linear system residuals and see if we have the quorum.
00522   int have = 0;
00523   ind_.resize( curLSIdx_.size() );
00524   std::vector<int> lclInd( curLSIdx_.size() );
00525   typename std::vector<int>::iterator p = curLSIdx_.begin();
00526   for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00527     // Check if this index is valid
00528     if (*p != -1) {     
00529       // Check if any of the residuals are larger than the tolerance.
00530       if (testvector_[ *p ] > currTolerance_) {
00531         // do nothing.
00532       } else if (testvector_[ *p ] <= currTolerance_) { 
00533         ind_[have] = *p;
00534         lclInd[have] = i;
00535         have++;
00536       } else {
00537         // Throw an std::exception if a NaN is found.
00538         status_ = Failed;
00539         TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestResNorm::checkStatus(): NaN has been detected.");
00540       }
00541     }
00542   } 
00543   ind_.resize(have);
00544   lclInd.resize(have);
00545   
00546   // Now check the exact residuals
00547   if (have) {
00548     Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
00549     curSoln_ = lp.updateSolution( cur_update );
00550     Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_) );
00551     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00552     tmp_resvector.resize( MVT::GetNumberVecs( *cur_res ) );
00553     std::vector<MagnitudeType> tmp_testvector( have );
00554     MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
00555     
00556     if ( scalevector_.size() > 0 ) {
00557       for (int i=0; i<have; ++i) {
00558   // Scale the vector accordingly
00559   if ( scalevector_[ ind_[i] ] != zero ) {
00560     // Don't intentionally divide by zero.
00561     tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevector_[ ind_[i] ] / scalevalue_;
00562   } else {
00563     tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00564   }
00565       }
00566     }
00567     else {
00568       for (int i=0; i<have; ++i) {
00569   tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00570       }
00571     } 
00572 
00573     // Check if we want to keep the linear system and try to reduce the residual more.
00574     int have2 = 0;
00575     for (int i=0; i<have; ++i) {
00576       MagnitudeType diff = Teuchos::ScalarTraits<MagnitudeType>::magnitude( testvector_[ ind_[i] ]-tmp_testvector[ i ] );
00577       if (tmp_testvector[ i ] <= currTolerance_) {
00578         ind_[have2] = ind_[i];
00579         have2++;
00580       }
00581       else if (diff > currTolerance_) {
00582         lossDetected_ = true;
00583         ind_[have2] = ind_[i];
00584         have2++;
00585       } 
00586       else {
00587         currTolerance_ = currTolerance_ - 1.5*diff;
00588         while (currTolerance_ < 0.0) currTolerance_ += 0.1*diff;
00589       }  
00590 
00591     }
00592     have = have2;
00593     ind_.resize(have);
00594   }
00595 
00596   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00597   status_ = (have >= need) ? Passed : Failed;
00598   
00599   // Return the current status
00600   return status_;
00601 }
00602 
00603 template <class ScalarType, class MV, class OP>
00604 void StatusTestImpResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00605 {
00606   for (int j = 0; j < indent; j ++)
00607     os << ' ';
00608   printStatus(os, status_);
00609   os << resFormStr();
00610   if (status_==Undefined)
00611     os << ", tol = " << tolerance_ << std::endl;
00612   else {
00613     os << std::endl;
00614     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00615       const MagnitudeType maxRelRes = *std::max_element(
00616         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00617         );
00618       for (int j = 0; j < indent + 13; j ++)
00619         os << ' ';
00620       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00621          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00622     }
00623     else {
00624       for ( int i=0; i<numrhs_; i++ ) {
00625         for (int j = 0; j < indent + 13; j ++)
00626           os << ' ';
00627         os << "residual [ " << i << " ] = " << testvector_[ i ];
00628         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00629       }
00630     }
00631   }
00632   os << std::endl;
00633 }
00634 
00635 template <class ScalarType, class MV, class OP>
00636 void StatusTestImpResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00637 {
00638   os << std::left << std::setw(13) << std::setfill('.');
00639   switch (type) {
00640   case  Passed:
00641     os << "Converged";
00642     break;
00643   case  Failed:
00644     if (lossDetected_)
00645       os << "Unconverged (LoA)";
00646     else
00647       os << "Unconverged";
00648     break;
00649   case  Undefined:
00650   default:
00651     os << "**";
00652     break;
00653   }
00654   os << std::left << std::setfill(' ');
00655     return;
00656 }
00657 
00658 template <class ScalarType, class MV, class OP>
00659 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00660 {
00661   int i;
00662   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00663   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00664   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00665   // Compute scaling term (done once for each block that's being solved)
00666   if (firstcallCheckStatus_) {
00667     //
00668     // Get some current solver information.
00669     //
00670     firstcallCheckStatus_ = false;
00671 
00672     if (scaletype_== NormOfRHS) {
00673       Teuchos::RCP<const MV> rhs = lp.getRHS();
00674       numrhs_ = MVT::GetNumberVecs( *rhs );
00675       scalevector_.resize( numrhs_ );
00676       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00677     }
00678     else if (scaletype_==NormOfInitRes) {
00679       Teuchos::RCP<const MV> init_res = lp.getInitResVec();
00680       numrhs_ = MVT::GetNumberVecs( *init_res );
00681       scalevector_.resize( numrhs_ );
00682       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00683     }
00684     else if (scaletype_==NormOfPrecInitRes) {
00685       Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
00686       numrhs_ = MVT::GetNumberVecs( *init_res );
00687       scalevector_.resize( numrhs_ );
00688       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00689     }
00690     else {
00691       numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
00692     }
00693 
00694     resvector_.resize( numrhs_ ); 
00695     testvector_.resize( numrhs_ );
00696 
00697     curLSNum_ = lp.getLSNumber();
00698     curLSIdx_ = lp.getLSIndex();
00699     curBlksz_ = (int)curLSIdx_.size();
00700     int validLS = 0;
00701     for (i=0; i<curBlksz_; ++i) {
00702       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00703         validLS++;
00704     }
00705     curNumRHS_ = validLS;
00706     //
00707     // Initialize the testvector.
00708     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00709 
00710     // Return an error if the scaling is zero.
00711     if (scalevalue_ == zero) {
00712       return Failed;
00713     }
00714   }
00715   return Undefined;
00716 }
00717 
00718 } // end namespace Belos
00719 
00720 #endif /* BELOS_STATUS_TEST_IMPRESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines