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