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