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