Belos Version of the Day
BelosStatusTestImpResNorm.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_IMPRESNORM_H
00043 #define BELOS_STATUS_TEST_IMPRESNORM_H
00044 
00050 #include "BelosStatusTestResNorm.hpp"
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosMultiVecTraits.hpp"
00053 #include "Teuchos_as.hpp"
00054 
00101 namespace Belos {
00102 
00103 template <class ScalarType, class MV, class OP>
00104 class StatusTestImpResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
00105 public:
00107   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00108 
00109 private:
00111 
00112   typedef Teuchos::ScalarTraits<ScalarType> STS;
00113   typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00114   typedef MultiVecTraits<ScalarType,MV> MVT;
00116 
00117 public:
00119 
00120 
00132   StatusTestImpResNorm (MagnitudeType Tolerance,
00133                         int quorum = -1,
00134                         bool showMaxResNormOnly = false);
00135 
00137   virtual ~StatusTestImpResNorm();
00138 
00140 
00141 
00142 
00144 
00150   int defineResForm( NormType TypeOfNorm);
00151 
00153 
00173   int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
00174 
00176 
00179   int setTolerance (MagnitudeType tolerance) {
00180     tolerance_ = tolerance;
00181     return 0;
00182   }
00183 
00186   int setQuorum (int quorum) {
00187     quorum_ = quorum;
00188     return 0;
00189   }
00190 
00192   int setShowMaxResNormOnly (bool showMaxResNormOnly) {
00193     showMaxResNormOnly_ = showMaxResNormOnly;
00194     return 0;
00195   }
00196 
00198 
00200 
00201 
00202 
00208   StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
00209 
00211   StatusType getStatus() const {return(status_);};
00213 
00215 
00216 
00218   void reset();
00219 
00221 
00223 
00224 
00226   void print(std::ostream& os, int indent = 0) const;
00227 
00229   void printStatus(std::ostream& os, StatusType type) const;
00231 
00233 
00234 
00236   Teuchos::RCP<MV> getSolution() { return curSoln_; }
00237 
00240   int getQuorum() const { return quorum_; }
00241 
00243   bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
00244 
00246   std::vector<int> convIndices() { return ind_; }
00247 
00255   MagnitudeType getTolerance () const {
00256     return tolerance_;
00257   }
00258 
00266   MagnitudeType getCurrTolerance () const {
00267     return currTolerance_;
00268   }
00269 
00271   const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
00272 
00274   const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
00275 
00277   const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
00278 
00280   bool getLOADetected() const { return lossDetected_; }
00282 
00283 
00286 
00292   StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
00294 
00297 
00299   std::string description() const
00300   {
00301     std::ostringstream oss;
00302     oss << "Belos::StatusTestImpResNorm<>: " << resFormStr();
00303     oss << ", tol = " << tolerance_;
00304     return oss.str();
00305   }
00307 
00308  protected:
00309 
00310  private:
00311 
00313 
00314 
00315   std::string resFormStr() const
00316   {
00317     std::ostringstream oss;
00318     oss << "(";
00319     oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00320     oss << " Res Vec) ";
00321 
00322     // If there is no residual scaling, return current string.
00323     if (scaletype_!=None)
00324     {
00325       // Insert division sign.
00326       oss << "/ ";
00327 
00328       // Determine output string for scaling, if there is any.
00329       if (scaletype_==UserProvided)
00330         oss << " (User Scale)";
00331       else {
00332         oss << "(";
00333         oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
00334         if (scaletype_==NormOfInitRes)
00335           oss << " Res0";
00336         else if (scaletype_==NormOfPrecInitRes)
00337           oss << " Prec Res0";
00338         else
00339           oss << " RHS ";
00340         oss << ")";
00341       }
00342     }
00343 
00344     return oss.str();
00345   }
00346 
00348 
00349 
00351   MagnitudeType tolerance_, currTolerance_;
00352 
00354   int quorum_;
00355 
00357   bool showMaxResNormOnly_;
00358 
00360   NormType resnormtype_;
00361 
00363   ScaleType scaletype_;
00364 
00366   NormType scalenormtype_;
00367 
00369   MagnitudeType scalevalue_;
00370 
00372   std::vector<MagnitudeType> scalevector_;
00373 
00375   std::vector<MagnitudeType> resvector_;
00376 
00378   std::vector<MagnitudeType> testvector_;
00379 
00381   Teuchos::RCP<MV> curSoln_;
00382 
00384   std::vector<int> ind_;
00385 
00387   StatusType status_;
00388 
00390   int curBlksz_;
00391 
00393   int curNumRHS_;
00394 
00396   std::vector<int> curLSIdx_;
00397 
00399   int curLSNum_;
00400 
00402   int numrhs_;
00403 
00405   bool firstcallCheckStatus_;
00406 
00408   bool firstcallDefineResForm_;
00409 
00411   bool firstcallDefineScaleForm_;
00412 
00414   bool lossDetected_;
00416 
00417 };
00418 
00419 template <class ScalarType, class MV, class OP>
00420 StatusTestImpResNorm<ScalarType,MV,OP>::
00421 StatusTestImpResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
00422   : tolerance_(Tolerance),
00423     currTolerance_(Tolerance),
00424     quorum_(quorum),
00425     showMaxResNormOnly_(showMaxResNormOnly),
00426     resnormtype_(TwoNorm),
00427     scaletype_(NormOfInitRes),
00428     scalenormtype_(TwoNorm),
00429     scalevalue_(1.0),
00430     status_(Undefined),
00431     curBlksz_(0),
00432     curLSNum_(0),
00433     numrhs_(0),
00434     firstcallCheckStatus_(true),
00435     firstcallDefineResForm_(true),
00436     firstcallDefineScaleForm_(true),
00437     lossDetected_(false)
00438 {
00439   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00440   // the implicit residual vector.
00441 }
00442 
00443 template <class ScalarType, class MV, class OP>
00444 StatusTestImpResNorm<ScalarType,MV,OP>::~StatusTestImpResNorm()
00445 {}
00446 
00447 template <class ScalarType, class MV, class OP>
00448 void StatusTestImpResNorm<ScalarType,MV,OP>::reset()
00449 {
00450   status_ = Undefined;
00451   curBlksz_ = 0;
00452   curLSNum_ = 0;
00453   curLSIdx_.resize(0);
00454   numrhs_ = 0;
00455   ind_.resize(0);
00456   currTolerance_ = tolerance_;
00457   firstcallCheckStatus_ = true;
00458   lossDetected_ = false;
00459   curSoln_ = Teuchos::null;
00460 }
00461 
00462 template <class ScalarType, class MV, class OP>
00463 int StatusTestImpResNorm<ScalarType,MV,OP>::defineResForm( NormType TypeOfNorm )
00464 {
00465   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00466         "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
00467   firstcallDefineResForm_ = false;
00468 
00469   resnormtype_ = TypeOfNorm;
00470 
00471   return(0);
00472 }
00473 
00474 template <class ScalarType, class MV, class OP>
00475 int StatusTestImpResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00476                                                          MagnitudeType ScaleValue )
00477 {
00478   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00479         "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
00480   firstcallDefineScaleForm_ = false;
00481 
00482   scaletype_ = TypeOfScaling;
00483   scalenormtype_ = TypeOfNorm;
00484   scalevalue_ = ScaleValue;
00485 
00486   return(0);
00487 }
00488 
00489 template <class ScalarType, class MV, class OP>
00490 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
00491 checkStatus (Iteration<ScalarType,MV,OP>* iSolver)
00492 {
00493   using Teuchos::as;
00494   using Teuchos::RCP;
00495 
00496   const MagnitudeType zero = STM::zero ();
00497   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem ();
00498 
00499   // Compute scaling term (done once for each block that's being solved)
00500   if (firstcallCheckStatus_) {
00501     StatusType status = firstCallCheckStatusSetup (iSolver);
00502     if (status == Failed) {
00503       status_ = Failed;
00504       return status_;
00505     }
00506   }
00507 
00508   // mfh 23 Apr 2012: I don't know exactly what this code does.  It
00509   // has something to do with picking the block of right-hand sides
00510   // which we're currently checking.
00511   if (curLSNum_ != lp.getLSNumber ()) {
00512     //
00513     // We have moved on to the next rhs block
00514     //
00515     curLSNum_ = lp.getLSNumber();
00516     curLSIdx_ = lp.getLSIndex();
00517     curBlksz_ = (int)curLSIdx_.size();
00518     int validLS = 0;
00519     for (int i=0; i<curBlksz_; ++i) {
00520       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00521         validLS++;
00522     }
00523     curNumRHS_ = validLS;
00524     curSoln_ = Teuchos::null;
00525   } else {
00526     //
00527     // We are in the same rhs block, return if we are converged
00528     //
00529     if (status_ == Passed) {
00530       return status_;
00531     }
00532   }
00533 
00534   //
00535   // Get the "native" residual norms from the solver for this block of
00536   // right-hand sides.  If the solver's getNativeResiduals() method
00537   // actually returns a multivector, compute the norms of the columns
00538   // of the multivector explicitly.  Otherwise, we assume that
00539   // resvector_ contains the norms.
00540   //
00541   // Note that "compute the norms explicitly" doesn't necessarily mean
00542   // the "explicit" residual norms (in the sense discussed in this
00543   // class' documentation).  These are just some vectors returned by
00544   // the solver.  Some Krylov methods, like CG, compute a residual
00545   // vector "recursively."  This is an "implicit residual" in the
00546   // sense of this class' documentation.  It equals the explicit
00547   // residual in exact arithmetic, but due to rounding error, it is
00548   // usually different than the explicit residual.
00549   //
00550   // FIXME (mfh 23 Apr 2012) This method does _not_ respect the
00551   // OrthoManager used by the solver.
00552   //
00553   std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00554   RCP<const MV> residMV = iSolver->getNativeResiduals (&tmp_resvector);
00555   if (! residMV.is_null ()) {
00556     // We got a multivector back.  Compute the norms explicitly.
00557     tmp_resvector.resize (MVT::GetNumberVecs (*residMV));
00558     MVT::MvNorm (*residMV, tmp_resvector, resnormtype_);
00559     typename std::vector<int>::iterator p = curLSIdx_.begin();
00560     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00561       // Check if this index is valid
00562       if (*p != -1) {
00563         resvector_[*p] = tmp_resvector[i];
00564       }
00565     }
00566   } else {
00567     typename std::vector<int>::iterator p = curLSIdx_.begin();
00568     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00569       // Check if this index is valid
00570       if (*p != -1) {
00571         resvector_[*p] = tmp_resvector[i];
00572       }
00573     }
00574   }
00575   //
00576   // Scale the unscaled residual norms we computed or obtained above.
00577   //
00578   if (scalevector_.size () > 0) {
00579     // There are per-vector scaling factors to apply.
00580     typename std::vector<int>::iterator p = curLSIdx_.begin();
00581     for (; p<curLSIdx_.end(); ++p) {
00582       // Check if this index is valid
00583       if (*p != -1) {
00584         // Scale the vector accordingly
00585         if ( scalevector_[ *p ] != zero ) {
00586           // Don't intentionally divide by zero.
00587           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00588         } else {
00589           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00590         }
00591       }
00592     }
00593   }
00594   else { // There are no per-vector scaling factors.
00595     typename std::vector<int>::iterator p = curLSIdx_.begin();
00596     for (; p<curLSIdx_.end(); ++p) {
00597       // Check if this index is valid
00598       if (*p != -1) {
00599         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00600       }
00601     }
00602   }
00603 
00604   // Count how many scaled residual norms (in testvector_) pass, using
00605   // the current tolerance (currTolerance_) rather than the original
00606   // tolerance (tolerance_).  If at least quorum_ of them pass, we
00607   // have a quorum for the whole test to pass.
00608   //
00609   // We also check here whether any of the scaled residual norms is
00610   // NaN, and throw an exception in that case.
00611   int have = 0;
00612   ind_.resize( curLSIdx_.size() );
00613   std::vector<int> lclInd( curLSIdx_.size() );
00614   typename std::vector<int>::iterator p = curLSIdx_.begin();
00615   for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00616     // Check if this index is valid
00617     if (*p != -1) {
00618       if (testvector_[ *p ] > currTolerance_) {
00619         // The current residual norm doesn't pass.  Do nothing.
00620       } else if (testvector_[ *p ] <= currTolerance_) {
00621         ind_[have] = *p;
00622         lclInd[have] = i;
00623         have++; // Yay, the current residual norm passes!
00624       } else {
00625         // Throw an std::exception if the current residual norm is
00626         // NaN.  We know that it's NaN because it is not less than,
00627         // equal to, or greater than the current tolerance.  This is
00628         // only possible if either the residual norm or the current
00629         // tolerance is NaN; we assume the former.  We also mark the
00630         // test as failed, in case you want to catch the exception.
00631         status_ = Failed;
00632         TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "Belos::"
00633           "StatusTestImpResNorm::checkStatus(): One or more of the current "
00634           "implicit residual norms is NaN.");
00635       }
00636     }
00637   }
00638   // "have" is the number of residual norms that passed.
00639   ind_.resize(have);
00640   lclInd.resize(have);
00641 
00642   // Now check the exact residuals
00643   if (have) { // At least one residual norm has converged.
00644     //
00645     // Compute the explicit residual norm(s) from the current solution update.
00646     //
00647     RCP<MV> cur_update = iSolver->getCurrentUpdate ();
00648     curSoln_ = lp.updateSolution (cur_update);
00649     RCP<MV> cur_res = MVT::Clone (*curSoln_, MVT::GetNumberVecs (*curSoln_));
00650     lp.computeCurrResVec (&*cur_res, &*curSoln_);
00651     tmp_resvector.resize (MVT::GetNumberVecs (*cur_res));
00652     std::vector<MagnitudeType> tmp_testvector (have);
00653     MVT::MvNorm (*cur_res, tmp_resvector, resnormtype_);
00654 
00655     // Scale the explicit residual norm(s), just like the implicit norm(s).
00656     if ( scalevector_.size() > 0 ) {
00657       for (int i=0; i<have; ++i) {
00658         // Scale the vector accordingly
00659         if ( scalevector_[ ind_[i] ] != zero ) {
00660           // Don't intentionally divide by zero.
00661           tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevector_[ ind_[i] ] / scalevalue_;
00662         } else {
00663           tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00664         }
00665       }
00666     }
00667     else {
00668       for (int i=0; i<have; ++i) {
00669         tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00670       }
00671     }
00672 
00673     //
00674     // Check whether the explicit residual norms also pass the
00675     // convergence test.  If not, check whether we want to try
00676     // iterating a little more to force both implicit and explicit
00677     // residual norms to pass.
00678     //
00679     int have2 = 0;
00680     for (int i = 0; i < have; ++i) {
00681       // testvector_ contains the implicit (i.e., recursive, computed
00682       // by the algorithm) (possibly scaled) residuals.  All of these
00683       // pass the convergence test.
00684       //
00685       // tmp_testvector contains the explicit (i.e., ||B-AX||)
00686       // (possibly scaled) residuals.  We're checking whether these
00687       // pass as well.  The explicit residual norms only have to meet
00688       // the _original_ tolerance (tolerance_), not the current
00689       // tolerance (currTolerance_).
00690       if (tmp_testvector[i] <= tolerance_) {
00691         ind_[have2] = ind_[i];
00692         have2++; // This right-hand side has converged.
00693       }
00694       else {
00695         // Absolute difference between the current explicit and
00696         // implicit residual norm.
00697         const MagnitudeType diff = STM::magnitude (testvector_[ind_[i]] - tmp_testvector[i]);
00698         if (diff > currTolerance_) {
00699           // If the above difference is bigger than the current
00700           // convergence tolerance, report a loss of accuracy, but
00701           // mark this right-hand side converged.  (The latter tells
00702           // users not to iterate further on this right-hand side,
00703           // since it probably won't do much good.)  Note that the
00704           // current tolerance may have been changed by the branch
00705           // below in a previous call to this method.
00706           lossDetected_ = true;
00707           ind_[have2] = ind_[i];
00708           have2++; // Count this right-hand side as converged.
00709         }
00710         else {
00711           // Otherwise, the explicit and implicit residuals are pretty
00712           // close together, and the implicit residual has converged,
00713           // but the explicit residual hasn't converged.  Reduce the
00714           // convergence tolerance by some formula related to the
00715           // difference, and keep iterating.
00716           //
00717           // mfh 23 Apr 2012: I have no idea why the currTolerance_
00718           // update formula is done the way it's done below.  It
00719           // doesn't make sense to me.  It definitely makes the
00720           // current tolerance smaller, though, which is what we want.
00721 
00722           // We define these constants in this way, rather than simply
00723           // writing 1.5 resp. 0.1, to ensure no rounding error from
00724           // translating from float to MagnitudeType.  Remember that
00725           // 0.1 doesn't have an exact representation in binary finite
00726           // floating-point arithmetic.
00727           const MagnitudeType onePointFive = as<MagnitudeType>(3) / as<MagnitudeType> (2);
00728           const MagnitudeType oneTenth = STM::one () / as<MagnitudeType> (10);
00729 
00730           currTolerance_ = currTolerance_ - onePointFive * diff;
00731           while (currTolerance_ < STM::zero ()) {
00732             currTolerance_ += oneTenth * diff;
00733           }
00734         }
00735       }
00736     }
00737     have = have2;
00738     ind_.resize(have);
00739   }
00740 
00741   // Check whether we've met the quorum of vectors necessary for the
00742   // whole test to pass.
00743   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00744   status_ = (have >= need) ? Passed : Failed;
00745 
00746   // Return the current status
00747   return status_;
00748 }
00749 
00750 template <class ScalarType, class MV, class OP>
00751 void StatusTestImpResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00752 {
00753   for (int j = 0; j < indent; j ++)
00754     os << ' ';
00755   printStatus(os, status_);
00756   os << resFormStr();
00757   if (status_==Undefined)
00758     os << ", tol = " << tolerance_ << std::endl;
00759   else {
00760     os << std::endl;
00761     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00762       const MagnitudeType maxRelRes = *std::max_element(
00763         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00764         );
00765       for (int j = 0; j < indent + 13; j ++)
00766         os << ' ';
00767       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00768          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00769     }
00770     else {
00771       for ( int i=0; i<numrhs_; i++ ) {
00772         for (int j = 0; j < indent + 13; j ++)
00773           os << ' ';
00774         os << "residual [ " << i << " ] = " << testvector_[ i ];
00775         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00776       }
00777     }
00778   }
00779   os << std::endl;
00780 }
00781 
00782 template <class ScalarType, class MV, class OP>
00783 void StatusTestImpResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const
00784 {
00785   os << std::left << std::setw(13) << std::setfill('.');
00786   switch (type) {
00787   case  Passed:
00788     os << "Converged";
00789     break;
00790   case  Failed:
00791     if (lossDetected_)
00792       os << "Unconverged (LoA)";
00793     else
00794       os << "Unconverged";
00795     break;
00796   case  Undefined:
00797   default:
00798     os << "**";
00799     break;
00800   }
00801   os << std::left << std::setfill(' ');
00802     return;
00803 }
00804 
00805 template <class ScalarType, class MV, class OP>
00806 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
00807 firstCallCheckStatusSetup (Iteration<ScalarType,MV,OP>* iSolver)
00808 {
00809   int i;
00810   const MagnitudeType zero = STM::zero ();
00811   const MagnitudeType one = STM::one ();
00812   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00813   // Compute scaling term (done once for each block that's being solved)
00814   if (firstcallCheckStatus_) {
00815     //
00816     // Get some current solver information.
00817     //
00818     firstcallCheckStatus_ = false;
00819 
00820     if (scaletype_== NormOfRHS) {
00821       Teuchos::RCP<const MV> rhs = lp.getRHS();
00822       numrhs_ = MVT::GetNumberVecs( *rhs );
00823       scalevector_.resize( numrhs_ );
00824       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00825     }
00826     else if (scaletype_==NormOfInitRes) {
00827       Teuchos::RCP<const MV> init_res = lp.getInitResVec();
00828       numrhs_ = MVT::GetNumberVecs( *init_res );
00829       scalevector_.resize( numrhs_ );
00830       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00831     }
00832     else if (scaletype_==NormOfPrecInitRes) {
00833       Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
00834       numrhs_ = MVT::GetNumberVecs( *init_res );
00835       scalevector_.resize( numrhs_ );
00836       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00837     }
00838     else {
00839       numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
00840     }
00841 
00842     resvector_.resize( numrhs_ );
00843     testvector_.resize( numrhs_ );
00844 
00845     curLSNum_ = lp.getLSNumber();
00846     curLSIdx_ = lp.getLSIndex();
00847     curBlksz_ = (int)curLSIdx_.size();
00848     int validLS = 0;
00849     for (i=0; i<curBlksz_; ++i) {
00850       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00851         validLS++;
00852     }
00853     curNumRHS_ = validLS;
00854     //
00855     // Initialize the testvector.
00856     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00857 
00858     // Return an error if the scaling is zero.
00859     if (scalevalue_ == zero) {
00860       return Failed;
00861     }
00862   }
00863   return Undefined;
00864 }
00865 
00866 } // end namespace Belos
00867 
00868 #endif /* BELOS_STATUS_TEST_IMPRESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines