Belos Package Browser (Single Doxygen Collection) Development
BelosStatusTestGenResNorm.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_GEN_RESNORM_H
00043 #define BELOS_STATUS_TEST_GEN_RESNORM_H
00044 
00050 #include "BelosStatusTestResNorm.hpp"
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosMultiVecTraits.hpp"
00053 
00076 namespace Belos {
00077 
00078 template <class ScalarType, class MV, class OP>
00079 class StatusTestGenResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00080 
00081  public:
00082 
00083   // Convenience typedefs
00084   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00085   typedef typename SCT::magnitudeType MagnitudeType;
00086   typedef MultiVecTraits<ScalarType,MV>  MVT;
00087 
00089 
00090 
00094   enum ResType {Implicit, 
00095     Explicit  
00097   };
00098 
00100 
00102 
00103 
00104 
00117   StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
00118 
00120   virtual ~StatusTestGenResNorm();
00122 
00124 
00125 
00127 
00134   int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
00135   
00137 
00157   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00158 
00160 
00163   int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
00164 
00167   int setQuorum(int quorum) {quorum_ = quorum; return(0);}
00168 
00170   int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
00171 
00173 
00175 
00176 
00177 
00183   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00184 
00186   StatusType getStatus() const {return(status_);};
00188 
00190 
00191  
00193   void reset();
00194 
00196 
00198 
00199 
00201   void print(std::ostream& os, int indent = 0) const;
00202 
00204   void printStatus(std::ostream& os, StatusType type) const; 
00206 
00208 
00209 
00213   Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
00214 
00217   int getQuorum() const { return quorum_; }
00218 
00220   bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
00221 
00223   std::vector<int> convIndices() { return ind_; }
00224 
00226   MagnitudeType getTolerance() const {return(tolerance_);};
00227   
00229   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00230 
00232   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00233 
00235   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00236 
00239   bool getLOADetected() const { return false; }
00240 
00242 
00243 
00246 
00252   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00254 
00257 
00259   std::string description() const
00260   {
00261     std::ostringstream oss;
00262     oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
00263     oss << ", tol = " << tolerance_;
00264     return oss.str();
00265   }
00267 
00268  protected:
00269 
00270  private:
00271 
00273 
00274 
00275   std::string resFormStr() const
00276   {
00277     std::ostringstream oss;
00278     oss << "(";
00279     oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00280     oss << ((restype_==Explicit) ? " Exp" : " Imp");
00281     oss << " Res Vec) ";
00282     if (scaletype_!=None)
00283       oss << "/ ";
00284     if (scaletype_==UserProvided)
00285       oss << " (User Scale)";
00286     else {
00287       oss << "(";
00288       oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00289       if (scaletype_==NormOfInitRes)
00290         oss << " Res0";
00291       else if (scaletype_==NormOfPrecInitRes)
00292         oss << " Prec Res0";
00293       else
00294         oss << " RHS ";
00295       oss << ")";
00296     }
00297     return oss.str();
00298   }
00299 
00301 
00303 
00304   
00306   MagnitudeType tolerance_;
00307 
00309   int quorum_;
00310 
00312   bool showMaxResNormOnly_;
00313  
00315   ResType restype_;
00316   
00318   NormType resnormtype_;
00319   
00321   ScaleType scaletype_;
00322   
00324   NormType scalenormtype_;
00325 
00327   MagnitudeType scalevalue_;
00328 
00330   std::vector<MagnitudeType> scalevector_;
00331   
00333   std::vector<MagnitudeType> resvector_;
00334 
00336   std::vector<MagnitudeType> testvector_;
00337 
00339   std::vector<int> ind_;
00340 
00342   Teuchos::RCP<MV> curSoln_;
00343   
00345   StatusType status_;
00346   
00348   int curBlksz_;
00349 
00351   int curNumRHS_;
00352 
00354   std::vector<int> curLSIdx_;
00355 
00357   int curLSNum_;
00358 
00360   int numrhs_;
00361 
00363   bool firstcallCheckStatus_;
00364 
00366   bool firstcallDefineResForm_;
00367 
00369   bool firstcallDefineScaleForm_;
00370 
00372 
00373 };
00374 
00375 template <class ScalarType, class MV, class OP>
00376 StatusTestGenResNorm<ScalarType,MV,OP>::StatusTestGenResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00377   : tolerance_(Tolerance),
00378     quorum_(quorum),
00379     showMaxResNormOnly_(showMaxResNormOnly),
00380     restype_(Implicit),
00381     resnormtype_(TwoNorm),  
00382     scaletype_(NormOfInitRes),
00383     scalenormtype_(TwoNorm),
00384     scalevalue_(1.0),
00385     status_(Undefined),
00386     curBlksz_(0),
00387     curLSNum_(0),
00388     numrhs_(0),
00389     firstcallCheckStatus_(true),
00390     firstcallDefineResForm_(true),
00391     firstcallDefineScaleForm_(true)
00392 {
00393   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00394   // the implicit residual std::vector.
00395 }
00396 
00397 template <class ScalarType, class MV, class OP>
00398 StatusTestGenResNorm<ScalarType,MV,OP>::~StatusTestGenResNorm() 
00399 {}
00400 
00401 template <class ScalarType, class MV, class OP>
00402 void StatusTestGenResNorm<ScalarType,MV,OP>::reset() 
00403 {
00404   status_ = Undefined;
00405   curBlksz_ = 0;
00406   curLSNum_ = 0;
00407   curLSIdx_.resize(0);
00408   numrhs_ = 0;
00409   ind_.resize(0);
00410   firstcallCheckStatus_ = true;
00411   curSoln_ = Teuchos::null;
00412 }
00413 
00414 template <class ScalarType, class MV, class OP>
00415 int StatusTestGenResNorm<ScalarType,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
00416 {    
00417   TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00418   "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
00419   firstcallDefineResForm_ = false;
00420     
00421   restype_ = TypeOfResidual;
00422   resnormtype_ = TypeOfNorm;
00423     
00424   return(0);
00425 }
00426 
00427 template <class ScalarType, class MV, class OP> 
00428 int StatusTestGenResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00429                                                          MagnitudeType ScaleValue )
00430 {
00431   TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00432   "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
00433   firstcallDefineScaleForm_ = false;
00434     
00435   scaletype_ = TypeOfScaling;
00436   scalenormtype_ = TypeOfNorm;
00437   scalevalue_ = ScaleValue;
00438     
00439   return(0);
00440 }
00441 
00442 template <class ScalarType, class MV, class OP>
00443 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00444 {
00445   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00446   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00447   // Compute scaling term (done once for each block that's being solved)
00448   if (firstcallCheckStatus_) {
00449     StatusType status = firstCallCheckStatusSetup(iSolver);
00450     if(status==Failed) {
00451       status_ = Failed;
00452       return(status_);
00453     }
00454   }
00455   //
00456   // This section computes the norm of the residual std::vector
00457   //
00458   if ( curLSNum_ != lp.getLSNumber() ) {
00459     //
00460     // We have moved on to the next rhs block
00461     //
00462     curLSNum_ = lp.getLSNumber();
00463     curLSIdx_ = lp.getLSIndex();
00464     curBlksz_ = (int)curLSIdx_.size();
00465     int validLS = 0;
00466     for (int i=0; i<curBlksz_; ++i) {
00467       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00468   validLS++;
00469     }
00470     curNumRHS_ = validLS; 
00471     curSoln_ = Teuchos::null;
00472     //
00473   } else {
00474     //
00475     // We are in the same rhs block, return if we are converged
00476     //
00477     if (status_==Passed) { return status_; }
00478   }
00479   if (restype_==Implicit) {
00480     //
00481     // get the native residual norms from the solver for this block of right-hand sides.
00482     // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00483     // Otherwise the native residual is assumed to be stored in the resvector_.
00484     //
00485     std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00486     RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00487     if ( residMV != Teuchos::null ) { 
00488       tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00489       MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00490       typename std::vector<int>::iterator p = curLSIdx_.begin();
00491       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00492         // Check if this index is valid
00493         if (*p != -1)  
00494           resvector_[*p] = tmp_resvector[i]; 
00495       }
00496     } else {
00497       typename std::vector<int>::iterator p = curLSIdx_.begin();
00498       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00499         // Check if this index is valid
00500         if (*p != -1)
00501           resvector_[*p] = tmp_resvector[i];
00502       }
00503     }
00504   }
00505   else if (restype_==Explicit) {
00506     //
00507     // Request the true residual for this block of right-hand sides.
00508     //
00509     RCP<MV> cur_update = iSolver->getCurrentUpdate();
00510     curSoln_ = lp.updateSolution( cur_update );
00511     RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
00512     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00513     std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
00514     MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
00515     typename std::vector<int>::iterator p = curLSIdx_.begin();
00516     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00517       // Check if this index is valid
00518       if (*p != -1)
00519         resvector_[*p] = tmp_resvector[i];
00520     }
00521   }
00522   //
00523   // Compute the new linear system residuals for testing.
00524   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00525   //
00526   if ( scalevector_.size() > 0 ) {
00527     typename std::vector<int>::iterator p = curLSIdx_.begin();
00528     for (; p<curLSIdx_.end(); ++p) {
00529       // Check if this index is valid
00530       if (*p != -1) {     
00531         // Scale the std::vector accordingly
00532         if ( scalevector_[ *p ] != zero ) {
00533           // Don't intentionally divide by zero.
00534           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00535         } else {
00536           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00537         }
00538       }
00539     }
00540   }
00541   else {
00542     typename std::vector<int>::iterator p = curLSIdx_.begin();
00543     for (; p<curLSIdx_.end(); ++p) {
00544       // Check if this index is valid
00545       if (*p != -1)     
00546         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00547     }
00548   } 
00549 
00550   // Check status of new linear system residuals and see if we have the quorum.
00551   int have = 0;
00552   ind_.resize( curLSIdx_.size() );
00553   typename std::vector<int>::iterator p = curLSIdx_.begin();
00554   for (; p<curLSIdx_.end(); ++p) {
00555     // Check if this index is valid
00556     if (*p != -1) {     
00557       // Check if any of the residuals are larger than the tolerance.
00558       if (testvector_[ *p ] > tolerance_) {
00559         // do nothing.
00560       } else if (testvector_[ *p ] <= tolerance_) { 
00561         ind_[have] = *p;
00562         have++;
00563       } else {
00564         // Throw an std::exception if a NaN is found.
00565         status_ = Failed;
00566         TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
00567       }
00568     }
00569   } 
00570   ind_.resize(have);
00571   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00572   status_ = (have >= need) ? Passed : Failed;
00573   
00574   // Return the current status
00575   return status_;
00576 }
00577 
00578 template <class ScalarType, class MV, class OP>
00579 void StatusTestGenResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00580 {
00581   for (int j = 0; j < indent; j ++)
00582     os << ' ';
00583   printStatus(os, status_);
00584   os << resFormStr();
00585   if (status_==Undefined)
00586     os << ", tol = " << tolerance_ << std::endl;
00587   else {
00588     os << std::endl;
00589     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00590       const MagnitudeType maxRelRes = *std::max_element(
00591         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00592         );
00593       for (int j = 0; j < indent + 13; j ++)
00594         os << ' ';
00595       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00596          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00597     }
00598     else {
00599       for ( int i=0; i<numrhs_; i++ ) {
00600         for (int j = 0; j < indent + 13; j ++)
00601           os << ' ';
00602         os << "residual [ " << i << " ] = " << testvector_[ i ];
00603         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00604       }
00605     }
00606   }
00607   os << std::endl;
00608 }
00609 
00610 template <class ScalarType, class MV, class OP>
00611 void StatusTestGenResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00612 {
00613   os << std::left << std::setw(13) << std::setfill('.');
00614   switch (type) {
00615   case  Passed:
00616     os << "Converged";
00617     break;
00618   case  Failed:
00619     os << "Unconverged";
00620     break;
00621   case  Undefined:
00622   default:
00623     os << "**";
00624     break;
00625   }
00626   os << std::left << std::setfill(' ');
00627     return;
00628 }
00629 
00630 template <class ScalarType, class MV, class OP>
00631 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00632 {
00633   int i;
00634   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00635   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00636   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00637   // Compute scaling term (done once for each block that's being solved)
00638   if (firstcallCheckStatus_) {
00639     //
00640     // Get some current solver information.
00641     //
00642     firstcallCheckStatus_ = false;
00643 
00644     if (scaletype_== NormOfRHS) {
00645       RCP<const MV> rhs = lp.getRHS();
00646       numrhs_ = MVT::GetNumberVecs( *rhs );
00647       scalevector_.resize( numrhs_ );
00648       resvector_.resize( numrhs_ ); 
00649       testvector_.resize( numrhs_ );
00650       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00651     }
00652     else if (scaletype_==NormOfInitRes) {
00653       RCP<const MV> init_res = lp.getInitResVec();
00654       numrhs_ = MVT::GetNumberVecs( *init_res );
00655       scalevector_.resize( numrhs_ );
00656       resvector_.resize( numrhs_ ); 
00657       testvector_.resize( numrhs_ );
00658       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00659     }
00660     else if (scaletype_==NormOfPrecInitRes) {
00661       RCP<const MV> init_res = lp.getInitPrecResVec();
00662       numrhs_ = MVT::GetNumberVecs( *init_res );
00663       scalevector_.resize( numrhs_ );
00664       resvector_.resize( numrhs_ ); 
00665       testvector_.resize( numrhs_ );
00666       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00667     }
00668 
00669     curLSNum_ = lp.getLSNumber();
00670     curLSIdx_ = lp.getLSIndex();
00671     curBlksz_ = (int)curLSIdx_.size();
00672     int validLS = 0;
00673     for (i=0; i<curBlksz_; ++i) {
00674       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00675         validLS++;
00676     }
00677     curNumRHS_ = validLS;
00678     //
00679     // Initialize the testvector.
00680     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00681 
00682     // Return an error if the scaling is zero.
00683     if (scalevalue_ == zero) {
00684       return Failed;
00685     }
00686   }
00687   return Undefined;
00688 }
00689 
00690 } // end namespace Belos
00691 
00692 #endif /* BELOS_STATUS_TEST_RESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines