Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_GeneralSolveCriteriaBelosStatusTest_def.hpp
Go to the documentation of this file.
00001 
00002 #ifndef THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
00003 #define THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
00004 
00005 #include "Thyra_GeneralSolveCriteriaBelosStatusTest.hpp"
00006 #include "Thyra_VectorSpaceBase.hpp"
00007 #include "Thyra_MultiVectorBase.hpp"
00008 #include "Thyra_VectorBase.hpp"
00009 #include "Thyra_MultiVectorStdOps.hpp"
00010 #include "Thyra_VectorStdOps.hpp"
00011 #include "BelosThyraAdapter.hpp"
00012 #include "Teuchos_DebugDefaultAsserts.hpp"
00013 #include "Teuchos_VerbosityLevel.hpp"
00014 #include "Teuchos_as.hpp"
00015 
00016 
00017 namespace Thyra {
00018 
00019 
00020 // Constructors/initializers/accessors
00021 
00022 
00023 template<class Scalar>
00024 GeneralSolveCriteriaBelosStatusTest<Scalar>::GeneralSolveCriteriaBelosStatusTest()
00025   :convergenceTestFrequency_(-1),
00026    compute_x_(false),
00027    compute_r_(false),
00028    lastCurrIter_(-1),
00029    lastRtnStatus_(Belos::Undefined)
00030 {
00031   GeneralSolveCriteriaBelosStatusTest<Scalar>::reset();
00032 }
00033 
00034 
00035 template<class Scalar>
00036 void GeneralSolveCriteriaBelosStatusTest<Scalar>::setSolveCriteria(
00037   const SolveCriteria<Scalar> &solveCriteria,
00038   const int convergenceTestFrequency
00039   )
00040 {
00041 
00042   using Teuchos::as;
00043   typedef ScalarTraits<ScalarMag> SMT;
00044 
00045   // Make sure the solve criteria is valid
00046 
00047   TEUCHOS_ASSERT_INEQUALITY(solveCriteria.requestedTol, >=, SMT::zero());
00048   TEUCHOS_ASSERT_INEQUALITY(solveCriteria.solveMeasureType.numerator, !=, SOLVE_MEASURE_ONE);
00049   TEUCHOS_ASSERT(nonnull(solveCriteria.numeratorReductionFunc) ||
00050     nonnull(solveCriteria.denominatorReductionFunc) );
00051 
00052   // Remember the solve criteria sett
00053 
00054   solveCriteria_ = solveCriteria;
00055   convergenceTestFrequency_ = convergenceTestFrequency;
00056 
00057   // Determine what needs to be computed on each check
00058 
00059   compute_r_ = solveCriteria.solveMeasureType.contains(SOLVE_MEASURE_NORM_RESIDUAL);
00060 
00061   compute_x_ = (compute_r_ ||
00062     solveCriteria.solveMeasureType.contains(SOLVE_MEASURE_NORM_SOLUTION));
00063 
00064 }
00065 
00066 
00067 template<class Scalar>
00068 ArrayView<const typename ScalarTraits<Scalar>::magnitudeType>
00069 GeneralSolveCriteriaBelosStatusTest<Scalar>::achievedTol() const
00070 {
00071   return lastAchievedTol_;
00072 }
00073 
00074 
00075 // Overridden public functions from Belos::StatusTest
00076 
00077 
00078 template <class Scalar>
00079 Belos::StatusType
00080 GeneralSolveCriteriaBelosStatusTest<Scalar>::checkStatus(
00081   Belos::Iteration<Scalar,MV,OP> *iSolver
00082   )
00083 {
00084 
00085   using Teuchos::null;
00086   typedef ScalarTraits<ScalarMag> SMT;
00087   
00088 #ifdef THYRA_DEBUG
00089   TEUCHOS_ASSERT(iSolver);
00090 #endif
00091 
00092   const int currIter = iSolver->getNumIters();
00093 
00094   if (currIter == 0 || currIter % convergenceTestFrequency_ != 0) {
00095     // We will skip this check!
00096     return Belos::Undefined;
00097   }
00098 
00099   const RCP<FancyOStream> out = this->getOStream();
00100   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00101 
00102   const Belos::LinearProblem<Scalar,MV,OP>& lp = iSolver->getProblem();
00103   const int numRhs = lp.getRHS()->domain()->dim();
00104   
00105   // Compute the rhs norm if requested and not already computed
00106   if (solveCriteria_.solveMeasureType.contains(SOLVE_MEASURE_NORM_RHS)) {
00107     TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||b||");
00108   }
00109 
00110   // Compute the initial residual norm if requested and not already computed
00111   if (solveCriteria_.solveMeasureType.contains(SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
00112     TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||r0||");
00113   }
00114 
00115   // Compute X if requested
00116   RCP<MultiVectorBase<Scalar> > X;
00117   if (compute_x_) {
00118     RCP<MV> X_update = iSolver->getCurrentUpdate();
00119     X = lp.updateSolution(X_update);
00120   }
00121 
00122   // Compute R if requested
00123   RCP<MultiVectorBase<Scalar> > R;
00124   if (compute_r_) {
00125     R = createMembers(lp.getOperator()->range(), X->domain());
00126     lp.computeCurrResVec(&*R, &*X);
00127   }
00128 
00129   // Form numerators and denominaotors gN(vN)/gD(vD) for each system
00130 
00131   lastNumerator_.resize(numRhs);
00132   lastDenominator_.resize(numRhs);
00133 
00134   for (int j = 0; j < numRhs; ++j) {
00135     const RCP<const VectorBase<Scalar> > x_j = (nonnull(X) ? X->col(j) : null);
00136     const RCP<const VectorBase<Scalar> > r_j = (nonnull(R) ? R->col(j) : null);
00137     lastNumerator_[j] = computeReductionFunctional(
00138       solveCriteria_.solveMeasureType.numerator,
00139       solveCriteria_.numeratorReductionFunc.ptr(),
00140       x_j.ptr(), r_j.ptr() );
00141     lastDenominator_[j] = computeReductionFunctional(
00142       solveCriteria_.solveMeasureType.denominator,
00143       solveCriteria_.denominatorReductionFunc.ptr(),
00144       x_j.ptr(), r_j.ptr() );
00145   }
00146 
00147   // Determine if convRatio <= requestedTol
00148 
00149   bool systemsAreConverged = true;
00150   lastAchievedTol_.resize(numRhs);
00151 
00152   for (int j = 0; j < numRhs; ++j) {
00153     const ScalarMag convRatio = lastNumerator_[j] / lastDenominator_[j]; 
00154     lastAchievedTol_[j] = convRatio;
00155     const bool sys_converged_j = (convRatio <= solveCriteria_.requestedTol);
00156     if (includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
00157       printRhsStatus(currIter, j, *out);
00158     }
00159     if (!sys_converged_j) {
00160       systemsAreConverged = false;
00161     }
00162   }
00163   
00164   lastRtnStatus_ = (systemsAreConverged ? Belos::Passed : Belos::Failed);
00165   lastCurrIter_ = currIter;
00166 
00167   return lastRtnStatus_;
00168 
00169 }
00170 
00171 template <class Scalar>
00172 Belos::StatusType
00173 GeneralSolveCriteriaBelosStatusTest<Scalar>::getStatus() const
00174 {
00175   return lastRtnStatus_;
00176 }
00177 
00178 
00179 template <class Scalar>
00180 void GeneralSolveCriteriaBelosStatusTest<Scalar>::reset()
00181 {
00182   r0_nrm_.clear();
00183   b_nrm_.clear();
00184   lastNumerator_.clear();
00185   lastDenominator_.clear();
00186   lastAchievedTol_.clear();
00187   lastCurrIter_ = -1;
00188   lastRtnStatus_ = Belos::Undefined;
00189 }
00190 
00191 
00192 template <class Scalar>
00193 void GeneralSolveCriteriaBelosStatusTest<Scalar>::print(
00194   std::ostream& os, int indent
00195   ) const
00196 {
00197   const int numRhs = lastNumerator_.size();
00198   for (int j = 0; j < numRhs; ++j) {
00199     printRhsStatus(lastCurrIter_, j, os, indent);
00200   }
00201 }
00202 
00203 
00204 // private
00205 
00206 
00207 template <class Scalar>
00208 typename GeneralSolveCriteriaBelosStatusTest<Scalar>::ScalarMag
00209 GeneralSolveCriteriaBelosStatusTest<Scalar>::computeReductionFunctional(
00210   ESolveMeasureNormType measureType,
00211   const Ptr<const ReductionFunctional<Scalar> > &reductFunc,
00212   const Ptr<const VectorBase<Scalar> > &x,
00213   const Ptr<const VectorBase<Scalar> > &r
00214   ) const
00215 {
00216   typedef ScalarTraits<ScalarMag> SMT;
00217   ScalarMag rtn = -SMT::one();
00218   Ptr<const VectorBase<Scalar> > v;
00219   switch(measureType) {
00220     case SOLVE_MEASURE_ONE:
00221       rtn = SMT::one();
00222       break;
00223     case SOLVE_MEASURE_NORM_RESIDUAL:
00224       v = r;
00225       break;
00226     case SOLVE_MEASURE_NORM_SOLUTION:
00227       v = x;
00228       break;
00229     case SOLVE_MEASURE_NORM_INIT_RESIDUAL:
00230       TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||r0||!)")
00231       break;
00232     case SOLVE_MEASURE_NORM_RHS:
00233       TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||b||!)");
00234       break;
00235     TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
00236    }
00237   if (rtn >= SMT::zero()) {
00238     // We already have what we need!
00239   }
00240   else if (nonnull(v) && rtn < SMT::zero()) {
00241     if (nonnull(reductFunc)) {
00242       rtn = reductFunc->reduce(*v);
00243     }
00244     else {
00245       rtn = norm(*v);
00246     }
00247   }
00248   TEUCHOS_IF_ELSE_DEBUG_ASSERT();
00249   return rtn;
00250 }
00251 
00252 
00253 template <class Scalar>
00254 void
00255 GeneralSolveCriteriaBelosStatusTest<Scalar>::printRhsStatus(
00256   const int currIter, const int j, std::ostream &out,
00257   int indent
00258   ) const
00259 {
00260   const ScalarMag convRatio = lastNumerator_[j] / lastDenominator_[j]; 
00261   const bool sys_converged_j = (convRatio <= solveCriteria_.requestedTol);
00262   for (int i = 0; i < indent; ++i) { out << " "; }
00263   out
00264     << "["<<currIter<<"] "
00265     << "gN(vN("<<j<<"))/gD(vD("<<j<<")) = "
00266     << lastNumerator_[j] << "/" << lastDenominator_[j] << " = "
00267     << convRatio << " <= " << solveCriteria_.requestedTol << " : "
00268     << (sys_converged_j ? " true" : "false")
00269     << "\n";
00270 }
00271 
00272 
00273 } // namespace Thyra
00274 
00275 
00276 #endif  // THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines