44 #ifndef OPTIPACK_NONLINEAR_CG_DEF_HPP 45 #define OPTIPACK_NONLINEAR_CG_DEF_HPP 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" 61 template<
typename Scalar>
63 NonlinearCG<Scalar>::solverType_validator_ = Teuchos::null;
69 template<
typename Scalar>
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),
86 template<
typename Scalar>
90 const int responseIndex,
95 model_ = model.assert_not_null();
96 paramIndex_ = paramIndex;
97 responseIndex_ = responseIndex;
98 linesearch_ = linesearch.assert_not_null();
102 template<
typename Scalar>
109 template<
typename Scalar>
117 template<
typename Scalar>
120 return alpha_reinit_;
124 template<
typename Scalar>
127 return and_conv_tests_;
131 template<
typename Scalar>
138 template<
typename Scalar>
145 template<
typename Scalar>
149 return g_reduct_tol_;
153 template<
typename Scalar>
161 template<
typename Scalar>
172 template<
typename Scalar>
177 namespace NCGU = NonlinearCGUtils;
178 using Teuchos::getParameter;
179 using Teuchos::getIntegralValue;
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);
196 Teuchos::readVerboseObjectSublist(&*paramList,
this);
197 setMyParamList(paramList);
201 template<
typename Scalar>
205 using Teuchos::tuple;
206 namespace NCGU = NonlinearCGUtils;
211 solverType_validator_ =
212 Teuchos::stringToIntegralParameterEntryValidator<NCGU::ESolverTypes>(
220 "Fletcher-Reeves Method",
221 "Polak-Ribiere Method",
222 "Fletcher-Reeves Polak-Ribiere Hybrid Method",
223 "Hestenes-Stiefel Method" 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
231 NCGU::solverType_name
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);
255 template<
typename Scalar>
256 NonlinearCGUtils::ESolveReturn
258 const Ptr<Thyra::VectorBase<Scalar> > &p_inout,
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;
288 using Thyra::eval_g_DgDp;
289 typedef Thyra::Ordinal Ordinal;
290 typedef Thyra::ModelEvaluatorBase MEB;
291 namespace NCGU = NonlinearCGUtils;
301 linesearch_->setOStream(out);
305 const bool compute_beta_PR =
307 solverType_ == NCGU::NONLINEAR_CG_PR_PLUS
309 solverType_ == NCGU::NONLINEAR_CG_FR_PR
312 const bool compute_beta_HS = (solverType_ == NCGU::NONLINEAR_CG_HS);
319 pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();
322 meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
323 model_, paramIndex_, responseIndex_ );
326 p_space = model_->get_p_space(paramIndex_),
327 g_space = model_->get_g_space(responseIndex_);
331 p_k = rcpFromPtr(p_inout),
332 p_kp1 = createMember(p_space),
333 g_vec = createMember(g_space),
334 g_grad_k = createMember(p_space),
335 d_k = createMember(p_space),
336 g_grad_k_diff_km1 = null;
343 alpha_km1 = SMT::zero(),
345 g_grad_km1_inner_g_grad_km1 = SMT::zero(),
346 g_grad_km1_inner_d_km1 = SMT::zero();
348 if (compute_beta_PR || compute_beta_HS) {
349 g_grad_km1 = createMember(p_space);
350 g_grad_k_diff_km1 = createMember(p_space);
353 if (compute_beta_HS) {
354 d_km1 = createMember(p_space);
361 *out <<
"\nStarting nonlinear CG iterations ...\n";
363 if (and_conv_tests_) {
364 *out <<
"\nNOTE: Using AND of convergence tests!\n";
367 *out <<
"\nNOTE: Using OR of convergence tests!\n";
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_ );
377 const Ordinal globalDim = p_space->dim();
379 bool foundSolution =
false;
380 bool fatalLinesearchFailure =
false;
382 int numConsecutiveLineSearchFailures = 0;
384 int numConsecutiveIters = 0;
386 for (numIters_ = 0; numIters_ < maxIters_; ++numIters_, ++numConsecutiveIters) {
390 *out <<
"\nNonlinear CG Iteration k = " << numIters_ <<
"\n";
399 *model_, paramIndex_, *p_k, responseIndex_,
400 numIters_ == 0 ? g_vec.
ptr() : null,
401 MEB::Derivative<Scalar>(g_grad_k, MEB::DERIV_MV_GRADIENT_FORM) );
403 const ScalarMag g_k = get_ele(*g_vec, 0);
412 bool g_reduct_converged =
false;
418 *out <<
"\ng_k - g_km1 = "<<g_reduct<<
"\n";
421 SMT::magnitude(g_reduct / SMT::magnitude(g_k + g_mag_));
423 g_reduct_converged = (g_reduct_err <= g_reduct_tol);
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";
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));
436 *out <<
"\n||g_grad_k|| = "<<norm_g_grad_k <<
"\n";
438 const ScalarMag g_grad_err = norm_g_grad_k / g_mag_;
440 const bool g_grad_converged = (g_grad_err <= g_grad_tol);
442 *out <<
"\nCheck convergence: ||g_grad_k|| / g_mag = "<<g_grad_err
443 << (g_grad_converged ?
" <= " :
" > ")
444 <<
"g_grad_tol = "<<g_grad_tol<<
"\n";
448 bool isConverged =
false;
449 if (and_conv_tests_) {
450 isConverged = g_reduct_converged && g_grad_converged;
453 isConverged = g_reduct_converged || g_grad_converged;
457 if (numIters_ < minIters_) {
458 *out <<
"\nnumIters="<<numIters_<<
" < minIters="<<minIters_
459 <<
", continuing on!\n";
462 *out <<
"\nFound solution, existing algorithm!\n";
463 foundSolution =
true;
467 *out <<
"\nNot converged!\n";
478 if (numConsecutiveIters == globalDim) {
480 *out <<
"\nThe number of consecutive iterations exceeds the" 481 <<
" global dimension so restarting!\n";
489 *out <<
"\nResetting search direction back to steppest descent!\n";
492 V_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
500 if (!
is_null(g_grad_k_diff_km1)) {
501 V_VmV( g_grad_k_diff_km1.ptr(), *g_grad_k, *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";
512 Scalar beta_PR = ST::zero();
513 if (compute_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";
521 Scalar beta_HS = ST::zero();
522 if (compute_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";
528 Scalar beta_k = ST::zero();
529 switch(solverType_) {
530 case NCGU::NONLINEAR_CG_FR: {
534 case NCGU::NONLINEAR_CG_PR_PLUS: {
535 beta_k = max(beta_PR, ST::zero());
538 case NCGU::NONLINEAR_CG_FR_PR: {
540 if (numConsecutiveIters < 2) {
543 else if (beta_PR < -beta_FR)
545 else if (ST::magnitude(beta_PR) <= beta_FR)
551 case NCGU::NONLINEAR_CG_HS: {
558 *out <<
"\nbeta_k = " << beta_k <<
"\n";
562 V_StV( d_k.ptr(), beta_k, *d_km1 );
564 Vt_S( d_k.ptr(), beta_k );
565 Vp_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
575 Scalar alpha_k = as<Scalar>(-1.0);
577 if (numIters_ == 0) {
578 alpha_k = alpha_init;
582 alpha_k = alpha_init;
593 pointEvaluator->initialize(tuple<
RCP<
const VectorBase<Scalar> > >(p_k, d_k)());
595 ScalarMag g_grad_k_inner_d_k = ST::zero();
598 meritFunc->setEvaluationQuantities(pointEvaluator, p_kp1, g_vec, null);
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;
606 ScalarMag g_kp1 = computeValue(*meritFunc, alpha_k);
609 PointEval1D<ScalarMag> point_kp1(alpha_k, g_kp1);
611 const bool linesearchResult = linesearch_->doLineSearch(
612 *meritFunc, point_k, inOutArg(point_kp1), null );
614 alpha_k = point_kp1.alpha;
615 g_kp1 = point_kp1.phi;
617 if (linesearchResult) {
618 numConsecutiveLineSearchFailures = 0;
621 if (numConsecutiveLineSearchFailures==0) {
622 *out <<
"\nLine search failure, resetting the search direction!\n";
625 if (numConsecutiveLineSearchFailures==1) {
626 *out <<
"\nLine search failure on last iteration also, terminating algorithm!\n";
627 fatalLinesearchFailure =
true;
629 ++numConsecutiveLineSearchFailures;
632 if (fatalLinesearchFailure) {
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);
646 std::swap(g_grad_km1, g_grad_k);
648 std::swap(d_k, d_km1);
652 V_S(g_grad_k.ptr(), ST::nan());
653 V_S(p_kp1.ptr(), ST::nan());
663 *g_opt_out = get_ele(*g_vec, 0);
666 V_V( p_inout, *p_k );
669 *numIters_out = numIters_;
672 if (numIters_ == maxIters_) {
673 *out <<
"\nMax nonlinear CG iterations exceeded!\n";
677 return NonlinearCGUtils::SOLVE_SOLUTION_FOUND;
679 else if(fatalLinesearchFailure) {
680 return NonlinearCGUtils::SOLVE_LINSEARCH_FAILURE;
684 return NonlinearCGUtils::SOLVE_MAX_ITERS_EXCEEDED;
692 #endif // OPTIPACK_NONLINEAR_CG_DEF_HPP bool get_and_conv_tests() const
ScalarMag get_alpha_init() const
bool is_null(const boost::shared_ptr< T > &p)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
bool get_alpha_reinit() const
const Ptr< T > & assert_not_null() const
void setParameterList(RCP< ParameterList > const ¶mList)
RCP< const ParameterList > getValidParameters() const
NonlinearCG()
Construct with default parameters.
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
ScalarMag get_g_reduct_tol() const
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
TypeTo as(const TypeFrom &t)
ScalarMag get_g_grad_tol() const
ScalarMag get_g_mag() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
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