OptiPack Package Browser (Single Doxygen Collection) Version of the Day
NonlinearCG_UnitTests.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    OptiPack: Collection of simple Thyra-based Optimization ANAs
00006 //                 Copyright (2009) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (rabartl@sandia.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 
00045 #include "Teuchos_UnitTestHarness.hpp"
00046 
00047 #include "Thyra_TestingTools.hpp"
00048 
00049 #include "OptiPack_NonlinearCG.hpp"
00050 #include "GlobiPack_ArmijoPolyInterpLineSearch.hpp"
00051 #include "GlobiPack_BrentsLineSearch.hpp"
00052 #include "Thyra_DiagonalQuadraticResponseOnlyModelEvaluator.hpp"
00053 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00054 #include "Thyra_ModelEvaluatorHelpers.hpp"
00055 #include "Thyra_VectorStdOps.hpp"
00056 #include "Thyra_SpmdVectorSpaceBase.hpp"
00057 #include "RTOpPack_RTOpTHelpers.hpp"
00058 #include "Teuchos_DefaultComm.hpp"
00059 
00060 
00061 namespace {
00062 
00063 
00064 //
00065 // Helper code and declarations
00066 //
00067 
00068 
00069 using Teuchos::as;
00070 using Teuchos::null;
00071 using Teuchos::RCP;
00072 using Teuchos::Ptr;
00073 using Teuchos::rcpFromRef;
00074 using Teuchos::rcp_dynamic_cast;
00075 using Teuchos::ArrayView;
00076 using Teuchos::outArg;
00077 using Teuchos::ParameterList;
00078 using Teuchos::parameterList;
00079 using Teuchos::ScalarTraits;
00080 using Teuchos::FancyOStream;
00081 using Thyra::createMember;
00082 using RTOpPack::ReductTarget;
00083 using RTOpPack::ConstSubVectorView;
00084 using RTOpPack::SubVectorView;
00085 typedef Thyra::Ordinal Ordinal;
00086 using Thyra::VectorSpaceBase;
00087 using Thyra::VectorBase;
00088 using Thyra::applyOp;
00089 using Thyra::DiagonalQuadraticResponseOnlyModelEvaluator;
00090 using Thyra::diagonalQuadraticResponseOnlyModelEvaluator;
00091 using GlobiPack::ArmijoPolyInterpLineSearch;
00092 using GlobiPack::armijoQuadraticLineSearch;
00093 using GlobiPack::BrentsLineSearch;
00094 using GlobiPack::brentsLineSearch;
00095 using OptiPack::NonlinearCG;
00096 using OptiPack::nonlinearCG;
00097 namespace NCGU = OptiPack::NonlinearCGUtils;
00098 
00099 
00100 std::string g_solverType = "FR";
00101 
00102 int g_globalDim = 16;
00103 
00104 double g_solve_tol_scale = 10.0;
00105 
00106 double g_error_tol_scale = 1000.0;
00107 
00108 int g_nonlin_max_iters = 20;
00109 
00110 double g_nonlin_term_factor = 1e-2;
00111 
00112 double g_nonlin_solve_tol = 1e-4;
00113 
00114 double g_nonlin_error_tol = 1e-3;
00115 
00116 
00117 TEUCHOS_STATIC_SETUP()
00118 {
00119   Teuchos::UnitTestRepository::getCLP().setOption(
00120     "solver-type", &g_solverType,
00121     "Type type of nonlinear solver.  Just pass in blah to see valid options." );
00122   Teuchos::UnitTestRepository::getCLP().setOption(
00123     "global-dim", &g_globalDim,
00124     "Number of global vector elements over all processes" );
00125   Teuchos::UnitTestRepository::getCLP().setOption(
00126     "solve-tol-scale", &g_solve_tol_scale,
00127     "Floating point tolerance for nonlinear CG solve for linear CG tests" );
00128   Teuchos::UnitTestRepository::getCLP().setOption(
00129     "error-tol-scale", &g_error_tol_scale,
00130     "Floating point tolerance for error checks for linear CG tests" );
00131   Teuchos::UnitTestRepository::getCLP().setOption(
00132     "nonlin-max-iters", &g_nonlin_max_iters,
00133     "Max nubmer of CG iterations for general nonlinear problem" );
00134   Teuchos::UnitTestRepository::getCLP().setOption(
00135     "nonlin-term-factor", &g_nonlin_term_factor,
00136     "Scale factor for cubic term in objective for general nonlinear problem" );
00137   Teuchos::UnitTestRepository::getCLP().setOption(
00138     "nonlin-solve-tol", &g_nonlin_solve_tol,
00139     "Floating point tolerance for general nonlinear CG solve" );
00140   Teuchos::UnitTestRepository::getCLP().setOption(
00141     "nonlin-error-tol", &g_nonlin_error_tol,
00142     "Floating point tolerance for error checks for general nonlinear CG solve" );
00143 
00144 }
00145 
00146 
00147 template<class Scalar>
00148 const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> >
00149 createModel(
00150   const int globalDim,
00151   const typename ScalarTraits<Scalar>::magnitudeType &g_offset
00152   )
00153 {
00154 
00155   const RCP<const Teuchos::Comm<Thyra::Ordinal> > comm =
00156     Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
00157 
00158   const int numProcs = comm->getSize();
00159   TEUCHOS_TEST_FOR_EXCEPT_MSG( numProcs > globalDim,
00160     "Error, the number of processors can not be greater than the global"
00161     " dimension of the vectors!." );
00162   const int localDim = globalDim / numProcs;
00163   const int localDimRemainder = globalDim % numProcs;
00164   TEUCHOS_TEST_FOR_EXCEPT_MSG( localDimRemainder != 0,
00165     "Error, the number of processors must divide into the global number"
00166     " of elements exactly for now!." );
00167 
00168   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00169     diagonalQuadraticResponseOnlyModelEvaluator<Scalar>(localDim);
00170   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00171   const RCP<VectorBase<Scalar> > ps = createMember(p_space);
00172   const Scalar ps_val = 2.0;
00173   V_S(ps.ptr(), ps_val);
00174   model->setSolutionVector(ps);
00175   model->setScalarOffset(g_offset);
00176 
00177   return model;
00178 
00179 }
00180 
00181 
00182 template<class Scalar>
00183 const RCP<NonlinearCG<Scalar> >
00184 createNonlinearCGSolver(
00185   const RCP<const Thyra::ModelEvaluator<Scalar> > &model,
00186   const RCP<FancyOStream> &out
00187   )
00188 {
00189  
00190   // Set up a quadratic interploation line search that will do just one
00191   // iteration and should exactly minimize a quadratic function.
00192   const RCP<ArmijoPolyInterpLineSearch<Scalar> > linesearch =
00193     armijoQuadraticLineSearch<Scalar>();
00194   const RCP<ParameterList> lsPL = parameterList();
00195   lsPL->set("Armijo Slope Fraction", 0.0);
00196   lsPL->set("Min Backtrack Fraction", 0.0);
00197   lsPL->set("Max Backtrack Fraction", 1e+50);
00198   lsPL->set("Min Num Iterations", 1);
00199   lsPL->set("Max Num Iterations", 2);
00200   linesearch->setParameterList(lsPL);
00201 
00202   const RCP<NonlinearCG<Scalar> > cgSolver =
00203     nonlinearCG<Scalar>(model, 0, 0, linesearch);
00204 
00205   const RCP<ParameterList> pl = parameterList();
00206   pl->set("Solver Type", g_solverType);
00207   cgSolver->setParameterList(pl);
00208 
00209   cgSolver->setOStream(out);
00210 
00211   return cgSolver;
00212 
00213 }
00214 
00215 
00216 //
00217 // RTOp to Assign elements z[i] = i + 1, i = 0...n-1
00218 //
00219 
00220 template<class Scalar>
00221 class TOpAssignValToGlobalIndex : public RTOpPack::RTOpT<Scalar> {
00222 public:
00223   TOpAssignValToGlobalIndex(const Teuchos::Range1D &range = Teuchos::Range1D())
00224     :range_(range)
00225     {}
00226 protected:
00227   bool coord_invariant_impl() const
00228     {
00229       return true;
00230     }
00231   void apply_op_impl(
00232     const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
00233     const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
00234     const Ptr<ReductTarget> &reduct_obj
00235     ) const
00236     {
00237       typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
00238       validate_apply_op( *this, 0, 1, false, sub_vecs, targ_sub_vecs, reduct_obj );
00239       const SubVectorView<Scalar> &z = targ_sub_vecs[0];
00240       const Ordinal z_global_offset = z.globalOffset();
00241       const Ordinal z_sub_dim = z.subDim();
00242       iter_t z_val = z.values().begin();
00243       const ptrdiff_t z_val_s = z.stride();
00244 
00245       for ( int i = 0; i < z_sub_dim; ++i, z_val += z_val_s ) {
00246         const Ordinal global_i = z_global_offset + i;
00247         if (range_.in_range(global_i)) {
00248           *z_val = as<Scalar>(global_i + 1);
00249         }
00250       }
00251     }
00252 private:
00253   Teuchos::Range1D range_;
00254 };
00255 
00256 
00257 
00258 //
00259 // Unit tests
00260 //
00261 
00262 
00263 //
00264 // Check that internal default parameters are set correctly
00265 //
00266 
00267 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, defaultParams, Scalar )
00268 {
00269   typedef ScalarTraits<Scalar> ST;
00270   typedef typename ST::magnitudeType ScalarMag;
00271   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00272   namespace NCGU = OptiPack::NonlinearCGUtils;
00273   const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
00274   TEST_EQUALITY(cgSolver->get_solverType(), as<ScalarMag>(NCGU::solverType_default_integral_val));
00275   TEST_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(NCGU::alpha_init_default));
00276   TEST_EQUALITY(cgSolver->get_alpha_reinit(), NCGU::alpha_reinit_default);
00277   TEST_EQUALITY(cgSolver->get_minIters(), NCGU::minIters_default);
00278   TEST_EQUALITY(cgSolver->get_maxIters(), NCGU::maxIters_default);
00279   TEST_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(NCGU::g_reduct_tol_default));
00280   TEST_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(NCGU::g_grad_tol_default));
00281   TEST_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(NCGU::g_mag_default));
00282 }
00283 
00284 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, defaultParams )
00285 
00286 
00287 //
00288 // Check that internal default parameters are set correctly when gotten off of
00289 // an empty parameter list
00290 //
00291 
00292 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, parseParamsDefaultParams, Scalar )
00293 {
00294   typedef ScalarTraits<Scalar> ST;
00295   typedef typename ST::magnitudeType ScalarMag;
00296   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00297   namespace NCGU = OptiPack::NonlinearCGUtils;
00298   const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
00299   const RCP<ParameterList> pl = parameterList();
00300   cgSolver->setParameterList(pl);
00301   TEST_EQUALITY(cgSolver->get_solverType(), as<ScalarMag>(NCGU::solverType_default_integral_val));
00302   TEST_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(NCGU::alpha_init_default));
00303   TEST_EQUALITY(cgSolver->get_alpha_reinit(), NCGU::alpha_reinit_default);
00304   TEST_EQUALITY(cgSolver->get_minIters(), NCGU::minIters_default);
00305   TEST_EQUALITY(cgSolver->get_maxIters(), NCGU::maxIters_default);
00306   TEST_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(NCGU::g_reduct_tol_default));
00307   TEST_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(NCGU::g_grad_tol_default));
00308   TEST_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(NCGU::g_mag_default));
00309 }
00310 
00311 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, parseParamsDefaultParams )
00312 
00313 
00314 //
00315 // Check that parameter list is parsed correctly
00316 //
00317 
00318 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, parseParams, Scalar )
00319 {
00320   typedef ScalarTraits<Scalar> ST;
00321   typedef typename ST::magnitudeType ScalarMag;
00322   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00323   namespace NCGU = OptiPack::NonlinearCGUtils;
00324   const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
00325   NCGU::ESolverTypes solverType = NCGU::NONLINEAR_CG_PR_PLUS;
00326   const std::string solverTypeStrVal = "PR+";
00327   const double alpha_init = 0.9;
00328   const bool alpha_reinit = true;
00329   const int minIters = 92;
00330   const int maxIters = 99;
00331   const double g_reduct_tol = 2.5;
00332   const double g_grad_tol = 2.8;
00333   const double g_mag = 3.1;
00334   TEST_INEQUALITY( solverType, NCGU::solverType_default_integral_val ); // different!
00335   TEST_INEQUALITY( alpha_reinit, NCGU::alpha_reinit_default ); // different!
00336   const RCP<ParameterList> pl = parameterList();
00337   pl->set("Solver Type", solverTypeStrVal);
00338   pl->set("Initial Linesearch Step Length", alpha_init);
00339   pl->set("Reinitlaize Linesearch Step Length", alpha_reinit);
00340   pl->set("Min Num Iterations", minIters);
00341   pl->set("Max Num Iterations", maxIters);
00342   pl->set("Objective Reduction Tol", g_reduct_tol);
00343   pl->set("Objective Gradient Tol", g_grad_tol);
00344   pl->set("Objective Magnitude", g_mag);
00345   cgSolver->setParameterList(pl);
00346   const ScalarMag tol = SMT::eps();
00347   TEST_EQUALITY(cgSolver->get_solverType(), solverType);
00348   TEST_FLOATING_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(alpha_init), tol);
00349   TEST_EQUALITY(cgSolver->get_alpha_reinit(), alpha_reinit);
00350   TEST_EQUALITY(cgSolver->get_minIters(), minIters);
00351   TEST_EQUALITY(cgSolver->get_maxIters(), maxIters);
00352   TEST_FLOATING_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(g_reduct_tol), tol);
00353   TEST_FLOATING_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(g_grad_tol), tol);
00354   TEST_FLOATING_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(g_mag), tol);
00355 }
00356 
00357 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, parseParams )
00358 
00359 
00360 //
00361 // Print valid parameters
00362 //
00363 
00364 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, printValidParams, Scalar )
00365 {
00366   const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
00367   const RCP<const ParameterList> validPL = cgSolver->getValidParameters();
00368   typedef Teuchos::ParameterList::PrintOptions PO;
00369   out << "\nvalidPL:\n";
00370   validPL->print(out, PO().indent(2).showTypes(true).showFlags(true).showDoc(true));
00371 }
00372 
00373 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, printValidParams )
00374 
00375 
00376 //
00377 // Test basic convergence in one iteration for one eignvalue
00378 //
00379 
00380 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, oneEigenVal, Scalar )
00381 {
00382 
00383   using Teuchos::optInArg;
00384   typedef ScalarTraits<Scalar> ST;
00385   typedef typename ST::magnitudeType ScalarMag;
00386   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00387 
00388   const ScalarMag g_offset = as<ScalarMag>(5.0);
00389   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00390     createModel<Scalar>(g_globalDim, g_offset);
00391   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00392   const int dim = p_space->dim();
00393 
00394   const RCP<NonlinearCG<Scalar> > cgSolver =
00395     createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
00396 
00397   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00398   V_S( p.ptr(), ST::zero() );
00399   
00400   ScalarMag g_opt = -1.0;
00401   const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
00402   const ScalarMag alpha_init = 10.0;
00403   int numIters = -1;
00404     
00405   const NCGU::ESolveReturn solveResult =
00406     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00407       optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
00408 
00409   out << "\n";
00410  
00411   const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
00412   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00413   TEST_EQUALITY( numIters, 1);
00414   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00415   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00416     "p", *p,
00417     "ps", *model->getSolutionVector(),
00418     "err_tol", err_tol,
00419     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00420     &out
00421     );
00422   if (!result) success = false;
00423   
00424 }
00425 
00426 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, oneEigenVal )
00427 
00428 
00429 //
00430 // Test convergence for partially unique eigen values
00431 //
00432 
00433 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, partialEigenVal, Scalar )
00434 {
00435 
00436   using Teuchos::Range1D;
00437   using Teuchos::optInArg;
00438   typedef ScalarTraits<Scalar> ST;
00439   typedef typename ST::magnitudeType ScalarMag;
00440   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00441 
00442   const ScalarMag g_offset = as<ScalarMag>(5.0);
00443   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00444     createModel<Scalar>(g_globalDim, g_offset);
00445 
00446   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00447   const int dim = p_space->dim();
00448 
00449   const int numUniqueEigenVals = 3;
00450   {
00451     const RCP<VectorBase<Scalar> > diag = createMember(p_space);
00452     V_S(diag.ptr(), ST::one());
00453     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(Range1D(0,numUniqueEigenVals-1)),
00454       null, tuple(diag.ptr())(), null );
00455     out << "diag =\n" << *diag;
00456     model->setDiagonalVector(diag);
00457   }
00458 
00459   const RCP<NonlinearCG<Scalar> > cgSolver =
00460     createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
00461 
00462   const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
00463   const int minIters = numUniqueEigenVals;
00464   const int maxIters = minIters + 1;
00465   //pl->set("Min Num Iterations", minIters);
00466   pl->set("Max Num Iterations", maxIters);
00467   cgSolver->setParameterList(pl);
00468 
00469   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00470   V_S( p.ptr(), ST::zero() );
00471   
00472   ScalarMag g_opt = -1.0;
00473   const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
00474   const ScalarMag alpha_init = 10.0;
00475   int numIters = -1;
00476   const NCGU::ESolveReturn solveResult =
00477     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00478       optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
00479 
00480   out << "\n";
00481  
00482   const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
00483   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00484   TEST_COMPARE( numIters, <=, maxIters );
00485   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00486   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00487     "p", *p,
00488     "ps", *model->getSolutionVector(),
00489     "err_tol", err_tol,
00490     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00491     &out
00492     );
00493   if (!result) success = false;
00494   
00495 }
00496 
00497 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, partialEigenVal )
00498 
00499 
00500 //
00501 // Test convergence in full iterations for unique eigen values
00502 //
00503 
00504 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, fullEigenVal, Scalar )
00505 {
00506 
00507   using Teuchos::optInArg;
00508   typedef ScalarTraits<Scalar> ST;
00509   typedef typename ST::magnitudeType ScalarMag;
00510   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00511 
00512   const ScalarMag g_offset = as<ScalarMag>(5.0);
00513   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00514     createModel<Scalar>(g_globalDim, g_offset);
00515   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00516   const int dim = p_space->dim();
00517   {
00518     const RCP<VectorBase<Scalar> > diag = createMember(p_space);
00519     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
00520       null, tuple(diag.ptr())(), null );
00521     out << "diag =\n" << *diag;
00522     model->setDiagonalVector(diag);
00523   }
00524 
00525   const RCP<NonlinearCG<Scalar> > cgSolver =
00526     createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
00527 
00528   const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
00529   const int minIters = dim;
00530   const int maxIters = minIters+2;
00531   //pl->set("Min Num Iterations", minIters);
00532   pl->set("Max Num Iterations", maxIters);
00533   cgSolver->setParameterList(pl);
00534 
00535   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00536   V_S( p.ptr(), ST::zero() );
00537   
00538   ScalarMag g_opt = -1.0;
00539   const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
00540   const ScalarMag alpha_init = 10.0;
00541   int numIters = -1;
00542   const NCGU::ESolveReturn solveResult =
00543     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00544       optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
00545 
00546   out << "\n";
00547  
00548   const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
00549   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00550   TEST_COMPARE( numIters, <=, maxIters );
00551   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00552   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00553     "p", *p,
00554     "ps", *model->getSolutionVector(),
00555     "err_tol", err_tol,
00556     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00557     &out
00558     );
00559   if (!result) success = false;
00560   
00561 }
00562 
00563 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, fullEigenVal )
00564 
00565 
00566 //
00567 // Test convergence for all unique eigen values but using a scalar product
00568 // that has the effect of clustering the eignvalues seen by the nonlinear CG
00569 // ANA.
00570 //
00571 
00572 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, fullEigenValScalarProd, Scalar )
00573 {
00574 
00575   using Teuchos::Range1D;
00576   using Teuchos::optInArg;
00577   typedef ScalarTraits<Scalar> ST;
00578   typedef typename ST::magnitudeType ScalarMag;
00579   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00580 
00581   const ScalarMag g_offset = as<ScalarMag>(5.0);
00582   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00583     createModel<Scalar>(g_globalDim, g_offset);
00584 
00585   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00586   const int dim = p_space->dim();
00587 
00588   {
00589     const RCP<VectorBase<Scalar> > diag = createMember(p_space);
00590     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
00591       null, tuple(diag.ptr())(), null );
00592     out << "diag =\n" << *diag;
00593     model->setDiagonalVector(diag);
00594   }
00595 
00596   const int numUniqueEigenVals = 3;
00597   {
00598     const RCP<VectorBase<Scalar> > diag_bar = createMember(p_space);
00599     V_S(diag_bar.ptr(), ST::one());
00600     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(Range1D(0, numUniqueEigenVals-1)),
00601       null, tuple(diag_bar.ptr())(), null );
00602     //applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
00603     //  null, tuple(diag_bar.ptr())(), null );
00604     out << "diag_bar =\n" << *diag_bar;
00605     model->setDiagonalBarVector(diag_bar);
00606   }
00607 
00608   const RCP<NonlinearCG<Scalar> > cgSolver =
00609     createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
00610 
00611   const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
00612   const int minIters = numUniqueEigenVals;
00613   const int maxIters = minIters + 2;
00614   pl->set("Max Num Iterations", maxIters);
00615   cgSolver->setParameterList(pl);
00616 
00617   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00618   V_S( p.ptr(), ST::zero() );
00619   
00620   ScalarMag g_opt = -1.0;
00621   const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
00622   const ScalarMag alpha_init = 10.0;
00623   int numIters = -1;
00624   const NCGU::ESolveReturn solveResult =
00625     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00626       optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
00627 
00628   out << "\n";
00629  
00630   const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
00631   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00632   TEST_COMPARE( numIters, <=, maxIters );
00633   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00634   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00635     "p", *p,
00636     "ps", *model->getSolutionVector(),
00637     "err_tol", err_tol,
00638     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00639     &out
00640     );
00641   if (!result) success = false;
00642   
00643 }
00644 
00645 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, fullEigenValScalarProd )
00646 
00647 
00648 //
00649 // Test general convergence for a general nonlinear objective
00650 //
00651 
00652 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, generalNonlinearProblem, Scalar )
00653 {
00654 
00655   using Teuchos::optInArg;
00656   typedef ScalarTraits<Scalar> ST;
00657   typedef typename ST::magnitudeType ScalarMag;
00658   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00659 
00660   const ScalarMag g_offset = as<ScalarMag>(5.0);
00661   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00662     createModel<Scalar>(g_globalDim, g_offset);
00663 
00664   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00665 
00666   {
00667     const RCP<VectorBase<Scalar> > diag = createMember(p_space);
00668     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
00669       null, tuple(diag.ptr())(), null );
00670     out << "diag =\n" << *diag;
00671     model->setDiagonalVector(diag);
00672   }
00673 
00674   const ScalarMag nonlinearTermFactor = as<ScalarMag>(g_nonlin_term_factor);
00675   model->setNonlinearTermFactor(nonlinearTermFactor);
00676   
00677   RCP<BrentsLineSearch<Scalar> > linesearch = brentsLineSearch<Scalar>();
00678 
00679   const RCP<NonlinearCG<Scalar> > cgSolver =
00680     nonlinearCG<Scalar>(model, 0, 0, linesearch);
00681 
00682   cgSolver->setOStream(rcpFromRef(out));
00683 
00684   const RCP<ParameterList> pl = parameterList();
00685   pl->set("Solver Type", g_solverType);
00686   pl->set("Reinitlaize Linesearch Step Length", false);
00687   pl->set("Max Num Iterations", g_nonlin_max_iters);
00688   cgSolver->setParameterList(pl);
00689 
00690   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00691   V_S( p.ptr(), ST::zero() );
00692   
00693   ScalarMag g_opt = -1.0;
00694   const ScalarMag tol = as<ScalarMag>(g_nonlin_solve_tol);
00695   const ScalarMag alpha_init = 5.0;
00696   int numIters = -1;
00697   const NCGU::ESolveReturn solveResult =
00698     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00699       optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
00700 
00701   out << "\n";
00702  
00703   const ScalarMag err_tol = as<ScalarMag>(g_nonlin_error_tol);
00704   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00705   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00706   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00707     "p", *p,
00708     "ps", *model->getSolutionVector(),
00709     "err_tol", err_tol,
00710     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00711     &out
00712     );
00713   if (!result) success = false;
00714   
00715 }
00716 
00717 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG,
00718   generalNonlinearProblem )
00719 
00720 
00721 //
00722 // Test general convergence for a general nonlinear objective passing all
00723 // control options through the PL
00724 //
00725 
00726 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, generalNonlinearProblem_PL, Scalar )
00727 {
00728 
00729   using Teuchos::optInArg;
00730   typedef ScalarTraits<Scalar> ST;
00731   typedef typename ST::magnitudeType ScalarMag;
00732   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00733 
00734   const ScalarMag g_offset = as<ScalarMag>(5.0);
00735   const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
00736     createModel<Scalar>(g_globalDim, g_offset);
00737 
00738   const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
00739 
00740   {
00741     const RCP<VectorBase<Scalar> > diag = createMember(p_space);
00742     applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
00743       null, tuple(diag.ptr())(), null );
00744     out << "diag =\n" << *diag;
00745     model->setDiagonalVector(diag);
00746   }
00747 
00748   const ScalarMag nonlinearTermFactor = as<ScalarMag>(g_nonlin_term_factor);
00749   model->setNonlinearTermFactor(nonlinearTermFactor);
00750    
00751   RCP<BrentsLineSearch<Scalar> > linesearch = brentsLineSearch<Scalar>();
00752 
00753   const RCP<NonlinearCG<Scalar> > cgSolver =
00754     nonlinearCG<Scalar>(model, 0, 0, linesearch);
00755 
00756   cgSolver->setOStream(rcpFromRef(out));
00757 
00758   const RCP<VectorBase<Scalar> > p = createMember(p_space);
00759   V_S( p.ptr(), ST::zero() );
00760   
00761   const double tol = as<double>(g_nonlin_solve_tol);
00762   const double alpha_init = as<double>(5.0);
00763   const RCP<ParameterList> pl = parameterList();
00764   pl->set("Solver Type", g_solverType);
00765   pl->set("Initial Linesearch Step Length", alpha_init);
00766   pl->set("Reinitlaize Linesearch Step Length", false);
00767   pl->set("Max Num Iterations", g_nonlin_max_iters);
00768   pl->set("Objective Reduction Tol", tol);
00769   pl->set("Objective Gradient Tol", tol);
00770   cgSolver->setParameterList(pl);
00771 
00772   ScalarMag g_opt = -1.0;
00773   int numIters = -1;
00774   const NCGU::ESolveReturn solveResult =
00775     cgSolver->doSolve( p.ptr(), outArg(g_opt),
00776       null, null, null, outArg(numIters) );
00777   
00778   out << "\n";
00779  
00780   const ScalarMag err_tol = as<ScalarMag>(g_nonlin_error_tol);
00781   TEST_EQUALITY(solveResult, NCGU::SOLVE_SOLUTION_FOUND);
00782   TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
00783   const bool result = Thyra::testRelNormDiffErr<Scalar>(
00784     "p", *p,
00785     "ps", *model->getSolutionVector(),
00786     "err_tol", err_tol,
00787     "2*err_tol", as<ScalarMag>(2.0)*err_tol,
00788     &out
00789     );
00790   if (!result) success = false;
00791   
00792 }
00793 
00794 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG,
00795   generalNonlinearProblem_PL )
00796 
00797 
00798 } // namespace
00799 
00800 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends