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