Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_Belos_StatusTest_UnitTests.cpp
Go to the documentation of this file.
00001 #include "Thyra_BelosLinearOpWithSolveFactory.hpp"
00002 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00003 #include "Thyra_MultiVectorStdOps.hpp"
00004 #include "Thyra_VectorBase.hpp"
00005 #include "Thyra_VectorStdOps.hpp"
00006 #include "Thyra_EpetraLinearOp.hpp"
00007 #include "EpetraExt_readEpetraLinearSystem.h"
00008 #include "Epetra_SerialComm.h"
00009 #include "Teuchos_XMLParameterListHelpers.hpp"
00010 #include "Teuchos_toString.hpp"
00011 
00012 #include "Teuchos_UnitTestHarness.hpp"
00013 
00014 
00015 namespace Thyra {
00016 
00017 
00018 //
00019 // Helper code
00020 //
00021 
00022 
00023 const std::string matrixFileName = "nos1.mtx";
00024 
00025 
00026 RCP<const LinearOpBase<double> > getFwdLinearOp()
00027 {
00028   static RCP<const LinearOpBase<double> > fwdLinearOp;
00029   if (is_null(fwdLinearOp)) {
00030     Teuchos::RCP<Epetra_CrsMatrix> epetraCrsMatrix;
00031     EpetraExt::readEpetraLinearSystem( matrixFileName, Epetra_SerialComm(), &epetraCrsMatrix );
00032     fwdLinearOp = epetraLinearOp(epetraCrsMatrix);
00033   }
00034   return fwdLinearOp;
00035 }
00036 
00037 
00039 template<class Scalar>
00040 class MockNormInfReductionFunctional : public ReductionFunctional<Scalar> {
00041 protected:
00042 
00045 
00047   virtual typename ScalarTraits<Scalar>::magnitudeType
00048   reduceImpl(const VectorBase<Scalar> &v) const
00049     { return norm_inf(v); }
00050 
00052   virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const
00053     { return true; }
00054 
00056 
00057 };
00058 
00059 
00064 template<class Scalar>
00065 RCP<MockNormInfReductionFunctional<Scalar> >
00066 createMockNormReductionFunctional()
00067 {
00068   return Teuchos::rcp(new MockNormInfReductionFunctional<Scalar>());
00069 }
00070 
00071 
00074 template<class Scalar>
00075 class MockMaxNormInfEpsReductionFunctional : public ReductionFunctional<Scalar> {
00076 protected:
00077 
00080 
00082   virtual typename ScalarTraits<Scalar>::magnitudeType
00083   reduceImpl(const VectorBase<Scalar> &v) const
00084     {
00085       typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
00086       return std::max(norm_inf(v), ScalarTraits<ScalarMag>::eps());
00087     }
00088 
00090   virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const
00091     { return true; }
00092 
00094 
00095 };
00096 
00097 
00102 template<class Scalar>
00103 RCP<MockMaxNormInfEpsReductionFunctional<Scalar> >
00104 createMockMaxNormInfEpsReductionFunctional()
00105 {
00106   return Teuchos::rcp(new MockMaxNormInfEpsReductionFunctional<Scalar>());
00107 }
00108 
00109 
00110 template<class Scalar>
00111 void runGeneralSolveCriteriaBelosStatusTestCase(
00112   const SolveCriteria<Scalar> &solveCriteria,
00113   const Ptr<RCP<const VectorBase<Scalar> > > &x_out,
00114   const Ptr<RCP<const VectorBase<Scalar> > > &r_out,
00115   bool &success,
00116   FancyOStream &out
00117   )
00118 {
00119 
00120   using Teuchos::describe; using Teuchos::optInArg; using Teuchos::rcpFromRef;
00121   using Teuchos::toString;
00122 
00123   typedef ScalarTraits<Scalar> ST;
00124   typedef typename ST::magnitudeType ScalarMag;
00125 
00126   // A) Set up the linear system
00127   
00128   const RCP<const LinearOpBase<Scalar> > fwdOp = getFwdLinearOp();
00129   out << "\nfwdOp = " << describe(*fwdOp, Teuchos::VERB_MEDIUM) << "\n";
00130   const RCP<VectorBase<Scalar> > b = createMember(fwdOp->range());
00131   V_S(b.ptr(), ST::one());
00132   const RCP<VectorBase<Scalar> > x = createMember(fwdOp->domain());
00133 
00134   // B) Print out the specialized SolveCriteria object
00135 
00136   out << "\nsolveCriteria:\n" << solveCriteria;
00137   
00138   // ToDo: Fill in the rest of the fields!
00139 
00140   // C) Solve the system with the given SolveCriteria object
00141 
00142   const int convergenceTestFrequency = 10;
00143 
00144   const RCP<ParameterList> pl = Teuchos::getParametersFromXmlString(
00145     "<ParameterList name=\"Belos\">"
00146     "  <Parameter name=\"Solver Type\" type=\"string\" value=\"Pseudo Block GMRES\"/>"
00147     "  <Parameter name=\"Convergence Test Frequency\" type=\"int\" value=\""+toString(convergenceTestFrequency)+"\"/>"
00148     "  <ParameterList name=\"Solver Types\">"
00149     "    <ParameterList name=\"Pseudo Block GMRES\">"
00150     "      <Parameter name=\"Block Size\" type=\"int\" value=\"1\"/>"
00151     "      <Parameter name=\"Convergence Tolerance\" type=\"double\" value=\"1e-13\"/>"
00152     "      <Parameter name=\"Output Frequency\" type=\"int\" value=\""+toString(convergenceTestFrequency)+"\"/>"
00153     "      <Parameter name=\"Show Maximum Residual Norm Only\" type=\"bool\" value=\"1\"/>"
00154     "      <Parameter name=\"Maximum Iterations\" type=\"int\" value=\"400\"/>"
00155     "      <Parameter name=\"Verbosity\" type=\"int\" value=\"1\"/>"
00156     "    </ParameterList>"
00157     "  </ParameterList>"
00158     "</ParameterList>"
00159     );
00160   out << "\n\npl:\n" << *pl;
00161 
00162   Thyra::BelosLinearOpWithSolveFactory<Scalar> lowsFactory;
00163   lowsFactory.setParameterList(pl);
00164   lowsFactory.setOStream(rcpFromRef(out));
00165   //lowsFactory.setVerbLevel(Teuchos::VERB_MEDIUM);
00166   lowsFactory.setVerbLevel(Teuchos::VERB_HIGH);
00167 
00168   // NOTE: To get Belos ouptut to be quite, turn down the Belos "Verbosity"
00169   // option above.  To get just the StatusTest VerboseObject output, turn up
00170   // lowsFactory output level.
00171 
00172 
00173   const RCP<LinearOpWithSolveBase<Scalar> > lows = linearOpWithSolve<Scalar>(
00174     lowsFactory, fwdOp);
00175 
00176   V_S(x.ptr(), ST::zero());
00177   SolveStatus<Scalar> solveStatus = solve<Scalar>(*lows, NOTRANS, *b, x.ptr(),
00178     optInArg(solveCriteria));
00179   out << "\nsolveStatus:\n" << solveStatus;
00180 
00181   TEST_COMPARE( solveStatus.achievedTol, <=, solveCriteria.requestedTol );
00182 
00183   // D) Compute the actual residual and return x and r
00184   
00185   const RCP<VectorBase<Scalar> > r = b->clone_v();
00186   fwdOp->apply(NOTRANS, *x, r.ptr(), ST::one(), -ST::one());
00187 
00188   *x_out = x;
00189   *r_out = r;
00190 
00191 }
00192   
00193 
00194 
00195 //
00196 // GeneralSolveCriteriaBelosStatusTest Unit Tests
00197 //
00198 
00199 
00200 TEUCHOS_UNIT_TEST( GeneralSolveCriteriaBelosStatusTest, norm_inf_r_over_norm_inf_r0 )
00201 {
00202   
00203   using Teuchos::outArg;
00204 
00205   typedef double Scalar;
00206   typedef ScalarTraits<Scalar> ST;
00207   typedef ST::magnitudeType ScalarMag;
00208 
00209   SolveCriteria<Scalar> solveCriteria;
00210   solveCriteria.solveMeasureType.numerator = SOLVE_MEASURE_NORM_RESIDUAL;
00211   solveCriteria.numeratorReductionFunc = createMockNormReductionFunctional<Scalar>();
00212   solveCriteria.solveMeasureType.denominator = SOLVE_MEASURE_NORM_SOLUTION;
00213   solveCriteria.denominatorReductionFunc = createMockMaxNormInfEpsReductionFunctional<Scalar>();
00214   solveCriteria.requestedTol = 0.9;
00215 
00216   RCP<const VectorBase<Scalar> > x, r;
00217   runGeneralSolveCriteriaBelosStatusTestCase(solveCriteria, outArg(x), outArg(r),
00218     success, out);
00219 
00220   out << "\nChecking convergence ...\n\n";
00221   
00222   const ScalarMag r_nrm_inf = norm_inf(*r);
00223   const ScalarMag x_nrm_inf = norm_inf(*x);
00224   
00225   out << "||r||inf = " << r_nrm_inf << "\n";
00226   out << "||x||inf = " << x_nrm_inf << "\n";
00227   
00228   TEST_COMPARE( r_nrm_inf / x_nrm_inf, <=, solveCriteria.requestedTol );
00229 
00230 }
00231 
00232 
00233 TEUCHOS_UNIT_TEST( GeneralSolveCriteriaBelosStatusTest, norm_inf_r_over_1 )
00234 {
00235   
00236   using Teuchos::outArg;
00237 
00238   typedef double Scalar;
00239   typedef ScalarTraits<Scalar> ST;
00240   typedef ST::magnitudeType ScalarMag;
00241 
00242   SolveCriteria<Scalar> solveCriteria;
00243   solveCriteria.solveMeasureType.numerator = SOLVE_MEASURE_NORM_RESIDUAL;
00244   solveCriteria.numeratorReductionFunc = createMockNormReductionFunctional<Scalar>();
00245   solveCriteria.solveMeasureType.denominator = SOLVE_MEASURE_ONE;
00246   solveCriteria.requestedTol = 0.9;
00247 
00248   RCP<const VectorBase<Scalar> > x, r;
00249   runGeneralSolveCriteriaBelosStatusTestCase(solveCriteria, outArg(x), outArg(r),
00250     success, out);
00251 
00252   out << "\nChecking convergence ...\n\n";
00253   
00254   const ScalarMag r_nrm_inf = norm_inf(*r);
00255   const ScalarMag x_nrm_inf = norm_inf(*x);
00256   
00257   out << "||r||inf = " << r_nrm_inf << "\n";
00258   out << "||x||inf = " << x_nrm_inf << "\n";
00259   
00260   TEST_COMPARE( r_nrm_inf, <=, solveCriteria.requestedTol );
00261 
00262 }
00263 
00264 
00265 
00266 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines