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 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_GEN_RESNORM_H
00030 #define BELOS_STATUS_TEST_GEN_RESNORM_H
00031 
00037 #include "BelosStatusTestResNorm.hpp"
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosMultiVecTraits.hpp"
00040 
00063 namespace Belos {
00064 
00065 template <class ScalarType, class MV, class OP>
00066 class StatusTestGenResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00067 
00068  public:
00069 
00070   // Convenience typedefs
00071   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00072   typedef typename SCT::magnitudeType MagnitudeType;
00073   typedef MultiVecTraits<ScalarType,MV>  MVT;
00074 
00076 
00077 
00081   enum ResType {Implicit, 
00082     Explicit  
00084   };
00085 
00087 
00089 
00090 
00091 
00104   StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
00105 
00107   virtual ~StatusTestGenResNorm();
00109 
00111 
00112 
00114 
00121   int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
00122   
00124 
00144   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00145 
00147 
00150   int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
00151 
00154   int setQuorum(int quorum) {quorum_ = quorum; return(0);}
00155 
00157   int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
00158 
00160 
00162 
00163 
00164 
00170   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00171 
00173   StatusType getStatus() const {return(status_);};
00175 
00177 
00178  
00180   void reset();
00181 
00183 
00185 
00186 
00188   void print(std::ostream& os, int indent = 0) const;
00189 
00191   void printStatus(std::ostream& os, StatusType type) const; 
00193 
00195 
00196 
00200   Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
00201 
00204   int getQuorum() const { return quorum_; }
00205 
00207   std::vector<int> convIndices() { return ind_; }
00208 
00210   MagnitudeType getTolerance() const {return(tolerance_);};
00211   
00213   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00214 
00216   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00217 
00219   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00220 
00223   bool getLOADetected() const { return false; }
00224 
00226 
00227 
00230 
00236   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00237 
00238  protected:
00239 
00240  private:
00241 
00243 
00244   
00246   MagnitudeType tolerance_;
00247 
00249   int quorum_;
00250 
00252   bool showMaxResNormOnly_;
00253  
00255   ResType restype_;
00256   
00258   NormType resnormtype_;
00259   
00261   ScaleType scaletype_;
00262   
00264   NormType scalenormtype_;
00265 
00267   MagnitudeType scalevalue_;
00268 
00270   std::vector<MagnitudeType> scalevector_;
00271   
00273   std::vector<MagnitudeType> resvector_;
00274 
00276   std::vector<MagnitudeType> testvector_;
00277 
00279   std::vector<int> ind_;
00280 
00282   Teuchos::RCP<MV> curSoln_;
00283   
00285   StatusType status_;
00286   
00288   int curBlksz_;
00289 
00291   int curNumRHS_;
00292 
00294   std::vector<int> curLSIdx_;
00295 
00297   int curLSNum_;
00298 
00300   int numrhs_;
00301 
00303   bool firstcallCheckStatus_;
00304 
00306   bool firstcallDefineResForm_;
00307 
00309   bool firstcallDefineScaleForm_;
00310 
00312 
00313 };
00314 
00315 template <class ScalarType, class MV, class OP>
00316 StatusTestGenResNorm<ScalarType,MV,OP>::StatusTestGenResNorm( MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly )
00317   : tolerance_(Tolerance),
00318     quorum_(quorum),
00319     showMaxResNormOnly_(showMaxResNormOnly),
00320     restype_(Implicit),
00321     resnormtype_(TwoNorm),  
00322     scaletype_(NormOfInitRes),
00323     scalenormtype_(TwoNorm),
00324     scalevalue_(1.0),
00325     status_(Undefined),
00326     curBlksz_(0),
00327     curLSNum_(0),
00328     numrhs_(0),
00329     firstcallCheckStatus_(true),
00330     firstcallDefineResForm_(true),
00331     firstcallDefineScaleForm_(true)
00332 {
00333   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00334   // the implicit residual std::vector.
00335 }
00336 
00337 template <class ScalarType, class MV, class OP>
00338 StatusTestGenResNorm<ScalarType,MV,OP>::~StatusTestGenResNorm() 
00339 {}
00340 
00341 template <class ScalarType, class MV, class OP>
00342 void StatusTestGenResNorm<ScalarType,MV,OP>::reset() 
00343 {
00344   status_ = Undefined;
00345   curBlksz_ = 0;
00346   curLSNum_ = 0;
00347   curLSIdx_.resize(0);
00348   numrhs_ = 0;
00349   ind_.resize(0);
00350   firstcallCheckStatus_ = true;
00351   curSoln_ = Teuchos::null;
00352 }
00353 
00354 template <class ScalarType, class MV, class OP>
00355 int StatusTestGenResNorm<ScalarType,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
00356 {    
00357   TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00358   "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
00359   firstcallDefineResForm_ = false;
00360     
00361   restype_ = TypeOfResidual;
00362   resnormtype_ = TypeOfNorm;
00363     
00364   return(0);
00365 }
00366 
00367 template <class ScalarType, class MV, class OP> 
00368 int StatusTestGenResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00369                                                          MagnitudeType ScaleValue )
00370 {
00371   TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00372   "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
00373   firstcallDefineScaleForm_ = false;
00374     
00375   scaletype_ = TypeOfScaling;
00376   scalenormtype_ = TypeOfNorm;
00377   scalevalue_ = ScaleValue;
00378     
00379   return(0);
00380 }
00381 
00382 template <class ScalarType, class MV, class OP>
00383 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
00384 {
00385   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00386   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00387   // Compute scaling term (done once for each block that's being solved)
00388   if (firstcallCheckStatus_) {
00389     StatusType status = firstCallCheckStatusSetup(iSolver);
00390     if(status==Failed) {
00391       status_ = Failed;
00392       return(status_);
00393     }
00394   }
00395   //
00396   // This section computes the norm of the residual std::vector
00397   //
00398   if ( curLSNum_ != lp.getLSNumber() ) {
00399     //
00400     // We have moved on to the next rhs block
00401     //
00402     curLSNum_ = lp.getLSNumber();
00403     curLSIdx_ = lp.getLSIndex();
00404     curBlksz_ = (int)curLSIdx_.size();
00405     int validLS = 0;
00406     for (int i=0; i<curBlksz_; ++i) {
00407       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_) 
00408   validLS++;
00409     }
00410     curNumRHS_ = validLS; 
00411     curSoln_ = Teuchos::null;
00412     //
00413   } else {
00414     //
00415     // We are in the same rhs block, return if we are converged
00416     //
00417     if (status_==Passed) { return status_; }
00418   }
00419   if (restype_==Implicit) {
00420     //
00421     // get the native residual norms from the solver for this block of right-hand sides.
00422     // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00423     // Otherwise the native residual is assumed to be stored in the resvector_.
00424     //
00425     std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00426     RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );     
00427     if ( residMV != Teuchos::null ) { 
00428       tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00429       MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );    
00430       typename std::vector<int>::iterator p = curLSIdx_.begin();
00431       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00432         // Check if this index is valid
00433         if (*p != -1)  
00434           resvector_[*p] = tmp_resvector[i]; 
00435       }
00436     } else {
00437       typename std::vector<int>::iterator p = curLSIdx_.begin();
00438       for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00439         // Check if this index is valid
00440         if (*p != -1)
00441           resvector_[*p] = tmp_resvector[i];
00442       }
00443     }
00444   }
00445   else if (restype_==Explicit) {
00446     //
00447     // Request the true residual for this block of right-hand sides.
00448     //
00449     RCP<MV> cur_update = iSolver->getCurrentUpdate();
00450     curSoln_ = lp.updateSolution( cur_update );
00451     RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
00452     lp.computeCurrResVec( &*cur_res, &*curSoln_ );
00453     std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
00454     MVT::MvNorm( *cur_res, 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   }
00462   //
00463   // Compute the new linear system residuals for testing.
00464   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00465   //
00466   if ( scalevector_.size() > 0 ) {
00467     typename std::vector<int>::iterator p = curLSIdx_.begin();
00468     for (; p<curLSIdx_.end(); ++p) {
00469       // Check if this index is valid
00470       if (*p != -1) {     
00471         // Scale the std::vector accordingly
00472         if ( scalevector_[ *p ] != zero ) {
00473           // Don't intentionally divide by zero.
00474           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00475         } else {
00476           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00477         }
00478       }
00479     }
00480   }
00481   else {
00482     typename std::vector<int>::iterator p = curLSIdx_.begin();
00483     for (; p<curLSIdx_.end(); ++p) {
00484       // Check if this index is valid
00485       if (*p != -1)     
00486         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00487     }
00488   } 
00489 
00490   // Check status of new linear system residuals and see if we have the quorum.
00491   int have = 0;
00492   ind_.resize( curLSIdx_.size() );
00493   typename std::vector<int>::iterator p = curLSIdx_.begin();
00494   for (; p<curLSIdx_.end(); ++p) {
00495     // Check if this index is valid
00496     if (*p != -1) {     
00497       // Check if any of the residuals are larger than the tolerance.
00498       if (testvector_[ *p ] > tolerance_) {
00499         // do nothing.
00500       } else if (testvector_[ *p ] <= tolerance_) { 
00501         ind_[have] = *p;
00502         have++;
00503       } else {
00504         // Throw an std::exception if a NaN is found.
00505         status_ = Failed;
00506         TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
00507       }
00508     }
00509   } 
00510   ind_.resize(have);
00511   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00512   status_ = (have >= need) ? Passed : Failed;
00513   
00514   // Return the current status
00515   return status_;
00516 }
00517 
00518 template <class ScalarType, class MV, class OP>
00519 void StatusTestGenResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00520 {
00521   for (int j = 0; j < indent; j ++)
00522     os << ' ';
00523   printStatus(os, status_);
00524   os << "(";
00525   os << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00526   os << ((restype_==Explicit) ? " Exp" : " Imp");
00527   os << " Res Vec) ";
00528   if (scaletype_!=None)
00529     os << "/ ";
00530   if (scaletype_==UserProvided)
00531     os << " (User Scale)";
00532   else {
00533     os << "(";
00534     os << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00535     if (scaletype_==NormOfInitRes)
00536       os << " Res0";
00537     else if (scaletype_==NormOfPrecInitRes)
00538       os << " Prec Res0";
00539     else
00540       os << " RHS ";
00541     os << ")";
00542   }
00543   if (status_==Undefined)
00544     os << ", tol = " << tolerance_ << std::endl;
00545   else {
00546     os << std::endl;
00547     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00548       const MagnitudeType maxRelRes = *std::max_element(
00549         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00550         );
00551       for (int j = 0; j < indent + 13; j ++)
00552         os << ' ';
00553       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00554          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00555     }
00556     else {
00557       for ( int i=0; i<numrhs_; i++ ) {
00558         for (int j = 0; j < indent + 13; j ++)
00559           os << ' ';
00560         os << "residual [ " << i << " ] = " << testvector_[ i ];
00561         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00562       }
00563     }
00564   }
00565   os << std::endl;
00566 }
00567 
00568 template <class ScalarType, class MV, class OP>
00569 void StatusTestGenResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const 
00570 {
00571   os << std::left << std::setw(13) << std::setfill('.');
00572   switch (type) {
00573   case  Passed:
00574     os << "Converged";
00575     break;
00576   case  Failed:
00577     os << "Unconverged";
00578     break;
00579   case  Undefined:
00580   default:
00581     os << "**";
00582     break;
00583   }
00584   os << std::left << std::setfill(' ');
00585     return;
00586 }
00587 
00588 template <class ScalarType, class MV, class OP>
00589 StatusType StatusTestGenResNorm<ScalarType,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
00590 {
00591   int i;
00592   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00593   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00594   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00595   // Compute scaling term (done once for each block that's being solved)
00596   if (firstcallCheckStatus_) {
00597     //
00598     // Get some current solver information.
00599     //
00600     firstcallCheckStatus_ = false;
00601 
00602     if (scaletype_== NormOfRHS) {
00603       RCP<const MV> rhs = lp.getRHS();
00604       numrhs_ = MVT::GetNumberVecs( *rhs );
00605       scalevector_.resize( numrhs_ );
00606       resvector_.resize( numrhs_ ); 
00607       testvector_.resize( numrhs_ );
00608       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00609     }
00610     else if (scaletype_==NormOfInitRes) {
00611       RCP<const MV> init_res = lp.getInitResVec();
00612       numrhs_ = MVT::GetNumberVecs( *init_res );
00613       scalevector_.resize( numrhs_ );
00614       resvector_.resize( numrhs_ ); 
00615       testvector_.resize( numrhs_ );
00616       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00617     }
00618     else if (scaletype_==NormOfPrecInitRes) {
00619       RCP<const MV> init_res = lp.getInitPrecResVec();
00620       numrhs_ = MVT::GetNumberVecs( *init_res );
00621       scalevector_.resize( numrhs_ );
00622       resvector_.resize( numrhs_ ); 
00623       testvector_.resize( numrhs_ );
00624       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00625     }
00626 
00627     curLSNum_ = lp.getLSNumber();
00628     curLSIdx_ = lp.getLSIndex();
00629     curBlksz_ = (int)curLSIdx_.size();
00630     int validLS = 0;
00631     for (i=0; i<curBlksz_; ++i) {
00632       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00633         validLS++;
00634     }
00635     curNumRHS_ = validLS;
00636     //
00637     // Initialize the testvector.
00638     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00639 
00640     // Return an error if the scaling is zero.
00641     if (scalevalue_ == zero) {
00642       return Failed;
00643     }
00644   }
00645   return Undefined;
00646 }
00647 
00648 } // end namespace Belos
00649 
00650 #endif /* BELOS_STATUS_TEST_RESNORM_H */

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7