Belos Package Browser (Single Doxygen Collection) Development
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_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
00430     status_(Undefined),
00431     curBlksz_(0),
00432     curNumRHS_(0),
00433     curLSNum_(0),
00434     numrhs_(0),
00435     firstcallCheckStatus_(true),
00436     firstcallDefineResForm_(true),
00437     firstcallDefineScaleForm_(true),
00438     lossDetected_(false)
00439 {
00440   // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
00441   // the implicit residual vector.
00442 }
00443 
00444 template <class ScalarType, class MV, class OP>
00445 StatusTestImpResNorm<ScalarType,MV,OP>::~StatusTestImpResNorm()
00446 {}
00447 
00448 template <class ScalarType, class MV, class OP>
00449 void StatusTestImpResNorm<ScalarType,MV,OP>::reset()
00450 {
00451   status_ = Undefined;
00452   curBlksz_ = 0;
00453   curLSNum_ = 0;
00454   curLSIdx_.resize(0);
00455   numrhs_ = 0;
00456   ind_.resize(0);
00457   currTolerance_ = tolerance_;
00458   firstcallCheckStatus_ = true;
00459   lossDetected_ = false;
00460   curSoln_ = Teuchos::null;
00461 }
00462 
00463 template <class ScalarType, class MV, class OP>
00464 int StatusTestImpResNorm<ScalarType,MV,OP>::defineResForm( NormType TypeOfNorm )
00465 {
00466   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
00467         "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
00468   firstcallDefineResForm_ = false;
00469 
00470   resnormtype_ = TypeOfNorm;
00471 
00472   return(0);
00473 }
00474 
00475 template <class ScalarType, class MV, class OP>
00476 int StatusTestImpResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
00477                                                          MagnitudeType ScaleValue )
00478 {
00479   TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
00480         "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
00481   firstcallDefineScaleForm_ = false;
00482 
00483   scaletype_ = TypeOfScaling;
00484   scalenormtype_ = TypeOfNorm;
00485   scalevalue_ = ScaleValue;
00486 
00487   return(0);
00488 }
00489 
00490 template <class ScalarType, class MV, class OP>
00491 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
00492 checkStatus (Iteration<ScalarType,MV,OP>* iSolver)
00493 {
00494   using Teuchos::as;
00495   using Teuchos::RCP;
00496 
00497   const MagnitudeType zero = STM::zero ();
00498   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem ();
00499 
00500   // Compute scaling term (done once for each block that's being solved)
00501   if (firstcallCheckStatus_) {
00502     StatusType status = firstCallCheckStatusSetup (iSolver);
00503     if (status == Failed) {
00504       status_ = Failed;
00505       return status_;
00506     }
00507   }
00508 
00509   // mfh 23 Apr 2012: I don't know exactly what this code does.  It
00510   // has something to do with picking the block of right-hand sides
00511   // which we're currently checking.
00512   if (curLSNum_ != lp.getLSNumber ()) {
00513     //
00514     // We have moved on to the next rhs block
00515     //
00516     curLSNum_ = lp.getLSNumber();
00517     curLSIdx_ = lp.getLSIndex();
00518     curBlksz_ = (int)curLSIdx_.size();
00519     int validLS = 0;
00520     for (int i=0; i<curBlksz_; ++i) {
00521       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00522         validLS++;
00523     }
00524     curNumRHS_ = validLS;
00525     curSoln_ = Teuchos::null;
00526   } else {
00527     //
00528     // We are in the same rhs block, return if we are converged
00529     //
00530     if (status_ == Passed) {
00531       return status_;
00532     }
00533   }
00534 
00535   //
00536   // Get the "native" residual norms from the solver for this block of
00537   // right-hand sides.  If the solver's getNativeResiduals() method
00538   // actually returns a multivector, compute the norms of the columns
00539   // of the multivector explicitly.  Otherwise, we assume that
00540   // resvector_ contains the norms.
00541   //
00542   // Note that "compute the norms explicitly" doesn't necessarily mean
00543   // the "explicit" residual norms (in the sense discussed in this
00544   // class' documentation).  These are just some vectors returned by
00545   // the solver.  Some Krylov methods, like CG, compute a residual
00546   // vector "recursively."  This is an "implicit residual" in the
00547   // sense of this class' documentation.  It equals the explicit
00548   // residual in exact arithmetic, but due to rounding error, it is
00549   // usually different than the explicit residual.
00550   //
00551   // FIXME (mfh 23 Apr 2012) This method does _not_ respect the
00552   // OrthoManager used by the solver.
00553   //
00554   std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
00555   RCP<const MV> residMV = iSolver->getNativeResiduals (&tmp_resvector);
00556   if (! residMV.is_null ()) {
00557     // We got a multivector back.  Compute the norms explicitly.
00558     tmp_resvector.resize (MVT::GetNumberVecs (*residMV));
00559     MVT::MvNorm (*residMV, tmp_resvector, resnormtype_);
00560     typename std::vector<int>::iterator p = curLSIdx_.begin();
00561     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00562       // Check if this index is valid
00563       if (*p != -1) {
00564         resvector_[*p] = tmp_resvector[i];
00565       }
00566     }
00567   } else {
00568     typename std::vector<int>::iterator p = curLSIdx_.begin();
00569     for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00570       // Check if this index is valid
00571       if (*p != -1) {
00572         resvector_[*p] = tmp_resvector[i];
00573       }
00574     }
00575   }
00576   //
00577   // Scale the unscaled residual norms we computed or obtained above.
00578   //
00579   if (scalevector_.size () > 0) {
00580     // There are per-vector scaling factors to apply.
00581     typename std::vector<int>::iterator p = curLSIdx_.begin();
00582     for (; p<curLSIdx_.end(); ++p) {
00583       // Check if this index is valid
00584       if (*p != -1) {
00585         // Scale the vector accordingly
00586         if ( scalevector_[ *p ] != zero ) {
00587           // Don't intentionally divide by zero.
00588           testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
00589         } else {
00590           testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00591         }
00592       }
00593     }
00594   }
00595   else { // There are no per-vector scaling factors.
00596     typename std::vector<int>::iterator p = curLSIdx_.begin();
00597     for (; p<curLSIdx_.end(); ++p) {
00598       // Check if this index is valid
00599       if (*p != -1) {
00600         testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
00601       }
00602     }
00603   }
00604 
00605   // Count how many scaled residual norms (in testvector_) pass, using
00606   // the current tolerance (currTolerance_) rather than the original
00607   // tolerance (tolerance_).  If at least quorum_ of them pass, we
00608   // have a quorum for the whole test to pass.
00609   //
00610   // We also check here whether any of the scaled residual norms is
00611   // NaN, and throw an exception in that case.
00612   int have = 0;
00613   ind_.resize( curLSIdx_.size() );
00614   std::vector<int> lclInd( curLSIdx_.size() );
00615   typename std::vector<int>::iterator p = curLSIdx_.begin();
00616   for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
00617     // Check if this index is valid
00618     if (*p != -1) {
00619       if (testvector_[ *p ] > currTolerance_) {
00620         // The current residual norm doesn't pass.  Do nothing.
00621       } else if (testvector_[ *p ] <= currTolerance_) {
00622         ind_[have] = *p;
00623         lclInd[have] = i;
00624         have++; // Yay, the current residual norm passes!
00625       } else {
00626         // Throw an std::exception if the current residual norm is
00627         // NaN.  We know that it's NaN because it is not less than,
00628         // equal to, or greater than the current tolerance.  This is
00629         // only possible if either the residual norm or the current
00630         // tolerance is NaN; we assume the former.  We also mark the
00631         // test as failed, in case you want to catch the exception.
00632         status_ = Failed;
00633         TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "Belos::"
00634           "StatusTestImpResNorm::checkStatus(): One or more of the current "
00635           "implicit residual norms is NaN.");
00636       }
00637     }
00638   }
00639   // "have" is the number of residual norms that passed.
00640   ind_.resize(have);
00641   lclInd.resize(have);
00642 
00643   // Now check the exact residuals
00644   if (have) { // At least one residual norm has converged.
00645     //
00646     // Compute the explicit residual norm(s) from the current solution update.
00647     //
00648     RCP<MV> cur_update = iSolver->getCurrentUpdate ();
00649     curSoln_ = lp.updateSolution (cur_update);
00650     RCP<MV> cur_res = MVT::Clone (*curSoln_, MVT::GetNumberVecs (*curSoln_));
00651     lp.computeCurrResVec (&*cur_res, &*curSoln_);
00652     tmp_resvector.resize (MVT::GetNumberVecs (*cur_res));
00653     std::vector<MagnitudeType> tmp_testvector (have);
00654     MVT::MvNorm (*cur_res, tmp_resvector, resnormtype_);
00655 
00656     // Scale the explicit residual norm(s), just like the implicit norm(s).
00657     if ( scalevector_.size() > 0 ) {
00658       for (int i=0; i<have; ++i) {
00659         // Scale the vector accordingly
00660         if ( scalevector_[ ind_[i] ] != zero ) {
00661           // Don't intentionally divide by zero.
00662           tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevector_[ ind_[i] ] / scalevalue_;
00663         } else {
00664           tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00665         }
00666       }
00667     }
00668     else {
00669       for (int i=0; i<have; ++i) {
00670         tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
00671       }
00672     }
00673 
00674     //
00675     // Check whether the explicit residual norms also pass the
00676     // convergence test.  If not, check whether we want to try
00677     // iterating a little more to force both implicit and explicit
00678     // residual norms to pass.
00679     //
00680     int have2 = 0;
00681     for (int i = 0; i < have; ++i) {
00682       // testvector_ contains the implicit (i.e., recursive, computed
00683       // by the algorithm) (possibly scaled) residuals.  All of these
00684       // pass the convergence test.
00685       //
00686       // tmp_testvector contains the explicit (i.e., ||B-AX||)
00687       // (possibly scaled) residuals.  We're checking whether these
00688       // pass as well.  The explicit residual norms only have to meet
00689       // the _original_ tolerance (tolerance_), not the current
00690       // tolerance (currTolerance_).
00691       if (tmp_testvector[i] <= tolerance_) {
00692         ind_[have2] = ind_[i];
00693         have2++; // This right-hand side has converged.
00694       }
00695       else {
00696         // Absolute difference between the current explicit and
00697         // implicit residual norm.
00698         const MagnitudeType diff = STM::magnitude (testvector_[ind_[i]] - tmp_testvector[i]);
00699         if (diff > currTolerance_) {
00700           // If the above difference is bigger than the current
00701           // convergence tolerance, report a loss of accuracy, but
00702           // mark this right-hand side converged.  (The latter tells
00703           // users not to iterate further on this right-hand side,
00704           // since it probably won't do much good.)  Note that the
00705           // current tolerance may have been changed by the branch
00706           // below in a previous call to this method.
00707           lossDetected_ = true;
00708           ind_[have2] = ind_[i];
00709           have2++; // Count this right-hand side as converged.
00710         }
00711         else {
00712           // Otherwise, the explicit and implicit residuals are pretty
00713           // close together, and the implicit residual has converged,
00714           // but the explicit residual hasn't converged.  Reduce the
00715           // convergence tolerance by some formula related to the
00716           // difference, and keep iterating.
00717           //
00718           // mfh 23 Apr 2012: I have no idea why the currTolerance_
00719           // update formula is done the way it's done below.  It
00720           // doesn't make sense to me.  It definitely makes the
00721           // current tolerance smaller, though, which is what we want.
00722 
00723           // We define these constants in this way, rather than simply
00724           // writing 1.5 resp. 0.1, to ensure no rounding error from
00725           // translating from float to MagnitudeType.  Remember that
00726           // 0.1 doesn't have an exact representation in binary finite
00727           // floating-point arithmetic.
00728           const MagnitudeType onePointFive = as<MagnitudeType>(3) / as<MagnitudeType> (2);
00729           const MagnitudeType oneTenth = STM::one () / as<MagnitudeType> (10);
00730 
00731           currTolerance_ = currTolerance_ - onePointFive * diff;
00732           while (currTolerance_ < STM::zero ()) {
00733             currTolerance_ += oneTenth * diff;
00734           }
00735         }
00736       }
00737     }
00738     have = have2;
00739     ind_.resize(have);
00740   }
00741 
00742   // Check whether we've met the quorum of vectors necessary for the
00743   // whole test to pass.
00744   int need = (quorum_ == -1) ? curNumRHS_: quorum_;
00745   status_ = (have >= need) ? Passed : Failed;
00746 
00747   // Return the current status
00748   return status_;
00749 }
00750 
00751 template <class ScalarType, class MV, class OP>
00752 void StatusTestImpResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00753 {
00754   for (int j = 0; j < indent; j ++)
00755     os << ' ';
00756   printStatus(os, status_);
00757   os << resFormStr();
00758   if (status_==Undefined)
00759     os << ", tol = " << tolerance_ << std::endl;
00760   else {
00761     os << std::endl;
00762     if(showMaxResNormOnly_ && curBlksz_ > 1) {
00763       const MagnitudeType maxRelRes = *std::max_element(
00764         testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
00765         );
00766       for (int j = 0; j < indent + 13; j ++)
00767         os << ' ';
00768       os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
00769          << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
00770     }
00771     else {
00772       for ( int i=0; i<numrhs_; i++ ) {
00773         for (int j = 0; j < indent + 13; j ++)
00774           os << ' ';
00775         os << "residual [ " << i << " ] = " << testvector_[ i ];
00776         os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << std::endl;
00777       }
00778     }
00779   }
00780   os << std::endl;
00781 }
00782 
00783 template <class ScalarType, class MV, class OP>
00784 void StatusTestImpResNorm<ScalarType,MV,OP>::printStatus(std::ostream& os, StatusType type) const
00785 {
00786   os << std::left << std::setw(13) << std::setfill('.');
00787   switch (type) {
00788   case  Passed:
00789     os << "Converged";
00790     break;
00791   case  Failed:
00792     if (lossDetected_)
00793       os << "Unconverged (LoA)";
00794     else
00795       os << "Unconverged";
00796     break;
00797   case  Undefined:
00798   default:
00799     os << "**";
00800     break;
00801   }
00802   os << std::left << std::setfill(' ');
00803     return;
00804 }
00805 
00806 template <class ScalarType, class MV, class OP>
00807 StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
00808 firstCallCheckStatusSetup (Iteration<ScalarType,MV,OP>* iSolver)
00809 {
00810   int i;
00811   const MagnitudeType zero = STM::zero ();
00812   const MagnitudeType one = STM::one ();
00813   const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
00814   // Compute scaling term (done once for each block that's being solved)
00815   if (firstcallCheckStatus_) {
00816     //
00817     // Get some current solver information.
00818     //
00819     firstcallCheckStatus_ = false;
00820 
00821     if (scaletype_== NormOfRHS) {
00822       Teuchos::RCP<const MV> rhs = lp.getRHS();
00823       numrhs_ = MVT::GetNumberVecs( *rhs );
00824       scalevector_.resize( numrhs_ );
00825       MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
00826     }
00827     else if (scaletype_==NormOfInitRes) {
00828       Teuchos::RCP<const MV> init_res = lp.getInitResVec();
00829       numrhs_ = MVT::GetNumberVecs( *init_res );
00830       scalevector_.resize( numrhs_ );
00831       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00832     }
00833     else if (scaletype_==NormOfPrecInitRes) {
00834       Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
00835       numrhs_ = MVT::GetNumberVecs( *init_res );
00836       scalevector_.resize( numrhs_ );
00837       MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
00838     }
00839     else {
00840       numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
00841     }
00842 
00843     resvector_.resize( numrhs_ );
00844     testvector_.resize( numrhs_ );
00845 
00846     curLSNum_ = lp.getLSNumber();
00847     curLSIdx_ = lp.getLSIndex();
00848     curBlksz_ = (int)curLSIdx_.size();
00849     int validLS = 0;
00850     for (i=0; i<curBlksz_; ++i) {
00851       if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
00852         validLS++;
00853     }
00854     curNumRHS_ = validLS;
00855     //
00856     // Initialize the testvector.
00857     for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
00858 
00859     // Return an error if the scaling is zero.
00860     if (scalevalue_ == zero) {
00861       return Failed;
00862     }
00863   }
00864   return Undefined;
00865 }
00866 
00867 } // end namespace Belos
00868 
00869 #endif /* BELOS_STATUS_TEST_IMPRESNORM_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines