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