BelosStatusTestResNorm.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_RESNORM_H
00030 #define BELOS_STATUS_TEST_RESNORM_H
00031 
00037 #include "BelosStatusTest.hpp"
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosMultiVecTraits.hpp"
00040 
00063 namespace Belos {
00064 
00065 template <class ScalarType, class MV, class OP>
00066 class StatusTestResNorm: public StatusTest<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 
00080   enum ResType {Implicit, 
00081     Explicit  
00083   };
00088   enum ScaleType {NormOfRHS,     
00089       NormOfInitRes, 
00090       NormOfPrecInitRes, 
00091       None,          
00092       UserProvided   
00094   };
00096 
00098 
00099 
00100 
00110   StatusTestResNorm( MagnitudeType Tolerance, bool showMaxResNormOnly = false );
00111 
00113   virtual ~StatusTestResNorm();
00115 
00117 
00118 
00120 
00127   int DefineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
00128   
00130 
00150   int DefineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00151 
00153 
00156   int ResetTolerance(MagnitudeType Tolerance) {tolerance_ = Tolerance; return(0);};
00157 
00159 
00161 
00162 
00163 
00169   StatusType CheckStatus(IterativeSolver<ScalarType,MV,OP>* iSolver);
00170 
00172   StatusType GetStatus() const {return(status_);};
00174 
00176 
00177  
00179   void Reset();
00180 
00182 
00184 
00185 
00187 
00195   bool ResidualVectorRequired() const { return(resvecrequired_); };
00196 
00198 
00200 
00201 
00203   ostream& Print(ostream& os, int indent = 0) const;
00205 
00207 
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 
00222 
00223 
00226 
00232   StatusType FirstcallCheckStatusSetup(IterativeSolver<ScalarType,MV,OP>* iSolver);
00233 
00234  protected:
00235 
00236  private:
00237 
00239 
00240   
00242   MagnitudeType tolerance_;
00243 
00245   bool showMaxResNormOnly_;
00246  
00248   ResType restype_;
00249   
00251   NormType resnormtype_;
00252   
00254   ScaleType scaletype_;
00255   
00257   NormType scalenormtype_;
00258 
00260   MagnitudeType scalevalue_;
00261 
00263   std::vector<MagnitudeType> scalevector_;
00264   
00266   std::vector<MagnitudeType> resvector_;
00267 
00269   std::vector<MagnitudeType> testvector_;
00270   
00272   StatusType status_;
00273   
00275   int cur_rhs_num_;
00276 
00278   int cur_blksz_;
00279 
00281   int cur_num_rhs_;
00282 
00284   int numrhs_;
00285 
00287   bool resvecrequired_;
00288 
00290   bool firstcallCheckStatus_;
00291 
00293   bool firstcallDefineResForm_;
00294 
00296   bool firstcallDefineScaleForm_;
00297 
00299 
00300 };
00301 
00302 template <class ScalarType, class MV, class OP>
00303 StatusTestResNorm<ScalarType,MV,OP>::StatusTestResNorm( MagnitudeType Tolerance, bool showMaxResNormOnly )
00304   : tolerance_(Tolerance),
00305     showMaxResNormOnly_(showMaxResNormOnly),
00306     restype_(Implicit),
00307     resnormtype_(TwoNorm),  
00308     scaletype_(NormOfInitRes),
00309     scalenormtype_(TwoNorm),
00310     scalevalue_(1.0),
00311     status_(Unchecked),
00312     cur_rhs_num_(0),
00313     cur_blksz_(0),
00314     numrhs_(0),
00315     resvecrequired_(false),
00316     firstcallCheckStatus_(true),
00317     firstcallDefineResForm_(true),
00318     firstcallDefineScaleForm_(true)
00319 {
00320   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00321   // the implicit residual vector.
00322 }
00323 
00324 template <class ScalarType, class MV, class OP>
00325 StatusTestResNorm<ScalarType,MV,OP>::~StatusTestResNorm() 
00326 {}
00327 
00328 template <class ScalarType, class MV, class OP>
00329 void StatusTestResNorm<ScalarType,MV,OP>::Reset() 
00330 {
00331   status_ = Unchecked;
00332   cur_rhs_num_ = 0;
00333   cur_blksz_ = 0;
00334   numrhs_ = 0;
00335   firstcallCheckStatus_ = true;
00336 }
00337 
00338 template <class ScalarType, class MV, class OP>
00339 int StatusTestResNorm<ScalarType,MV,OP>::DefineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
00340 {    
00341   assert( firstcallDefineResForm_ );
00342   firstcallDefineResForm_ = false;
00343     
00344   restype_ = TypeOfResidual;
00345   resnormtype_ = TypeOfNorm;
00346     
00347   // These conditions force the residual vector to be computed
00348   if (restype_==Explicit)
00349     resvecrequired_ = true;
00350     
00351   return(0);
00352 }
00353 
00354 template <class ScalarType, class MV, class OP> 
00355 int StatusTestResNorm<ScalarType,MV,OP>::DefineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00356                                                          MagnitudeType ScaleValue )
00357 {
00358     
00359   assert( firstcallDefineScaleForm_ );
00360   firstcallDefineScaleForm_ = false;
00361     
00362   scaletype_ = TypeOfScaling;
00363   scalenormtype_ = TypeOfNorm;
00364   scalevalue_ = ScaleValue;
00365     
00366   return(0);
00367 }
00368 
00369 template <class ScalarType, class MV, class OP>
00370 StatusType StatusTestResNorm<ScalarType,MV,OP>::CheckStatus( IterativeSolver<ScalarType,MV,OP>* iSolver )
00371 {
00372   int i;
00373   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00374   RefCountPtr<LinearProblem<ScalarType,MV,OP> > lp = iSolver->GetLinearProblem();
00375   // Compute scaling term (done once for each block that's being solved)
00376   if (firstcallCheckStatus_) {
00377     StatusType status = FirstcallCheckStatusSetup(iSolver);
00378     if(status==Failed) {
00379       status_ = Failed;
00380       return(status_);
00381     }
00382   }
00383   //
00384   // This section computes the norm of the residual vector
00385   //
00386   if ( cur_rhs_num_ != lp->GetRHSIndex() || cur_blksz_ != lp->GetCurrBlockSize() || cur_num_rhs_ != lp->GetNumToSolve() ) {
00387     //
00388     // We have moved on to the next rhs block
00389     //
00390     cur_rhs_num_ = lp->GetRHSIndex();
00391     cur_blksz_ = lp->GetCurrBlockSize();
00392     cur_num_rhs_ = lp->GetNumToSolve();
00393     //
00394   } else {
00395     //
00396     // We are in the same rhs block, return if we are converged
00397     //
00398     if (status_==Converged) { return status_; }
00399   }
00400   if (restype_==Implicit) {
00401     //
00402     // Get the native residual norms from the solver for this block of right-hand sides.
00403     // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
00404     // Otherwise the native residual is assumed to be stored in the resvector_.
00405     //
00406     std::vector<MagnitudeType> tmp_resvector( cur_blksz_ );
00407     RefCountPtr<const MV> residMV = iSolver->GetNativeResiduals( &tmp_resvector );     
00408     if ( residMV.get() != NULL ) { 
00409       tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
00410       MVT::MvNorm( *residMV, &tmp_resvector, resnormtype_ );    
00411       for (i=0; i<MVT::GetNumberVecs( *residMV ) && i<cur_num_rhs_; i++)
00412         resvector_[i+cur_rhs_num_] = tmp_resvector[i]; 
00413     } else {
00414       for (i=0; i<cur_num_rhs_; i++)
00415         resvector_[i+cur_rhs_num_] = tmp_resvector[i];
00416     }
00417   }
00418   else if (restype_==Explicit) {
00419     //
00420     // Request the true residual for this block of right-hand sides.
00421     // See if the linear problem manager has been updated before
00422     // asking for the true residual from the solver.
00423     //
00424     //
00425     if ( lp->IsSolutionUpdated() ) {
00426       const MV &cur_res = lp->GetCurrResVec();
00427       std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( cur_res ) );
00428       MVT::MvNorm( cur_res, &tmp_resvector, resnormtype_ );
00429       for (i=0; i<MVT::GetNumberVecs( cur_res ) && i<cur_num_rhs_; i++)
00430         resvector_[i+cur_rhs_num_] = tmp_resvector[i];
00431     } else {
00432       RefCountPtr<const MV> cur_soln = iSolver->GetCurrentSoln();
00433       const MV &cur_res = lp->GetCurrResVec( &*cur_soln );
00434       std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( cur_res ) );
00435       MVT::MvNorm( cur_res, &tmp_resvector, resnormtype_ );
00436       for (i=0; i<MVT::GetNumberVecs( cur_res ) && i<cur_num_rhs_; i++)
00437         resvector_[i+cur_rhs_num_] = tmp_resvector[i];      
00438     }
00439   }
00440   //
00441   // Compute the new linear system residuals for testing.
00442   // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
00443   //
00444   if ( scalevector_.size() > 0 ) {
00445     for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++) {
00446      
00447       // Scale the vector accordingly
00448       if ( scalevector_[i] != zero ) {
00449         // Don't intentionally divide by zero.
00450         testvector_[ i ] = resvector_[ i ] / scalevector_[ i ] / scalevalue_;
00451       } else {
00452         testvector_[ i ] = resvector_[ i ] / scalevalue_;
00453       }
00454     }
00455   }
00456   else {
00457     for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++)
00458       testvector_[ i ] = resvector_[ i ] / scalevalue_;
00459   } 
00460 
00461   // Check status of new linear system residuals
00462   status_ = Converged; // This may be set to unconverged or NaN.
00463   for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++) {
00464     // Check if any of the residuals are larger than the tolerance.
00465     if (testvector_[ i ] > tolerance_) {
00466       status_ = Unconverged;
00467       return(status_);
00468     } else if (testvector_[ i ] <= tolerance_) { 
00469       // do nothing.
00470     } else {
00471       status_ = NaN;            
00472       return(status_); // Return immediately if we detect a NaN.
00473     }
00474   } 
00475   
00476   // Return the current status
00477   return status_;
00478 }
00479 
00480 template <class ScalarType, class MV, class OP>
00481 ostream& StatusTestResNorm<ScalarType,MV,OP>::Print(ostream& os, int indent) const
00482 {
00483   for (int j = 0; j < indent; j ++)
00484     os << ' ';
00485   this->PrintStatus(os, status_);
00486   os << "(";
00487   os << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00488   os << ((restype_==Explicit) ? " Exp" : " Imp");
00489   os << " Res Vec) ";
00490   if (scaletype_!=None)
00491     os << "/ ";
00492   if (scaletype_==UserProvided)
00493     os << " (User Scale)";
00494   else {
00495     os << "(";
00496     os << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00497     if (scaletype_==NormOfInitRes)
00498       os << " Res0";
00499     else if (scaletype_==NormOfPrecInitRes)
00500       os << " Prec Res0";
00501     else
00502       os << " RHS ";
00503     os << ")";
00504   }
00505   if (status_==Unchecked)
00506     os << " Unchecked ( tol = " << tolerance_ << " ) "<<endl;
00507   else {
00508     os << endl;
00509     if(showMaxResNormOnly_) {
00510       const MagnitudeType maxRelRes = *std::max_element(
00511         testvector_.begin()+cur_rhs_num_,testvector_.begin()+cur_rhs_num_+cur_num_rhs_
00512         );
00513       for (int j = 0; j < indent + 13; j ++)
00514         os << ' ';
00515       os << "max{residual["<<cur_rhs_num_<<"..."<<cur_rhs_num_+cur_num_rhs_-1<<"]} = " << maxRelRes
00516          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << endl;
00517     }
00518     else {
00519       for ( int i=0; i<numrhs_; i++ ) {
00520         for (int j = 0; j < indent + 13; j ++)
00521           os << ' ';
00522         os << "residual [ " << i << " ] = " << testvector_[ i ];
00523         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << endl;
00524       }
00525     }
00526   }
00527   os << endl;
00528   return os;
00529 }
00530 
00531 template <class ScalarType, class MV, class OP>
00532 StatusType StatusTestResNorm<ScalarType,MV,OP>::FirstcallCheckStatusSetup( IterativeSolver<ScalarType,MV,OP>* iSolver )
00533 {
00534   int i;
00535   MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00536   MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00537   RefCountPtr<LinearProblem<ScalarType,MV,OP> > lp = iSolver->GetLinearProblem();
00538   // Compute scaling term (done once for each block that's being solved)
00539   if (firstcallCheckStatus_) {
00540     //
00541     // Get some current solver information.
00542     //
00543     firstcallCheckStatus_ = false;
00544     cur_rhs_num_ = lp->GetRHSIndex();
00545     cur_blksz_ = lp->GetCurrBlockSize();
00546     cur_num_rhs_ = lp->GetNumToSolve();
00547     //
00548     if (scaletype_== NormOfRHS) {
00549       const MV& rhs = *(lp->GetRHS());
00550       numrhs_ = MVT::GetNumberVecs( rhs );
00551       scalevector_.resize( numrhs_ );
00552       resvector_.resize( numrhs_ ); 
00553       testvector_.resize( numrhs_ );
00554       MVT::MvNorm( rhs, &scalevector_, scalenormtype_ );
00555     }
00556     else if (scaletype_==NormOfInitRes) {
00557       const MV &init_res = lp->GetInitResVec();
00558       numrhs_ = MVT::GetNumberVecs( init_res );
00559       scalevector_.resize( numrhs_ );
00560       resvector_.resize( numrhs_ ); 
00561       testvector_.resize( numrhs_ );
00562       MVT::MvNorm( init_res, &scalevector_, scalenormtype_ );
00563     }
00564     else if (scaletype_==NormOfPrecInitRes) {
00565       const MV& init_res = lp->GetInitResVec();
00566       numrhs_ = MVT::GetNumberVecs( init_res );
00567       scalevector_.resize( numrhs_ );
00568       resvector_.resize( numrhs_ ); 
00569       testvector_.resize( numrhs_ );
00570       RefCountPtr<MV> prec_init_res = MVT::Clone( init_res, numrhs_ );
00571       if (lp->ApplyLeftPrec( init_res, *prec_init_res ) != Undefined)
00572         MVT::MvNorm( *prec_init_res, &scalevector_, scalenormtype_ );
00573       else 
00574         MVT::MvNorm( init_res, &scalevector_, scalenormtype_ );
00575     }
00576 
00577     // Initialize the testvector.
00578     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00579 
00580     // Return an error if the scaling is zero.
00581     if (scalevalue_ == zero) {
00582       return Failed;
00583     }
00584   }
00585   return Unchecked;
00586 }
00587 
00588 } // end namespace Belos
00589 
00590 #endif /* BELOS_STATUS_TEST_RESNORM_H */

Generated on Thu Sep 18 12:30:12 2008 for Belos by doxygen 1.3.9.1