OptiPack Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
NonlinearCG_UnitTests.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // OptiPack: Collection of simple Thyra-based Optimization ANAs
6 // Copyright (2009) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 
45 #include "Teuchos_UnitTestHarness.hpp"
46 
47 #include "Thyra_TestingTools.hpp"
48 
49 #include "OptiPack_NonlinearCG.hpp"
50 #include "GlobiPack_ArmijoPolyInterpLineSearch.hpp"
51 #include "GlobiPack_BrentsLineSearch.hpp"
52 #include "Thyra_DiagonalQuadraticResponseOnlyModelEvaluator.hpp"
53 #include "Thyra_DefaultSpmdVectorSpace.hpp"
54 #include "Thyra_ModelEvaluatorHelpers.hpp"
55 #include "Thyra_VectorStdOps.hpp"
56 #include "Thyra_SpmdVectorSpaceBase.hpp"
57 #include "RTOpPack_RTOpTHelpers.hpp"
58 #include "Teuchos_DefaultComm.hpp"
59 
60 
61 namespace {
62 
63 
64 //
65 // Helper code and declarations
66 //
67 
68 
69 using Teuchos::as;
70 using Teuchos::null;
71 using Teuchos::RCP;
72 using Teuchos::Ptr;
73 using Teuchos::rcpFromRef;
74 using Teuchos::rcp_dynamic_cast;
75 using Teuchos::ArrayView;
76 using Teuchos::outArg;
78 using Teuchos::parameterList;
81 using Thyra::createMember;
85 typedef Thyra::Ordinal Ordinal;
86 using Thyra::VectorSpaceBase;
87 using Thyra::VectorBase;
88 using Thyra::applyOp;
89 using Thyra::DiagonalQuadraticResponseOnlyModelEvaluator;
90 using Thyra::diagonalQuadraticResponseOnlyModelEvaluator;
92 using GlobiPack::armijoQuadraticLineSearch;
94 using GlobiPack::brentsLineSearch;
96 using OptiPack::nonlinearCG;
97 namespace NCGU = OptiPack::NonlinearCGUtils;
98 
99 
100 std::string g_solverType = "FR";
101 
102 int g_globalDim = 16;
103 
104 double g_solve_tol_scale = 10.0;
105 
106 double g_error_tol_scale = 1000.0;
107 
108 int g_nonlin_max_iters = 20;
109 
110 double g_nonlin_term_factor = 1e-2;
111 
112 double g_nonlin_solve_tol = 1e-4;
113 
114 double g_nonlin_error_tol = 1e-3;
115 
116 
118 {
120  "solver-type", &g_solverType,
121  "Type type of nonlinear solver. Just pass in blah to see valid options." );
123  "global-dim", &g_globalDim,
124  "Number of global vector elements over all processes" );
126  "solve-tol-scale", &g_solve_tol_scale,
127  "Floating point tolerance for nonlinear CG solve for linear CG tests" );
129  "error-tol-scale", &g_error_tol_scale,
130  "Floating point tolerance for error checks for linear CG tests" );
132  "nonlin-max-iters", &g_nonlin_max_iters,
133  "Max nubmer of CG iterations for general nonlinear problem" );
135  "nonlin-term-factor", &g_nonlin_term_factor,
136  "Scale factor for cubic term in objective for general nonlinear problem" );
138  "nonlin-solve-tol", &g_nonlin_solve_tol,
139  "Floating point tolerance for general nonlinear CG solve" );
141  "nonlin-error-tol", &g_nonlin_error_tol,
142  "Floating point tolerance for error checks for general nonlinear CG solve" );
143 
144 }
145 
146 
147 template<class Scalar>
148 const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> >
149 createModel(
150  const int globalDim,
151  const typename ScalarTraits<Scalar>::magnitudeType &g_offset
152  )
153 {
154 
155  const RCP<const Teuchos::Comm<Thyra::Ordinal> > comm =
157 
158  const int numProcs = comm->getSize();
159  TEUCHOS_TEST_FOR_EXCEPT_MSG( numProcs > globalDim,
160  "Error, the number of processors can not be greater than the global"
161  " dimension of the vectors!." );
162  const int localDim = globalDim / numProcs;
163  const int localDimRemainder = globalDim % numProcs;
164  TEUCHOS_TEST_FOR_EXCEPT_MSG( localDimRemainder != 0,
165  "Error, the number of processors must divide into the global number"
166  " of elements exactly for now!." );
167 
168  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
169  diagonalQuadraticResponseOnlyModelEvaluator<Scalar>(localDim);
170  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
171  const RCP<VectorBase<Scalar> > ps = createMember(p_space);
172  const Scalar ps_val = 2.0;
173  V_S(ps.ptr(), ps_val);
174  model->setSolutionVector(ps);
175  model->setScalarOffset(g_offset);
176 
177  return model;
178 
179 }
180 
181 
182 template<class Scalar>
183 const RCP<NonlinearCG<Scalar> >
184 createNonlinearCGSolver(
185  const RCP<const Thyra::ModelEvaluator<Scalar> > &model,
186  const RCP<FancyOStream> &out
187  )
188 {
189 
190  // Set up a quadratic interploation line search that will do just one
191  // iteration and should exactly minimize a quadratic function.
192  const RCP<ArmijoPolyInterpLineSearch<Scalar> > linesearch =
193  armijoQuadraticLineSearch<Scalar>();
194  const RCP<ParameterList> lsPL = parameterList();
195  lsPL->set("Armijo Slope Fraction", 0.0);
196  lsPL->set("Min Backtrack Fraction", 0.0);
197  lsPL->set("Max Backtrack Fraction", 1e+50);
198  lsPL->set("Min Num Iterations", 1);
199  lsPL->set("Max Num Iterations", 2);
200  linesearch->setParameterList(lsPL);
201 
202  const RCP<NonlinearCG<Scalar> > cgSolver =
203  nonlinearCG<Scalar>(model, 0, 0, linesearch);
204 
205  const RCP<ParameterList> pl = parameterList();
206  pl->set("Solver Type", g_solverType);
207  cgSolver->setParameterList(pl);
208 
209  cgSolver->setOStream(out);
210 
211  return cgSolver;
212 
213 }
214 
215 
216 //
217 // RTOp to Assign elements z[i] = i + 1, i = 0...n-1
218 //
219 
220 template<class Scalar>
221 class TOpAssignValToGlobalIndex : public RTOpPack::RTOpT<Scalar> {
222 public:
223  TOpAssignValToGlobalIndex (const Teuchos::Range1D& theRange = Teuchos::Range1D ())
224  : range_ (theRange)
225  {}
226 protected:
227  bool coord_invariant_impl() const
228  {
229  return true;
230  }
231  void apply_op_impl(
232  const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
233  const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
234  const Ptr<ReductTarget> &reduct_obj
235  ) const
236  {
237  typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
238  validate_apply_op( *this, 0, 1, false, sub_vecs, targ_sub_vecs, reduct_obj );
239  const SubVectorView<Scalar> &z = targ_sub_vecs[0];
240  const Ordinal z_global_offset = z.globalOffset();
241  const Ordinal z_sub_dim = z.subDim();
242  iter_t z_val = z.values().begin();
243  const ptrdiff_t z_val_s = z.stride();
244 
245  for ( int i = 0; i < z_sub_dim; ++i, z_val += z_val_s ) {
246  const Ordinal global_i = z_global_offset + i;
247  if (range_.in_range(global_i)) {
248  *z_val = as<Scalar>(global_i + 1);
249  }
250  }
251  }
252 private:
253  Teuchos::Range1D range_;
254 };
255 
256 
257 
258 //
259 // Unit tests
260 //
261 
262 
263 //
264 // Check that internal default parameters are set correctly
265 //
266 
267 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, defaultParams, Scalar )
268 {
269  typedef ScalarTraits<Scalar> ST;
270  typedef typename ST::magnitudeType ScalarMag;
272  namespace NCGU = OptiPack::NonlinearCGUtils;
273  const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
274  TEST_EQUALITY(cgSolver->get_solverType(), as<ScalarMag>(NCGU::solverType_default_integral_val));
275  TEST_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(NCGU::alpha_init_default));
276  TEST_EQUALITY(cgSolver->get_alpha_reinit(), NCGU::alpha_reinit_default);
277  TEST_EQUALITY(cgSolver->get_minIters(), NCGU::minIters_default);
278  TEST_EQUALITY(cgSolver->get_maxIters(), NCGU::maxIters_default);
279  TEST_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(NCGU::g_reduct_tol_default));
280  TEST_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(NCGU::g_grad_tol_default));
281  TEST_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(NCGU::g_mag_default));
282 }
283 
284 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, defaultParams )
285 
286 
287 //
288 // Check that internal default parameters are set correctly when gotten off of
289 // an empty parameter list
290 //
291 
292 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, parseParamsDefaultParams, Scalar )
293 {
294  typedef ScalarTraits<Scalar> ST;
295  typedef typename ST::magnitudeType ScalarMag;
297  namespace NCGU = OptiPack::NonlinearCGUtils;
298  const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
299  const RCP<ParameterList> pl = parameterList();
300  cgSolver->setParameterList(pl);
301  TEST_EQUALITY(cgSolver->get_solverType(), as<ScalarMag>(NCGU::solverType_default_integral_val));
302  TEST_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(NCGU::alpha_init_default));
303  TEST_EQUALITY(cgSolver->get_alpha_reinit(), NCGU::alpha_reinit_default);
304  TEST_EQUALITY(cgSolver->get_minIters(), NCGU::minIters_default);
305  TEST_EQUALITY(cgSolver->get_maxIters(), NCGU::maxIters_default);
306  TEST_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(NCGU::g_reduct_tol_default));
307  TEST_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(NCGU::g_grad_tol_default));
308  TEST_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(NCGU::g_mag_default));
309 }
310 
311 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, parseParamsDefaultParams )
312 
313 
314 //
315 // Check that parameter list is parsed correctly
316 //
317 
318 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, parseParams, Scalar )
319 {
320  typedef ScalarTraits<Scalar> ST;
321  typedef typename ST::magnitudeType ScalarMag;
323  namespace NCGU = OptiPack::NonlinearCGUtils;
324  const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
326  const std::string solverTypeStrVal = "PR+";
327  const double alpha_init = 0.9;
328  const bool alpha_reinit = true;
329  const int minIters = 92;
330  const int maxIters = 99;
331  const double g_reduct_tol = 2.5;
332  const double g_grad_tol = 2.8;
333  const double g_mag = 3.1;
334  TEST_INEQUALITY( solverType, NCGU::solverType_default_integral_val ); // different!
335  TEST_INEQUALITY( alpha_reinit, NCGU::alpha_reinit_default ); // different!
336  const RCP<ParameterList> pl = parameterList();
337  pl->set("Solver Type", solverTypeStrVal);
338  pl->set("Initial Linesearch Step Length", alpha_init);
339  pl->set("Reinitlaize Linesearch Step Length", alpha_reinit);
340  pl->set("Min Num Iterations", minIters);
341  pl->set("Max Num Iterations", maxIters);
342  pl->set("Objective Reduction Tol", g_reduct_tol);
343  pl->set("Objective Gradient Tol", g_grad_tol);
344  pl->set("Objective Magnitude", g_mag);
345  cgSolver->setParameterList(pl);
346  const ScalarMag tol = SMT::eps();
347  TEST_EQUALITY(cgSolver->get_solverType(), solverType);
348  TEST_FLOATING_EQUALITY(cgSolver->get_alpha_init(), as<ScalarMag>(alpha_init), tol);
349  TEST_EQUALITY(cgSolver->get_alpha_reinit(), alpha_reinit);
350  TEST_EQUALITY(cgSolver->get_minIters(), minIters);
351  TEST_EQUALITY(cgSolver->get_maxIters(), maxIters);
352  TEST_FLOATING_EQUALITY(cgSolver->get_g_reduct_tol(), as<ScalarMag>(g_reduct_tol), tol);
353  TEST_FLOATING_EQUALITY(cgSolver->get_g_grad_tol(), as<ScalarMag>(g_grad_tol), tol);
354  TEST_FLOATING_EQUALITY(cgSolver->get_g_mag(), as<ScalarMag>(g_mag), tol);
355 }
356 
358 
359 
360 //
361 // Print valid parameters
362 //
363 
364 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, printValidParams, Scalar )
365 {
366  const RCP<NonlinearCG<Scalar> > cgSolver = nonlinearCG<Scalar>();
367  const RCP<const ParameterList> validPL = cgSolver->getValidParameters();
369  out << "\nvalidPL:\n";
370  validPL->print(out, PO().indent(2).showTypes(true).showFlags(true).showDoc(true));
371 }
372 
373 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, printValidParams )
374 
375 
376 //
377 // Test basic convergence in one iteration for one eignvalue
378 //
379 
380 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, oneEigenVal, Scalar )
381 {
382 
383  using Teuchos::optInArg;
384  typedef ScalarTraits<Scalar> ST;
385  typedef typename ST::magnitudeType ScalarMag;
387 
388  const ScalarMag g_offset = as<ScalarMag>(5.0);
389  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
390  createModel<Scalar>(g_globalDim, g_offset);
391  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
392  const int dim = p_space->dim();
393 
394  const RCP<NonlinearCG<Scalar> > cgSolver =
395  createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
396 
397  const RCP<VectorBase<Scalar> > p = createMember(p_space);
398  V_S( p.ptr(), ST::zero() );
399 
400  ScalarMag g_opt = -1.0;
401  const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
402  const ScalarMag alpha_init = 10.0;
403  int numIters = -1;
404 
405  const NCGU::ESolveReturn solveResult =
406  cgSolver->doSolve( p.ptr(), outArg(g_opt),
407  optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
408 
409  out << "\n";
410 
411  const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
413  TEST_EQUALITY( numIters, 1);
414  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
415  const bool result = Thyra::testRelNormDiffErr<Scalar>(
416  "p", *p,
417  "ps", *model->getSolutionVector(),
418  "err_tol", err_tol,
419  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
420  &out
421  );
422  if (!result) success = false;
423 
424 }
425 
427 
428 
429 //
430 // Test convergence for partially unique eigen values
431 //
432 
433 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, partialEigenVal, Scalar )
434 {
435 
436  using Teuchos::Range1D;
437  using Teuchos::optInArg;
438  typedef ScalarTraits<Scalar> ST;
439  typedef typename ST::magnitudeType ScalarMag;
441 
442  const ScalarMag g_offset = as<ScalarMag>(5.0);
443  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
444  createModel<Scalar>(g_globalDim, g_offset);
445 
446  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
447  const int dim = p_space->dim();
448 
449  const int numUniqueEigenVals = 3;
450  {
451  const RCP<VectorBase<Scalar> > diag = createMember(p_space);
452  V_S(diag.ptr(), ST::one());
453  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(Range1D(0,numUniqueEigenVals-1)),
454  null, tuple(diag.ptr())(), null );
455  out << "diag =\n" << *diag;
456  model->setDiagonalVector(diag);
457  }
458 
459  const RCP<NonlinearCG<Scalar> > cgSolver =
460  createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
461 
462  const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
463  const int minIters = numUniqueEigenVals;
464  const int maxIters = minIters + 1;
465  //pl->set("Min Num Iterations", minIters);
466  pl->set("Max Num Iterations", maxIters);
467  cgSolver->setParameterList(pl);
468 
469  const RCP<VectorBase<Scalar> > p = createMember(p_space);
470  V_S( p.ptr(), ST::zero() );
471 
472  ScalarMag g_opt = -1.0;
473  const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
474  const ScalarMag alpha_init = 10.0;
475  int numIters = -1;
476  const NCGU::ESolveReturn solveResult =
477  cgSolver->doSolve( p.ptr(), outArg(g_opt),
478  optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
479 
480  out << "\n";
481 
482  const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
484  TEST_COMPARE( numIters, <=, maxIters );
485  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
486  const bool result = Thyra::testRelNormDiffErr<Scalar>(
487  "p", *p,
488  "ps", *model->getSolutionVector(),
489  "err_tol", err_tol,
490  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
491  &out
492  );
493  if (!result) success = false;
494 
495 }
496 
497 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, partialEigenVal )
498 
499 
500 //
501 // Test convergence in full iterations for unique eigen values
502 //
503 
504 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, fullEigenVal, Scalar )
505 {
506 
507  using Teuchos::optInArg;
508  typedef ScalarTraits<Scalar> ST;
509  typedef typename ST::magnitudeType ScalarMag;
511 
512  const ScalarMag g_offset = as<ScalarMag>(5.0);
513  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
514  createModel<Scalar>(g_globalDim, g_offset);
515  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
516  const int dim = p_space->dim();
517  {
518  const RCP<VectorBase<Scalar> > diag = createMember(p_space);
519  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
520  null, tuple(diag.ptr())(), null );
521  out << "diag =\n" << *diag;
522  model->setDiagonalVector(diag);
523  }
524 
525  const RCP<NonlinearCG<Scalar> > cgSolver =
526  createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
527 
528  const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
529  const int minIters = dim;
530  const int maxIters = minIters+2;
531  //pl->set("Min Num Iterations", minIters);
532  pl->set("Max Num Iterations", maxIters);
533  cgSolver->setParameterList(pl);
534 
535  const RCP<VectorBase<Scalar> > p = createMember(p_space);
536  V_S( p.ptr(), ST::zero() );
537 
538  ScalarMag g_opt = -1.0;
539  const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
540  const ScalarMag alpha_init = 10.0;
541  int numIters = -1;
542  const NCGU::ESolveReturn solveResult =
543  cgSolver->doSolve( p.ptr(), outArg(g_opt),
544  optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
545 
546  out << "\n";
547 
548  const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
550  TEST_COMPARE( numIters, <=, maxIters );
551  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
552  const bool result = Thyra::testRelNormDiffErr<Scalar>(
553  "p", *p,
554  "ps", *model->getSolutionVector(),
555  "err_tol", err_tol,
556  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
557  &out
558  );
559  if (!result) success = false;
560 
561 }
562 
564 
565 
566 //
567 // Test convergence for all unique eigen values but using a scalar product
568 // that has the effect of clustering the eignvalues seen by the nonlinear CG
569 // ANA.
570 //
571 
572 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, fullEigenValScalarProd, Scalar )
573 {
574 
575  using Teuchos::Range1D;
576  using Teuchos::optInArg;
577  typedef ScalarTraits<Scalar> ST;
578  typedef typename ST::magnitudeType ScalarMag;
580 
581  const ScalarMag g_offset = as<ScalarMag>(5.0);
582  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
583  createModel<Scalar>(g_globalDim, g_offset);
584 
585  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
586  const int dim = p_space->dim();
587 
588  {
589  const RCP<VectorBase<Scalar> > diag = createMember(p_space);
590  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
591  null, tuple(diag.ptr())(), null );
592  out << "diag =\n" << *diag;
593  model->setDiagonalVector(diag);
594  }
595 
596  const int numUniqueEigenVals = 3;
597  {
598  const RCP<VectorBase<Scalar> > diag_bar = createMember(p_space);
599  V_S(diag_bar.ptr(), ST::one());
600  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(Range1D(0, numUniqueEigenVals-1)),
601  null, tuple(diag_bar.ptr())(), null );
602  //applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
603  // null, tuple(diag_bar.ptr())(), null );
604  out << "diag_bar =\n" << *diag_bar;
605  model->setDiagonalBarVector(diag_bar);
606  }
607 
608  const RCP<NonlinearCG<Scalar> > cgSolver =
609  createNonlinearCGSolver<Scalar>(model, rcpFromRef(out));
610 
611  const RCP<ParameterList> pl = cgSolver->getNonconstParameterList();
612  const int minIters = numUniqueEigenVals;
613  const int maxIters = minIters + 2;
614  pl->set("Max Num Iterations", maxIters);
615  cgSolver->setParameterList(pl);
616 
617  const RCP<VectorBase<Scalar> > p = createMember(p_space);
618  V_S( p.ptr(), ST::zero() );
619 
620  ScalarMag g_opt = -1.0;
621  const ScalarMag tol = as<Scalar>(g_solve_tol_scale * dim) * ST::eps();
622  const ScalarMag alpha_init = 10.0;
623  int numIters = -1;
624  const NCGU::ESolveReturn solveResult =
625  cgSolver->doSolve( p.ptr(), outArg(g_opt),
626  optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
627 
628  out << "\n";
629 
630  const ScalarMag err_tol = as<Scalar>(g_error_tol_scale * dim) * ST::eps();
632  TEST_COMPARE( numIters, <=, maxIters );
633  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
634  const bool result = Thyra::testRelNormDiffErr<Scalar>(
635  "p", *p,
636  "ps", *model->getSolutionVector(),
637  "err_tol", err_tol,
638  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
639  &out
640  );
641  if (!result) success = false;
642 
643 }
644 
645 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( NonlinearCG, fullEigenValScalarProd )
646 
647 
648 //
649 // Test general convergence for a general nonlinear objective
650 //
651 
652 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, generalNonlinearProblem, Scalar )
653 {
654 
655  using Teuchos::optInArg;
656  typedef ScalarTraits<Scalar> ST;
657  typedef typename ST::magnitudeType ScalarMag;
659 
660  const ScalarMag g_offset = as<ScalarMag>(5.0);
661  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
662  createModel<Scalar>(g_globalDim, g_offset);
663 
664  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
665 
666  {
667  const RCP<VectorBase<Scalar> > diag = createMember(p_space);
668  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
669  null, tuple(diag.ptr())(), null );
670  out << "diag =\n" << *diag;
671  model->setDiagonalVector(diag);
672  }
673 
674  const ScalarMag nonlinearTermFactor = as<ScalarMag>(g_nonlin_term_factor);
675  model->setNonlinearTermFactor(nonlinearTermFactor);
676 
677  RCP<BrentsLineSearch<Scalar> > linesearch = brentsLineSearch<Scalar>();
678 
679  const RCP<NonlinearCG<Scalar> > cgSolver =
680  nonlinearCG<Scalar>(model, 0, 0, linesearch);
681 
682  cgSolver->setOStream(rcpFromRef(out));
683 
684  const RCP<ParameterList> pl = parameterList();
685  pl->set("Solver Type", g_solverType);
686  pl->set("Reinitlaize Linesearch Step Length", false);
687  pl->set("Max Num Iterations", g_nonlin_max_iters);
688  cgSolver->setParameterList(pl);
689 
690  const RCP<VectorBase<Scalar> > p = createMember(p_space);
691  V_S( p.ptr(), ST::zero() );
692 
693  ScalarMag g_opt = -1.0;
694  const ScalarMag tol = as<ScalarMag>(g_nonlin_solve_tol);
695  const ScalarMag alpha_init = 5.0;
696  int numIters = -1;
697  const NCGU::ESolveReturn solveResult =
698  cgSolver->doSolve( p.ptr(), outArg(g_opt),
699  optInArg(tol), optInArg(tol), optInArg(alpha_init), outArg(numIters) );
700 
701  out << "\n";
702 
703  const ScalarMag err_tol = as<ScalarMag>(g_nonlin_error_tol);
705  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
706  const bool result = Thyra::testRelNormDiffErr<Scalar>(
707  "p", *p,
708  "ps", *model->getSolutionVector(),
709  "err_tol", err_tol,
710  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
711  &out
712  );
713  if (!result) success = false;
714 
715 }
716 
718  generalNonlinearProblem )
719 
720 
721 //
722 // Test general convergence for a general nonlinear objective passing all
723 // control options through the PL
724 //
725 
726 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( NonlinearCG, generalNonlinearProblem_PL, Scalar )
727 {
728 
729  using Teuchos::optInArg;
730  typedef ScalarTraits<Scalar> ST;
731  typedef typename ST::magnitudeType ScalarMag;
733 
734  const ScalarMag g_offset = as<ScalarMag>(5.0);
735  const RCP<DiagonalQuadraticResponseOnlyModelEvaluator<Scalar> > model =
736  createModel<Scalar>(g_globalDim, g_offset);
737 
738  const RCP<const VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
739 
740  {
741  const RCP<VectorBase<Scalar> > diag = createMember(p_space);
742  applyOp<Scalar>( TOpAssignValToGlobalIndex<Scalar>(),
743  null, tuple(diag.ptr())(), null );
744  out << "diag =\n" << *diag;
745  model->setDiagonalVector(diag);
746  }
747 
748  const ScalarMag nonlinearTermFactor = as<ScalarMag>(g_nonlin_term_factor);
749  model->setNonlinearTermFactor(nonlinearTermFactor);
750 
751  RCP<BrentsLineSearch<Scalar> > linesearch = brentsLineSearch<Scalar>();
752 
753  const RCP<NonlinearCG<Scalar> > cgSolver =
754  nonlinearCG<Scalar>(model, 0, 0, linesearch);
755 
756  cgSolver->setOStream(rcpFromRef(out));
757 
758  const RCP<VectorBase<Scalar> > p = createMember(p_space);
759  V_S( p.ptr(), ST::zero() );
760 
761  const double tol = as<double>(g_nonlin_solve_tol);
762  const double alpha_init = as<double>(5.0);
763  const RCP<ParameterList> pl = parameterList();
764  pl->set("Solver Type", g_solverType);
765  pl->set("Initial Linesearch Step Length", alpha_init);
766  pl->set("Reinitlaize Linesearch Step Length", false);
767  pl->set("Max Num Iterations", g_nonlin_max_iters);
768  pl->set("Objective Reduction Tol", tol);
769  pl->set("Objective Gradient Tol", tol);
770  cgSolver->setParameterList(pl);
771 
772  ScalarMag g_opt = -1.0;
773  int numIters = -1;
774  const NCGU::ESolveReturn solveResult =
775  cgSolver->doSolve( p.ptr(), outArg(g_opt),
776  null, null, null, outArg(numIters) );
777 
778  out << "\n";
779 
780  const ScalarMag err_tol = as<ScalarMag>(g_nonlin_error_tol);
782  TEST_FLOATING_EQUALITY(g_opt, g_offset, err_tol);
783  const bool result = Thyra::testRelNormDiffErr<Scalar>(
784  "p", *p,
785  "ps", *model->getSolutionVector(),
786  "err_tol", err_tol,
787  "2*err_tol", as<ScalarMag>(2.0)*err_tol,
788  &out
789  );
790  if (!result) success = false;
791 
792 }
793 
795  generalNonlinearProblem_PL )
796 
797 
798 } // namespace
799 
800 
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Teuchos_Conditions, NumberConditionSerialization, T)
#define TEST_INEQUALITY(v1, v2)
static CommandLineProcessor & getCLP()
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Concrete class implementing several nonlinear CG algorithms.
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
virtual void apply_op_impl(const ArrayView< const ConstSubVectorView< Scalar > > &sub_vecs, const ArrayView< const SubVectorView< Scalar > > &targ_sub_vecs, const Ptr< ReductTarget > &reduct_obj) const =0
virtual bool coord_invariant_impl() const
const ESolverTypes solverType_default_integral_val
TypeTo as(const TypeFrom &t)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES(TEST_GROUP, TEST_NAME)
void validate_apply_op(const RTOpT< Scalar > &op, const int allowed_num_sub_vecs, const int allowed_num_targ_sub_vecs, const bool expect_reduct_obj, const ArrayView< const ConstSubVectorView< Scalar > > &sub_vecs, const ArrayView< const SubVectorView< Scalar > > &targ_sub_vecs, const Ptr< const ReductTarget > &reduct_obj)
#define TEST_EQUALITY(v1, v2)
Teuchos_Ordinal Ordinal
TEUCHOS_STATIC_SETUP()
#define TEST_COMPARE(v1, comp, v2)