OptiPack  Version of the Day
OptiPack_NonlinearCG_def.hpp
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 #ifndef OPTIPACK_NONLINEAR_CG_DEF_HPP
45 #define OPTIPACK_NONLINEAR_CG_DEF_HPP
46 
47 
48 #include "OptiPack_NonlinearCG_decl.hpp"
49 #include "OptiPack_DefaultPolyLineSearchPointEvaluator.hpp"
50 #include "OptiPack_UnconstrainedOptMeritFunc1D.hpp"
51 #include "Thyra_ModelEvaluatorHelpers.hpp"
52 #include "Thyra_VectorStdOps.hpp"
53 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
54 #include "Teuchos_StandardParameterEntryValidators.hpp"
55 #include "Teuchos_Tuple.hpp"
56 
57 
58 namespace OptiPack {
59 
60 
61 template<typename Scalar>
62 RCP<Teuchos::ParameterEntryValidator>
63 NonlinearCG<Scalar>::solverType_validator_ = Teuchos::null;
64 
65 
66 // Constructor/Initializers/Accessors
67 
68 
69 template<typename Scalar>
71  : paramIndex_(-1),
72  responseIndex_(-1),
73  solverType_(NonlinearCGUtils::solverType_default_integral_val),
74  alpha_init_(NonlinearCGUtils::alpha_init_default),
75  alpha_reinit_(NonlinearCGUtils::alpha_reinit_default),
76  and_conv_tests_(NonlinearCGUtils::and_conv_tests_default),
77  minIters_(NonlinearCGUtils::minIters_default),
78  maxIters_(NonlinearCGUtils::maxIters_default),
79  g_reduct_tol_(NonlinearCGUtils::g_reduct_tol_default),
80  g_grad_tol_(NonlinearCGUtils::g_grad_tol_default),
81  g_mag_(NonlinearCGUtils::g_mag_default),
82  numIters_(0)
83 {}
84 
85 
86 template<typename Scalar>
88  const RCP<const Thyra::ModelEvaluator<Scalar> > &model,
89  const int paramIndex,
90  const int responseIndex,
91  const RCP<GlobiPack::LineSearchBase<Scalar> > &linesearch
92  )
93 {
94  // ToDo: Validate input objects!
95  model_ = model.assert_not_null();
96  paramIndex_ = paramIndex;
97  responseIndex_ = responseIndex;
98  linesearch_ = linesearch.assert_not_null();
99 }
100 
101 
102 template<typename Scalar>
103 NonlinearCGUtils::ESolverTypes NonlinearCG<Scalar>::get_solverType() const
104 {
105  return solverType_;
106 }
107 
108 
109 template<typename Scalar>
112 {
113  return alpha_init_;
114 }
115 
116 
117 template<typename Scalar>
119 {
120  return alpha_reinit_;
121 }
122 
123 
124 template<typename Scalar>
126 {
127  return and_conv_tests_;
128 }
129 
130 
131 template<typename Scalar>
133 {
134  return minIters_;
135 }
136 
137 
138 template<typename Scalar>
140 {
141  return maxIters_;
142 }
143 
144 
145 template<typename Scalar>
148 {
149  return g_reduct_tol_;
150 }
151 
152 
153 template<typename Scalar>
156 {
157  return g_grad_tol_;
158 }
159 
160 
161 template<typename Scalar>
164 {
165  return g_mag_;
166 }
167 
168 
169 // Overridden from ParameterListAcceptor (simple forwarding functions)
170 
171 
172 template<typename Scalar>
173 void NonlinearCG<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
174 {
175  typedef ScalarTraits<Scalar> ST;
176  typedef ScalarTraits<ScalarMag> SMT;
177  namespace NCGU = NonlinearCGUtils;
178  using Teuchos::getParameter;
179  using Teuchos::getIntegralValue;
180  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
181  solverType_ = getIntegralValue<NCGU::ESolverTypes>(*paramList, NCGU::solverType_name);
182  alpha_init_ = getParameter<double>(*paramList, NCGU::alpha_init_name);
183  alpha_reinit_ = getParameter<bool>(*paramList, NCGU::alpha_reinit_name);
184  and_conv_tests_ = getParameter<bool>(*paramList, NCGU::and_conv_tests_name);
185  minIters_ = getParameter<int>(*paramList, NCGU::minIters_name);
186  maxIters_ = getParameter<int>(*paramList, NCGU::maxIters_name);
187  g_reduct_tol_ = getParameter<double>(*paramList, NCGU::g_reduct_tol_name);
188  g_grad_tol_ = getParameter<double>(*paramList, NCGU::g_grad_tol_name);
189  g_mag_ = getParameter<double>(*paramList, NCGU::g_mag_name);
190  TEUCHOS_ASSERT_INEQUALITY( alpha_init_, >, SMT::zero() );
191  TEUCHOS_ASSERT_INEQUALITY( minIters_, >=, 0 );
192  TEUCHOS_ASSERT_INEQUALITY( minIters_, <, maxIters_ );
193  TEUCHOS_ASSERT_INEQUALITY( g_reduct_tol_, >=, SMT::zero() );
194  TEUCHOS_ASSERT_INEQUALITY( g_grad_tol_, >=, SMT::zero() );
195  TEUCHOS_ASSERT_INEQUALITY( g_mag_, >, SMT::zero() );
196  Teuchos::readVerboseObjectSublist(&*paramList, this);
197  setMyParamList(paramList);
198 }
199 
200 
201 template<typename Scalar>
202 RCP<const ParameterList>
204 {
205  using Teuchos::tuple;
206  namespace NCGU = NonlinearCGUtils;
207  static RCP<const ParameterList> validPL;
208  if (is_null(validPL)) {
209  RCP<Teuchos::ParameterList>
210  pl = Teuchos::rcp(new Teuchos::ParameterList());
211  solverType_validator_ =
212  Teuchos::stringToIntegralParameterEntryValidator<NCGU::ESolverTypes>(
213  tuple<std::string>(
214  "FR",
215  "PR+",
216  "FR-PR",
217  "HS"
218  ),
219  tuple<std::string>(
220  "Fletcher-Reeves Method",
221  "Polak-Ribiere Method",
222  "Fletcher-Reeves Polak-Ribiere Hybrid Method",
223  "Hestenes-Stiefel Method"
224  ),
225  tuple<NCGU::ESolverTypes>(
226  NCGU::NONLINEAR_CG_FR,
227  NCGU::NONLINEAR_CG_PR_PLUS,
228  NCGU::NONLINEAR_CG_FR_PR,
229  NCGU::NONLINEAR_CG_HS
230  ),
231  NCGU::solverType_name
232  );
233  pl->set( NCGU::solverType_name, NCGU::solverType_default,
234  "Set the type of nonlinear CG solver algorithm to use.",
235  solverType_validator_ );
236  pl->set( NCGU::alpha_init_name, NCGU::alpha_init_default );
237  pl->set( NCGU::alpha_reinit_name, NCGU::alpha_reinit_default );
238  pl->set( NCGU::and_conv_tests_name, NCGU::and_conv_tests_default );
239  pl->set( NCGU::minIters_name, NCGU::minIters_default );
240  pl->set( NCGU::maxIters_name, NCGU::maxIters_default );
241  pl->set( NCGU::g_reduct_tol_name, NCGU::g_reduct_tol_default );
242  pl->set( NCGU::g_grad_tol_name, NCGU::g_grad_tol_default );
243  pl->set( NCGU::g_mag_name, NCGU::g_mag_default );
244  Teuchos::setupVerboseObjectSublist(&*pl);
245  validPL = pl;
246  // ToDo: Add documentation for these parameters
247  }
248  return validPL;
249 }
250 
251 
252 // Solve
253 
254 
255 template<typename Scalar>
256 NonlinearCGUtils::ESolveReturn
258  const Ptr<Thyra::VectorBase<Scalar> > &p_inout,
259  const Ptr<ScalarMag> &g_opt_out,
260  const Ptr<const ScalarMag> &g_reduct_tol_in,
261  const Ptr<const ScalarMag> &g_grad_tol_in,
262  const Ptr<const ScalarMag> &alpha_init_in,
263  const Ptr<int> &numIters_out
264  )
265 {
266 
267  typedef ScalarTraits<Scalar> ST;
268  typedef ScalarTraits<ScalarMag> SMT;
269 
270  using Teuchos::null;
271  using Teuchos::as;
272  using Teuchos::tuple;
273  using Teuchos::rcpFromPtr;
274  using Teuchos::optInArg;
275  using Teuchos::inOutArg;
276  using GlobiPack::computeValue;
278  using Thyra::VectorSpaceBase;
279  using Thyra::VectorBase;
280  using Thyra::MultiVectorBase;
281  using Thyra::scalarProd;
282  using Thyra::createMember;
283  using Thyra::createMembers;
284  using Thyra::get_ele;
285  using Thyra::norm;
286  using Thyra::V_StV;
287  using Thyra::Vt_S;
288  using Thyra::eval_g_DgDp;
289  typedef Thyra::Ordinal Ordinal;
290  typedef Thyra::ModelEvaluatorBase MEB;
291  namespace NCGU = NonlinearCGUtils;
292  using std::max;
293 
294  // Validate input
295 
296  g_opt_out.assert_not_null();
297 
298  // Set streams
299 
300  const RCP<Teuchos::FancyOStream> out = this->getOStream();
301  linesearch_->setOStream(out);
302 
303  // Determine what step constants will be computed
304 
305  const bool compute_beta_PR =
306  (
307  solverType_ == NCGU::NONLINEAR_CG_PR_PLUS
308  ||
309  solverType_ == NCGU::NONLINEAR_CG_FR_PR
310  );
311 
312  const bool compute_beta_HS = (solverType_ == NCGU::NONLINEAR_CG_HS);
313 
314  //
315  // A) Set up the storage for the algorithm
316  //
317 
318  const RCP<DefaultPolyLineSearchPointEvaluator<Scalar> >
319  pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();
320 
321  const RCP<UnconstrainedOptMeritFunc1D<Scalar> >
322  meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
323  model_, paramIndex_, responseIndex_ );
324 
325  const RCP<const VectorSpaceBase<Scalar> >
326  p_space = model_->get_p_space(paramIndex_),
327  g_space = model_->get_g_space(responseIndex_);
328 
329  // Stoarge for current iteration
330  RCP<VectorBase<Scalar> >
331  p_k = rcpFromPtr(p_inout), // Current solution for p
332  p_kp1 = createMember(p_space), // Trial point for p (in line search)
333  g_vec = createMember(g_space), // Vector (size 1) form of objective g(p)
334  g_grad_k = createMember(p_space), // Gradient of g DgDp^T
335  d_k = createMember(p_space), // Search direction
336  g_grad_k_diff_km1 = null; // g_grad_k - g_grad_km1 (if needed)
337 
338  // Storage for previous iteration
339  RCP<VectorBase<Scalar> >
340  g_grad_km1 = null, // Will allocate if we need it!
341  d_km1 = null; // Will allocate if we need it!
342  ScalarMag
343  alpha_km1 = SMT::zero(),
344  g_km1 = SMT::zero(),
345  g_grad_km1_inner_g_grad_km1 = SMT::zero(),
346  g_grad_km1_inner_d_km1 = SMT::zero();
347 
348  if (compute_beta_PR || compute_beta_HS) {
349  g_grad_km1 = createMember(p_space);
350  g_grad_k_diff_km1 = createMember(p_space);
351  }
352 
353  if (compute_beta_HS) {
354  d_km1 = createMember(p_space);
355  }
356 
357  //
358  // B) Do the nonlinear CG iterations
359  //
360 
361  *out << "\nStarting nonlinear CG iterations ...\n";
362 
363  if (and_conv_tests_) {
364  *out << "\nNOTE: Using AND of convergence tests!\n";
365  }
366  else {
367  *out << "\nNOTE: Using OR of convergence tests!\n";
368  }
369 
370  const Scalar alpha_init =
371  ( !is_null(alpha_init_in) ? *alpha_init_in : alpha_init_ );
372  const Scalar g_reduct_tol =
373  ( !is_null(g_reduct_tol_in) ? *g_reduct_tol_in : g_reduct_tol_ );
374  const Scalar g_grad_tol =
375  ( !is_null(g_grad_tol_in) ? *g_grad_tol_in : g_grad_tol_ );
376 
377  const Ordinal globalDim = p_space->dim();
378 
379  bool foundSolution = false;
380  bool fatalLinesearchFailure = false;
381  bool restart = true;
382  int numConsecutiveLineSearchFailures = 0;
383 
384  int numConsecutiveIters = 0;
385 
386  for (numIters_ = 0; numIters_ < maxIters_; ++numIters_, ++numConsecutiveIters) {
387 
388  Teuchos::OSTab tab(out);
389 
390  *out << "\nNonlinear CG Iteration k = " << numIters_ << "\n";
391 
392  Teuchos::OSTab tab2(out);
393 
394  //
395  // B.1) Evaluate the point (on first iteration)
396  //
397 
398  eval_g_DgDp(
399  *model_, paramIndex_, *p_k, responseIndex_,
400  numIters_ == 0 ? g_vec.ptr() : null, // Only on first iteration
401  MEB::Derivative<Scalar>(g_grad_k, MEB::DERIV_MV_GRADIENT_FORM) );
402 
403  const ScalarMag g_k = get_ele(*g_vec, 0);
404  // Above: If numIters_ > 0, then g_vec was updated in meritFunc->eval(...).
405 
406  //
407  // B.2) Check for convergence
408  //
409 
410  // B.2.a) ||g_k - g_km1|| |g_k + g_mag| <= g_reduct_tol
411 
412  bool g_reduct_converged = false;
413 
414  if (numIters_ > 0) {
415 
416  const ScalarMag g_reduct = g_k - g_km1;
417 
418  *out << "\ng_k - g_km1 = "<<g_reduct<<"\n";
419 
420  const ScalarMag g_reduct_err =
421  SMT::magnitude(g_reduct / SMT::magnitude(g_k + g_mag_));
422 
423  g_reduct_converged = (g_reduct_err <= g_reduct_tol);
424 
425  *out << "\nCheck convergence: |g_k - g_km1| / |g_k + g_mag| = "<<g_reduct_err
426  << (g_reduct_converged ? " <= " : " > ")
427  << "g_reduct_tol = "<<g_reduct_tol<<"\n";
428 
429  }
430 
431  // B.2.b) ||g_grad_k|| g_mag <= g_grad_tol
432 
433  const Scalar g_grad_k_inner_g_grad_k = scalarProd<Scalar>(*g_grad_k, *g_grad_k);
434  const ScalarMag norm_g_grad_k = ST::magnitude(ST::squareroot(g_grad_k_inner_g_grad_k));
435 
436  *out << "\n||g_grad_k|| = "<<norm_g_grad_k << "\n";
437 
438  const ScalarMag g_grad_err = norm_g_grad_k / g_mag_;
439 
440  const bool g_grad_converged = (g_grad_err <= g_grad_tol);
441 
442  *out << "\nCheck convergence: ||g_grad_k|| / g_mag = "<<g_grad_err
443  << (g_grad_converged ? " <= " : " > ")
444  << "g_grad_tol = "<<g_grad_tol<<"\n";
445 
446  // B.2.c) Convergence status
447 
448  bool isConverged = false;
449  if (and_conv_tests_) {
450  isConverged = g_reduct_converged && g_grad_converged;
451  }
452  else {
453  isConverged = g_reduct_converged || g_grad_converged;
454  }
455 
456  if (isConverged) {
457  if (numIters_ < minIters_) {
458  *out << "\nnumIters="<<numIters_<<" < minIters="<<minIters_
459  << ", continuing on!\n";
460  }
461  else {
462  *out << "\nFound solution, existing algorithm!\n";
463  foundSolution = true;
464  }
465  }
466  else {
467  *out << "\nNot converged!\n";
468  }
469 
470  if (foundSolution) {
471  break;
472  }
473 
474  //
475  // B.3) Compute the search direction d_k
476  //
477 
478  if (numConsecutiveIters == globalDim) {
479 
480  *out << "\nThe number of consecutive iterations exceeds the"
481  << " global dimension so restarting!\n";
482 
483  restart = true;
484 
485  }
486 
487  if (restart) {
488 
489  *out << "\nResetting search direction back to steppest descent!\n";
490 
491  // d_k = -g_grad_k
492  V_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
493 
494  restart = false;
495 
496  }
497  else {
498 
499  // g_grad_k - g_grad_km1
500  if (!is_null(g_grad_k_diff_km1)) {
501  V_VmV( g_grad_k_diff_km1.ptr(), *g_grad_k, *g_grad_km1 );
502  }
503 
504  // beta_FR = inner(g_grad_k, g_grad_k) / inner(g_grad_km1, g_grad_km1)
505  const Scalar beta_FR =
506  g_grad_k_inner_g_grad_k / g_grad_km1_inner_g_grad_km1;
507  *out << "\nbeta_FR = " << beta_FR << "\n";
508  // NOTE: Computing beta_FR is free so we might as well just do it!
509 
510  // beta_PR = inner(g_grad_k, g_grad_k - g_grad_km1) /
511  // inner(g_grad_km1, g_grad_km1)
512  Scalar beta_PR = ST::zero();
513  if (compute_beta_PR) {
514  beta_PR =
515  inner(*g_grad_k, *g_grad_k_diff_km1) / g_grad_km1_inner_g_grad_km1;
516  *out << "\nbeta_PR = " << beta_PR << "\n";
517  }
518 
519  // beta_HS = inner(g_grad_k, g_grad_k - g_grad_km1) /
520  // inner(g_grad_k - g_grad_km1, d_km1)
521  Scalar beta_HS = ST::zero();
522  if (compute_beta_HS) {
523  beta_HS =
524  inner(*g_grad_k, *g_grad_k_diff_km1) / inner(*g_grad_k_diff_km1, *d_km1);
525  *out << "\nbeta_HS = " << beta_HS << "\n";
526  }
527 
528  Scalar beta_k = ST::zero();
529  switch(solverType_) {
530  case NCGU::NONLINEAR_CG_FR: {
531  beta_k = beta_FR;
532  break;
533  }
534  case NCGU::NONLINEAR_CG_PR_PLUS: {
535  beta_k = max(beta_PR, ST::zero());
536  break;
537  }
538  case NCGU::NONLINEAR_CG_FR_PR: {
539  // NOTE: This does not seem to be working :-(
540  if (numConsecutiveIters < 2) {
541  beta_k = beta_PR;
542  }
543  else if (beta_PR < -beta_FR)
544  beta_k = -beta_FR;
545  else if (ST::magnitude(beta_PR) <= beta_FR)
546  beta_k = beta_PR;
547  else // beta_PR > beta_FR
548  beta_k = beta_FR;
549  break;
550  }
551  case NCGU::NONLINEAR_CG_HS: {
552  beta_k = beta_HS;
553  break;
554  }
555  default:
556  TEUCHOS_TEST_FOR_EXCEPT(true);
557  }
558  *out << "\nbeta_k = " << beta_k << "\n";
559 
560  // d_k = beta_k * d_last + -g_grad_k
561  if (!is_null(d_km1))
562  V_StV( d_k.ptr(), beta_k, *d_km1 );
563  else
564  Vt_S( d_k.ptr(), beta_k );
565  Vp_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
566 
567  }
568 
569  //
570  // B.4) Perform the line search
571  //
572 
573  // B.4.a) Compute the initial step length
574 
575  Scalar alpha_k = as<Scalar>(-1.0);
576 
577  if (numIters_ == 0) {
578  alpha_k = alpha_init;
579  }
580  else {
581  if (alpha_reinit_) {
582  alpha_k = alpha_init;
583  }
584  else {
585  alpha_k = alpha_km1;
586  // ToDo: Implement better logic from Nocedal and Wright for selecting
587  // this step length after first iteration!
588  }
589  }
590 
591  // B.4.b) Perform the linesearch (computing updated quantities in process)
592 
593  pointEvaluator->initialize(tuple<RCP<const VectorBase<Scalar> > >(p_k, d_k)());
594 
595  ScalarMag g_grad_k_inner_d_k = ST::zero();
596 
597  // Set up the merit function to only compute the value
598  meritFunc->setEvaluationQuantities(pointEvaluator, p_kp1, g_vec, null);
599 
600  PointEval1D<ScalarMag> point_k(ST::zero(), g_k);
601  if (linesearch_->requiresBaseDeriv()) {
602  g_grad_k_inner_d_k = scalarProd(*g_grad_k, *d_k);
603  point_k.Dphi = g_grad_k_inner_d_k;
604  }
605 
606  ScalarMag g_kp1 = computeValue(*meritFunc, alpha_k);
607  // NOTE: The above call updates p_kp1 and g_vec as well!
608 
609  PointEval1D<ScalarMag> point_kp1(alpha_k, g_kp1);
610 
611  const bool linesearchResult = linesearch_->doLineSearch(
612  *meritFunc, point_k, inOutArg(point_kp1), null );
613 
614  alpha_k = point_kp1.alpha;
615  g_kp1 = point_kp1.phi;
616 
617  if (linesearchResult) {
618  numConsecutiveLineSearchFailures = 0;
619  }
620  else {
621  if (numConsecutiveLineSearchFailures==0) {
622  *out << "\nLine search failure, resetting the search direction!\n";
623  restart = true;
624  }
625  if (numConsecutiveLineSearchFailures==1) {
626  *out << "\nLine search failure on last iteration also, terminating algorithm!\n";
627  fatalLinesearchFailure = true;
628  }
629  ++numConsecutiveLineSearchFailures;
630  }
631 
632  if (fatalLinesearchFailure) {
633  break;
634  }
635 
636  //
637  // B.5) Transition to the next iteration
638  //
639 
640  alpha_km1 = alpha_k;
641  g_km1 = g_k;
642  g_grad_km1_inner_g_grad_km1 = g_grad_k_inner_g_grad_k;
643  g_grad_km1_inner_d_km1 = g_grad_k_inner_d_k;
644  std::swap(p_k, p_kp1);
645  if (!is_null(g_grad_km1))
646  std::swap(g_grad_km1, g_grad_k);
647  if (!is_null(d_km1))
648  std::swap(d_k, d_km1);
649 
650 #ifdef TEUCHOS_DEBUG
651  // Make sure we compute these correctly before they are used!
652  V_S(g_grad_k.ptr(), ST::nan());
653  V_S(p_kp1.ptr(), ST::nan());
654 #endif
655 
656  }
657 
658  //
659  // C) Final clean up
660  //
661 
662  // Get the most current value of g(p)
663  *g_opt_out = get_ele(*g_vec, 0);
664 
665  // Make sure that the final value for p has been copied in!
666  V_V( p_inout, *p_k );
667 
668  if (!is_null(numIters_out)) {
669  *numIters_out = numIters_;
670  }
671 
672  if (numIters_ == maxIters_) {
673  *out << "\nMax nonlinear CG iterations exceeded!\n";
674  }
675 
676  if (foundSolution) {
677  return NonlinearCGUtils::SOLVE_SOLUTION_FOUND;
678  }
679  else if(fatalLinesearchFailure) {
680  return NonlinearCGUtils::SOLVE_LINSEARCH_FAILURE;
681  }
682 
683  // Else, the max number of iterations was exceeded
684  return NonlinearCGUtils::SOLVE_MAX_ITERS_EXCEEDED;
685 
686 }
687 
688 
689 } // namespace OptiPack
690 
691 
692 #endif // OPTIPACK_NONLINEAR_CG_DEF_HPP
void setParameterList(RCP< ParameterList > const &paramList)
RCP< const ParameterList > getValidParameters() const
NonlinearCG()
Construct with default parameters.
void initialize(const RCP< const Thyra::ModelEvaluator< Scalar > > &model, const int paramIndex, const int responseIndex, const RCP< GlobiPack::LineSearchBase< Scalar > > &linesearch)
Initialize.
NonlinearCGUtils::ESolverTypes get_solverType() const
NonlinearCGUtils::ESolveReturn doSolve(const Ptr< Thyra::VectorBase< Scalar > > &p, const Ptr< ScalarMag > &g_opt, const Ptr< const ScalarMag > &g_reduct_tol=Teuchos::null, const Ptr< const ScalarMag > &g_grad_tol=Teuchos::null, const Ptr< const ScalarMag > &alpha_init=Teuchos::null, const Ptr< int > &numIters=Teuchos::null)
Perform a solve.
ScalarTraits< Scalar >::magnitudeType ScalarMag