Belos Version of the Day
BelosLSQRStatusTest.hpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //                 Belos: Block Linear Solvers Package
00006 //                  Copyright 2004 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #ifndef BELOS_LSQR_STATUS_TEST_HPP
00045 #define BELOS_LSQR_STATUS_TEST_HPP
00046 
00052 #include "BelosStatusTest.hpp"
00053 #include "BelosLSQRIter.hpp"
00054 
00061 namespace Belos {
00062 
00063 
00064 template <class ScalarType, class MV, class OP>
00065 class LSQRStatusTest: public Belos::StatusTest<ScalarType,MV,OP> {
00066 
00067 public:
00068 
00069   // Convenience typedefs
00070   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00071   typedef typename SCT::magnitudeType MagnitudeType;
00072   typedef Belos::MultiVecTraits<ScalarType,MV>  MVT;
00073 
00075 
00076 
00078 
00083   LSQRStatusTest( MagnitudeType condMax = 0.0, 
00084                   int term_iter_max = 1, 
00085                   MagnitudeType rel_rhs_err = 0.0, 
00086                   MagnitudeType rel_mat_err = 0.0 );
00087 
00089   virtual ~LSQRStatusTest();
00091 
00093 
00094 
00096 
00098   Belos::StatusType checkStatus(Belos::Iteration<ScalarType,MV,OP> *iSolver );
00099 
00101   Belos::StatusType getStatus() const {return(status_);}
00102 
00104 
00106 
00107 
00109   void reset();
00110 
00112   int setCondLim(MagnitudeType condMax) {
00113     condMax_ = condMax; 
00114     rcondMin_ = (condMax > 0) ? (Teuchos::ScalarTraits< MagnitudeType >::one() / condMax) : Teuchos::ScalarTraits< MagnitudeType >::eps(); 
00115     return(0);}
00116 
00117   int setTermIterMax(int term_iter_max) {
00118     term_iter_max_ = term_iter_max;
00119     if (term_iter_max_ < 1)
00120       term_iter_max_ = 1;
00121     return(0);}
00122 
00123   int setRelRhsErr(MagnitudeType rel_rhs_err) {
00124     rel_rhs_err_ = rel_rhs_err;
00125     return(0);}
00126 
00127   int setRelMatErr(MagnitudeType rel_mat_err) {
00128     rel_mat_err_ = rel_mat_err;
00129     return(0);}
00130 
00132 
00134 
00135 
00137   MagnitudeType getCondMaxLim() const {return(condMax_);}
00138 
00140   int getTermIterMax() const {return(term_iter_max_);}
00141 
00143   MagnitudeType getRelRhsErr() const {return(rel_rhs_err_);}
00144 
00146   MagnitudeType getMatErr() const {return(rel_mat_err_);}
00147 
00149   MagnitudeType getMatCondNum() const {return(matCondNum_);}
00150 
00152   MagnitudeType getMatNorm() const {return(matNorm_);}
00153 
00155   int getTermIter() const { return term_iter_; }
00156 
00158   MagnitudeType getResidNorm() const {return(resNorm_);}
00159 
00161   MagnitudeType getLSResidNorm() const {return(matResNorm_);}
00163 
00164 
00166 
00167 
00169   void print(std::ostream& os, int indent = 0) const;
00170 
00172   void printStatus(std::ostream& os, Belos::StatusType type) const;
00173 
00175 
00178 
00181   Belos::StatusType firstCallCheckStatusSetup(Belos::Iteration<ScalarType,MV,OP>* iSolver);
00183 
00186 
00188   std::string description() const
00189   {
00190     std::ostringstream oss;
00191     oss << "LSQRStatusTest<>: [ limit of condition number = " << condMax_ << " ]";
00192     return oss.str();
00193   }
00195 
00196 private:
00197 
00199 
00200 
00202   MagnitudeType condMax_;
00203 
00205   int term_iter_max_;
00206 
00208   MagnitudeType rel_rhs_err_;
00209 
00211   MagnitudeType rel_mat_err_;
00212 
00214   MagnitudeType rcondMin_;
00215 
00217   Belos::StatusType status_;
00218 
00219   // term_iter_ records the number of consecutive "successful" iterations.
00220   // convergence requires that term_iter_max consecutive iterates satisfy the other convergence tests
00221   int term_iter_;
00222 
00223   // condition number of the operator 
00224   MagnitudeType matCondNum_;
00225 
00226   // Frobenius norm of the operator 
00227   MagnitudeType matNorm_;
00228 
00229   // residual norm for the linear system
00230   MagnitudeType resNorm_;
00231 
00232   // least squares residual, operator^Transpose * residual
00233   MagnitudeType matResNorm_;
00234 
00236 
00237 };
00238 
00239 template <class ScalarType, class MV, class OP>
00240 LSQRStatusTest<ScalarType,MV,OP>::LSQRStatusTest( MagnitudeType condMax /* = 0 */, int term_iter_max /* = 1 */, MagnitudeType rel_rhs_err /* = 0 */, MagnitudeType rel_mat_err /* = 0 */)
00241   : condMax_(condMax),
00242     term_iter_max_ (term_iter_max),
00243     rel_rhs_err_ (rel_rhs_err),
00244     rel_mat_err_ (rel_mat_err),
00245     status_ (Belos::Undefined),
00246     term_iter_ (0),
00247     matCondNum_ ( Teuchos::ScalarTraits<MagnitudeType>::one() ),
00248     matNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
00249     resNorm_  ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
00250     matResNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() )
00251 {}
00252 
00253 template <class ScalarType, class MV, class OP>
00254 LSQRStatusTest<ScalarType,MV,OP>::~LSQRStatusTest()
00255 {}
00256 
00257 template <class ScalarType, class MV, class OP>
00258 void LSQRStatusTest<ScalarType,MV,OP>::reset()
00259 {
00260   status_ = Belos::Undefined;
00261 }
00262 
00263 template <class ScalarType, class MV, class OP>
00264 Belos::StatusType LSQRStatusTest<ScalarType,MV,OP>::checkStatus( Belos::Iteration<ScalarType,MV,OP>* iSolver) 
00265 {
00266   const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00267   const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00268   if (condMax_ > MTzero ) 
00269     {
00270   rcondMin_ = MTone / condMax_;
00271     }
00272   else 
00273     {
00274   rcondMin_ = Teuchos::ScalarTraits< MagnitudeType >::eps();
00275     }
00276 
00277   bool termIterFlag = false;
00278   LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
00279   LSQRIterationState< ScalarType, MV > state = solver->getState();
00280   //
00281   //   LSQR solves a least squares problem.  A converged preconditioned residual norm
00282   // suffices for convergence, but is not necessary.  LSQR sometimes returns a larger
00283   // relative residual norm than what would have been returned by a linear solver.
00284   // This section evaluates three stopping criteria.  In the Solver Manager, this test
00285   // is combined with a generic number of iteration test.
00286   //   If the linear system includes a preconditioner, then the least squares problem
00287   // is solved for the preconditioned linear system.  Preconditioning changes the least
00288   // squares problem (in the sense of changing the norms), and the solution depends
00289   // on the preconditioner in this sense.
00290   //   In the context of Linear Least Squares problems, preconditioning refers
00291   // to the regularization matrix.  Here the regularization matrix is always a scalar
00292   // multiple of the identity (standard form least squres).
00293   //   The "loss of accuracy" concept is not yet implemented here, becuase it is unclear
00294   // what this means for linear least squares.  LSQR solves an inconsistent system
00295   // in a least-squares sense.  "Loss of accuracy" would correspond to
00296   // the difference between the preconditioned residual and the unpreconditioned residual.
00297   //
00298 
00299   std::cout << " X " << state.sol_norm 
00300             << "  b-AX " << state.resid_norm 
00301             << "  Atr  " << state.mat_resid_norm 
00302             << "  A " << state.frob_mat_norm 
00303             << "  cond  " << state.mat_cond_num
00304             << "  relResNorm " << state.resid_norm/state.bnorm
00305             << "  LS " << state.mat_resid_norm /( state.resid_norm * state.frob_mat_norm )
00306             << std::endl;
00307 
00308   const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00309   const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00310   ScalarType stop_crit_1 = zero; // b = 0, done
00311   if( state.bnorm  > zero )
00312     {
00313       stop_crit_1 = state.resid_norm / state.bnorm;
00314     }
00315   ScalarType stop_crit_2 = zero;
00316   if( state.frob_mat_norm  > zero  && state.resid_norm > zero )
00317     {
00318       stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
00319     }
00320   else
00321     { 
00322      if( state.resid_norm == zero )
00323        {
00324          stop_crit_2 = zero; 
00325        }
00326      else 
00327        {
00328          stop_crit_2 = one; // Initial mat_norm always vanishes
00329        }
00330     }
00331   ScalarType stop_crit_3 = one / state.mat_cond_num;
00332   ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.frob_mat_norm * state.sol_norm / state.bnorm;
00333   ScalarType resid_tol_mach = Teuchos::ScalarTraits< MagnitudeType >::eps() + Teuchos::ScalarTraits< MagnitudeType >::eps() * state.frob_mat_norm * state.sol_norm / state.bnorm;
00334 
00335   // The expected use case for our users is that the linear system will almost
00336   // always be compatible, but occasionally may not be.  However, some users
00337   // may use LSQR for more general cases.  This is why we include the full
00338   // suite of tests, for both compatible and incompatible systems.
00339   //
00340   // Users will have to be educated that sometimes they will get an answer X
00341   // that does _not_ satisfy the linear system AX=B, but _does_ satisfy the
00342   // corresponding least-squares problem.  Perhaps the solution manager should
00343   // provide them with a way to find out.
00344 
00345   // stop_crit_1 is for compatible linear systems.
00346   // stop_crit_2 is for incompatible linear systems.
00347   // stop_crit_3 is for either compatible or incompatible linear systems.
00348 
00349   // Have we met any of the stopping criteria?
00350   if (stop_crit_1 <= resid_tol || stop_crit_2 <= rel_mat_err_ || stop_crit_3 <= rcondMin_ || stop_crit_1 <= resid_tol_mach || stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() || stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps()) {
00351     termIterFlag = true;
00352 
00353     if (stop_crit_1 <= resid_tol )
00354       std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol " << resid_tol << std::endl;
00355 
00356     if (stop_crit_1 <=  resid_tol_mach )
00357       std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol_mach " << resid_tol_mach << std::endl;
00358 
00359     if (stop_crit_2 <= rel_mat_err_ )
00360       std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " rel_mat_err " << rel_mat_err_ << std::endl;
00361 
00362     if (stop_crit_2 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
00363       std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;
00364 
00365     if (stop_crit_3 <= rcondMin_ )
00366       std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " rcondMin_ " << rcondMin_ << std::endl;
00367 
00368     if (stop_crit_3 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
00369       std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;
00370   }
00371 
00372   // update number of consecutive successful iterations
00373   if (!termIterFlag) {
00374     term_iter_ = 0;
00375   } else {
00376     term_iter_++;
00377   }
00378   status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;
00379 
00380   matCondNum_ = state.mat_cond_num; // information that defined convergence
00381   matNorm_ = state.frob_mat_norm;   // in accessible variables 
00382   resNorm_  = state.resid_norm;     
00383   matResNorm_ = state.mat_resid_norm;
00384 
00385   return status_;
00386 }
00387 
00388 template <class ScalarType, class MV, class OP>
00389 void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00390 {
00391   for (int j = 0; j < indent; j++)
00392     os << ' ';
00393   printStatus(os, status_);
00394   os << "limit of condition number = " << condMax_ << std::endl;
00395   os << "limit of condition number = " << condMax_ << std::endl;
00396 }
00397 
00398 template <class ScalarType, class MV, class OP>
00399 void LSQRStatusTest<ScalarType,MV,OP>::printStatus(std::ostream&os, Belos::StatusType type) const
00400 {
00401   os << std::left << std::setw(13) << std::setfill('.');
00402   switch (type) {
00403   case Belos::Passed:
00404     os << "Passed";
00405     break;
00406   case Belos::Failed:
00407     os << "Failed";
00408     break;
00409   case Belos::Undefined:
00410   default:
00411     os << "Undefined";
00412     break;
00413   }
00414   os << std::left << std::setfill(' ');
00415   return;
00416 }
00417 
00418 } // end Belos namespace
00419 
00420 
00421 #endif /* BELOS_LSQR_STATUS_TEST_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines