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 cond_lim = 0.0, int term_iter_max = 1, MagnitudeType rel_rhs_err = 0.0, MagnitudeType rel_mat_err = 0.0 );
00084 
00086   virtual ~LSQRStatusTest();
00088 
00090 
00091 
00093 
00096   int setCondLim(MagnitudeType cond_lim) {
00097     cond_lim_ = cond_lim; 
00098     cond_tol_ = (cond_lim > 0) ? (Teuchos::ScalarTraits< MagnitudeType >::one() / cond_lim) : Teuchos::ScalarTraits< MagnitudeType >::eps(); 
00099     return(0);}
00100 
00101   int setTermIterMax(int term_iter_max) {
00102     term_iter_max_ = term_iter_max;
00103     if (term_iter_max_ < 1)
00104       term_iter_max_ = 1;
00105     return(0);}
00106 
00107   int setRelRhsErr(MagnitudeType rel_rhs_err) {
00108     rel_rhs_err_ = rel_rhs_err;
00109     return(0);}
00110 
00111   int setRelMatErr(MagnitudeType rel_mat_err) {
00112     rel_mat_err_ = rel_mat_err;
00113     return(0);}
00114 
00116 
00118 
00119 
00121 
00123   Belos::StatusType checkStatus(Belos::Iteration<ScalarType,MV,OP> *iSolver );
00124 
00126   Belos::StatusType getStatus() const {return(status_);}
00127 
00129 
00131 
00132 
00134   void reset();
00135 
00137 
00139 
00140 
00142   void print(std::ostream& os, int indent = 0) const;
00143 
00145   void printStatus(std::ostream& os, Belos::StatusType type) const;
00146 
00148 
00150 
00151 
00153   MagnitudeType getCondLim() const {return(cond_lim_);};
00154 
00156   int getTermIterMax() const {return(term_iter_max_);};
00157 
00159   MagnitudeType getRelRhsErr() const {return(rel_rhs_err_);};
00160 
00162   MagnitudeType getMatErr() const {return(rel_mat_err_);};
00163 
00165 
00168 
00171   Belos::StatusType firstCallCheckStatusSetup(Belos::Iteration<ScalarType,MV,OP>* iSolver);
00173 
00176 
00178   std::string description() const
00179   {
00180     std::ostringstream oss;
00181     oss << "LSQRStatusTest<>: [ limit of condition number = " << cond_lim_ << " ]";
00182     return oss.str();
00183   }
00185 
00186 private:
00187 
00189 
00190 
00192   MagnitudeType cond_lim_;
00193 
00195   int term_iter_max_;
00196 
00198   MagnitudeType rel_rhs_err_;
00199 
00201   MagnitudeType rel_mat_err_;
00202 
00204   MagnitudeType cond_tol_;
00205 
00207   Belos::StatusType status_;
00208 
00210   bool firstcallCheckStatus_;
00211 
00213   int term_iter_;
00214 
00216 
00217 };
00218 
00219 template <class ScalarType, class MV, class OP>
00220 LSQRStatusTest<ScalarType,MV,OP>::LSQRStatusTest( MagnitudeType cond_lim /* = 0 */, int term_iter_max /* = 1 */, MagnitudeType rel_rhs_err /* = 0 */, MagnitudeType rel_mat_err /* = 0 */)
00221   : cond_lim_(cond_lim),
00222     term_iter_max_(term_iter_max),
00223     rel_rhs_err_(rel_rhs_err),
00224     rel_mat_err_(rel_mat_err),
00225     status_(Belos::Undefined),
00226     firstcallCheckStatus_(true)
00227 {}
00228 
00229 template <class ScalarType, class MV, class OP>
00230 LSQRStatusTest<ScalarType,MV,OP>::~LSQRStatusTest()
00231 {}
00232 
00233 template <class ScalarType, class MV, class OP>
00234 void LSQRStatusTest<ScalarType,MV,OP>::reset()
00235 {
00236   status_ = Belos::Undefined;
00237   firstcallCheckStatus_ = true;
00238 }
00239 
00240 template <class ScalarType, class MV, class OP>
00241 Belos::StatusType LSQRStatusTest<ScalarType,MV,OP>::firstCallCheckStatusSetup(Belos::Iteration<ScalarType,MV,OP>* iSolver)
00242 {
00243   if (firstcallCheckStatus_) {
00244       firstcallCheckStatus_ = false;
00245       const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00246       const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00247       term_iter_ = 0;
00248       if (cond_lim_ > MTzero ) 
00249   cond_tol_ = MTone / cond_lim_;
00250       else
00251   cond_tol_ = Teuchos::ScalarTraits< MagnitudeType >::eps();
00252     }
00253     return Belos::Undefined;
00254 }
00255 
00256 template <class ScalarType, class MV, class OP>
00257 Belos::StatusType LSQRStatusTest<ScalarType,MV,OP>::checkStatus( Belos::Iteration<ScalarType,MV,OP>* iSolver) 
00258 {
00259   if (firstcallCheckStatus_) {
00260     Belos::StatusType status = firstCallCheckStatusSetup(iSolver);
00261     if(status==Belos::Failed) {
00262       status_ = Belos::Failed;
00263       return(status_);
00264     }
00265   }
00266 
00267   const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00268   const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00269 
00270   bool flag = false;
00271   LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
00272   LSQRIterationState< ScalarType, MV > state = solver->getState();
00273 
00274   // This section computes the three stopping criteria
00275   ScalarType stop_crit_1 = state.resid_norm / state.bnorm;
00276   ScalarType stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
00277   ScalarType stop_crit_3 = one / state.mat_cond_num;
00278   ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.mat_resid_norm * state.sol_norm / state.bnorm;
00279   ScalarType resid_tol_mach = Teuchos::ScalarTraits< MagnitudeType >::eps() + Teuchos::ScalarTraits< MagnitudeType >::eps() * state.mat_resid_norm * state.sol_norm / state.bnorm;
00280 
00281   // This section checks if any stopping criteria have been met
00282   if (stop_crit_1 <= resid_tol || stop_crit_2 <= rel_mat_err_ || stop_crit_3 <= cond_tol_ || stop_crit_1 <= resid_tol_mach || stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() || stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps()) {
00283     flag = true;
00284   }
00285 
00286 
00287   // history requirement:
00288   // converged if thresholds met for user specified number of consecutive iterations
00289   if (!flag) {
00290     term_iter_ = -1;
00291   }
00292   term_iter_++;
00293   status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;
00294   return status_;
00295 }
00296 
00297 template <class ScalarType, class MV, class OP>
00298 void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00299 {
00300   for (int j = 0; j < indent; j++)
00301     os << ' ';
00302   printStatus(os, status_);
00303   os << "limit of condition number = " << cond_lim_ << std::endl;
00304 }
00305 
00306 template <class ScalarType, class MV, class OP>
00307 void LSQRStatusTest<ScalarType,MV,OP>::printStatus(std::ostream&os, Belos::StatusType type) const
00308 {
00309   os << std::left << std::setw(13) << std::setfill('.');
00310   switch (type) {
00311   case Belos::Passed:
00312     os << "OK";
00313     break;
00314   case Belos::Failed:
00315     os << "Failed";
00316     break;
00317   case Belos::Undefined:
00318   default:
00319     os << "**";
00320     break;
00321   }
00322   os << std::left << std::setfill(' ');
00323   return;
00324 }
00325 
00326 } // end Belos namespace
00327 
00328 
00329 #endif /* BELOS_LSQR_STATUS_TEST_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines