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