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>::
00241 LSQRStatusTest (MagnitudeType condMax /* = 0 */,
00242                 int term_iter_max /* = 1 */,
00243                 MagnitudeType rel_rhs_err /* = 0 */,
00244                 MagnitudeType rel_mat_err /* = 0 */)
00245   : condMax_(condMax),
00246     term_iter_max_ (term_iter_max),
00247     rel_rhs_err_ (rel_rhs_err),
00248     rel_mat_err_ (rel_mat_err),
00249     rcondMin_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
00250     status_ (Belos::Undefined),
00251     term_iter_ (0),
00252     matCondNum_ ( Teuchos::ScalarTraits<MagnitudeType>::one() ),
00253     matNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
00254     resNorm_  ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
00255     matResNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() )
00256 {}
00257 
00258 template <class ScalarType, class MV, class OP>
00259 LSQRStatusTest<ScalarType,MV,OP>::~LSQRStatusTest()
00260 {}
00261 
00262 template <class ScalarType, class MV, class OP>
00263 void LSQRStatusTest<ScalarType,MV,OP>::reset()
00264 {
00265   status_ = Belos::Undefined;
00266 }
00267 
00268 template <class ScalarType, class MV, class OP>
00269 Belos::StatusType LSQRStatusTest<ScalarType,MV,OP>::checkStatus( Belos::Iteration<ScalarType,MV,OP>* iSolver)
00270 {
00271   const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00272   const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00273   if (condMax_ > MTzero )
00274     {
00275         rcondMin_ = MTone / condMax_;
00276     }
00277   else
00278     {
00279         rcondMin_ = Teuchos::ScalarTraits< MagnitudeType >::eps();
00280     }
00281 
00282   bool termIterFlag = false;
00283   LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
00284   TEUCHOS_ASSERT(solver != NULL);
00285   LSQRIterationState< ScalarType, MV > state = solver->getState();
00286   //
00287   //   LSQR solves a least squares problem.  A converged preconditioned residual norm
00288   // suffices for convergence, but is not necessary.  LSQR sometimes returns a larger
00289   // relative residual norm than what would have been returned by a linear solver.
00290   // This section evaluates three stopping criteria.  In the Solver Manager, this test
00291   // is combined with a generic number of iteration test.
00292   //   If the linear system includes a preconditioner, then the least squares problem
00293   // is solved for the preconditioned linear system.  Preconditioning changes the least
00294   // squares problem (in the sense of changing the norms), and the solution depends
00295   // on the preconditioner in this sense.
00296   //   In the context of Linear Least Squares problems, preconditioning refers
00297   // to the regularization matrix.  Here the regularization matrix is always a scalar
00298   // multiple of the identity (standard form least squres).
00299   //   The "loss of accuracy" concept is not yet implemented here, becuase it is unclear
00300   // what this means for linear least squares.  LSQR solves an inconsistent system
00301   // in a least-squares sense.  "Loss of accuracy" would correspond to
00302   // the difference between the preconditioned residual and the unpreconditioned residual.
00303   //
00304 
00305   std::cout << " X " << state.sol_norm
00306             << "  b-AX " << state.resid_norm
00307             << "  Atr  " << state.mat_resid_norm
00308             << "  A " << state.frob_mat_norm
00309             << "  cond  " << state.mat_cond_num
00310             << "  relResNorm " << state.resid_norm/state.bnorm
00311             << "  LS " << state.mat_resid_norm /( state.resid_norm * state.frob_mat_norm )
00312             << std::endl;
00313 
00314   const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00315   const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00316   ScalarType stop_crit_1 = zero; // b = 0, done
00317   if( state.bnorm  > zero )
00318     {
00319       stop_crit_1 = state.resid_norm / state.bnorm;
00320     }
00321   ScalarType stop_crit_2 = zero;
00322   if( state.frob_mat_norm  > zero  && state.resid_norm > zero )
00323     {
00324       stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
00325     }
00326   else
00327     {
00328      if( state.resid_norm == zero )
00329        {
00330          stop_crit_2 = zero;
00331        }
00332      else
00333        {
00334          stop_crit_2 = one; // Initial mat_norm always vanishes
00335        }
00336     }
00337   ScalarType stop_crit_3 = one / state.mat_cond_num;
00338   ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.frob_mat_norm * state.sol_norm / state.bnorm;
00339   ScalarType resid_tol_mach = Teuchos::ScalarTraits< MagnitudeType >::eps() + Teuchos::ScalarTraits< MagnitudeType >::eps() * state.frob_mat_norm * state.sol_norm / state.bnorm;
00340 
00341   // The expected use case for our users is that the linear system will almost
00342   // always be compatible, but occasionally may not be.  However, some users
00343   // may use LSQR for more general cases.  This is why we include the full
00344   // suite of tests, for both compatible and incompatible systems.
00345   //
00346   // Users will have to be educated that sometimes they will get an answer X
00347   // that does _not_ satisfy the linear system AX=B, but _does_ satisfy the
00348   // corresponding least-squares problem.  Perhaps the solution manager should
00349   // provide them with a way to find out.
00350 
00351   // stop_crit_1 is for compatible linear systems.
00352   // stop_crit_2 is for incompatible linear systems.
00353   // stop_crit_3 is for either compatible or incompatible linear systems.
00354 
00355   // Have we met any of the stopping criteria?
00356   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()) {
00357     termIterFlag = true;
00358 
00359     if (stop_crit_1 <= resid_tol )
00360       std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol " << resid_tol << std::endl;
00361 
00362     if (stop_crit_1 <=  resid_tol_mach )
00363       std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol_mach " << resid_tol_mach << std::endl;
00364 
00365     if (stop_crit_2 <= rel_mat_err_ )
00366       std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " rel_mat_err " << rel_mat_err_ << std::endl;
00367 
00368     if (stop_crit_2 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
00369       std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;
00370 
00371     if (stop_crit_3 <= rcondMin_ )
00372       std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " rcondMin_ " << rcondMin_ << std::endl;
00373 
00374     if (stop_crit_3 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
00375       std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;
00376   }
00377 
00378   // update number of consecutive successful iterations
00379   if (!termIterFlag) {
00380     term_iter_ = 0;
00381   } else {
00382     term_iter_++;
00383   }
00384   status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;
00385 
00386   matCondNum_ = state.mat_cond_num; // information that defined convergence
00387   matNorm_ = state.frob_mat_norm;   // in accessible variables
00388   resNorm_  = state.resid_norm;
00389   matResNorm_ = state.mat_resid_norm;
00390 
00391   return status_;
00392 }
00393 
00394 template <class ScalarType, class MV, class OP>
00395 void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00396 {
00397   for (int j = 0; j < indent; j++)
00398     os << ' ';
00399   printStatus(os, status_);
00400   os << "limit of condition number = " << condMax_ << std::endl;
00401   os << "limit of condition number = " << condMax_ << std::endl;
00402 }
00403 
00404 template <class ScalarType, class MV, class OP>
00405 void LSQRStatusTest<ScalarType,MV,OP>::printStatus(std::ostream&os, Belos::StatusType type) const
00406 {
00407   os << std::left << std::setw(13) << std::setfill('.');
00408   switch (type) {
00409   case Belos::Passed:
00410     os << "Passed";
00411     break;
00412   case Belos::Failed:
00413     os << "Failed";
00414     break;
00415   case Belos::Undefined:
00416   default:
00417     os << "Undefined";
00418     break;
00419   }
00420   os << std::left << std::setfill(' ');
00421   return;
00422 }
00423 
00424 } // end Belos namespace
00425 
00426 
00427 #endif /* BELOS_LSQR_STATUS_TEST_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines