44 #ifndef OPTIPACK_NONLINEAR_CG_DEF_HPP 45 #define OPTIPACK_NONLINEAR_CG_DEF_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>
62 RCP<Teuchos::ParameterEntryValidator>
69 template<
typename Scalar>
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>
176 typedef ScalarTraits<ScalarMag> SMT;
177 namespace NCGU = NonlinearCGUtils;
178 using Teuchos::getParameter;
179 using Teuchos::getIntegralValue;
180 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
185 minIters_ = getParameter<int>(*paramList, NCGU::minIters_name);
186 maxIters_ = getParameter<int>(*paramList, NCGU::maxIters_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);
201 template<
typename Scalar>
202 RCP<const ParameterList>
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>(
220 "Fletcher-Reeves Method",
221 "Polak-Ribiere Method",
222 "Fletcher-Reeves Polak-Ribiere Hybrid Method",
223 "Hestenes-Stiefel Method" 225 tuple<NCGU::ESolverTypes>(
234 "Set the type of nonlinear CG solver algorithm to use.",
235 solverType_validator_ );
239 pl->set( NCGU::minIters_name, NCGU::minIters_default );
240 pl->set( NCGU::maxIters_name, NCGU::maxIters_default );
244 Teuchos::setupVerboseObjectSublist(&*pl);
255 template<
typename Scalar>
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
267 typedef ScalarTraits<Scalar> ST;
268 typedef ScalarTraits<ScalarMag> SMT;
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;
296 g_opt_out.assert_not_null();
300 const RCP<Teuchos::FancyOStream> out = this->getOStream();
301 linesearch_->setOStream(out);
305 const bool compute_beta_PR =
318 const RCP<DefaultPolyLineSearchPointEvaluator<Scalar> >
319 pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();
321 const RCP<UnconstrainedOptMeritFunc1D<Scalar> >
322 meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
323 model_, paramIndex_, responseIndex_ );
325 const RCP<const VectorSpaceBase<Scalar> >
326 p_space = model_->get_p_space(paramIndex_),
327 g_space = model_->get_g_space(responseIndex_);
330 RCP<VectorBase<Scalar> >
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;
339 RCP<VectorBase<Scalar> >
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) {
388 Teuchos::OSTab tab(out);
390 *out <<
"\nNonlinear CG Iteration k = " << numIters_ <<
"\n";
392 Teuchos::OSTab tab2(out);
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_) {
535 beta_k = max(beta_PR, ST::zero());
540 if (numConsecutiveIters < 2) {
543 else if (beta_PR < -beta_FR)
545 else if (ST::magnitude(beta_PR) <= beta_FR)
556 TEUCHOS_TEST_FOR_EXCEPT(
true);
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);
645 if (!is_null(g_grad_km1))
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 );
668 if (!is_null(numIters_out)) {
669 *numIters_out = numIters_;
672 if (numIters_ == maxIters_) {
673 *out <<
"\nMax nonlinear CG iterations exceeded!\n";
679 else if(fatalLinesearchFailure) {
692 #endif // OPTIPACK_NONLINEAR_CG_DEF_HPP bool get_and_conv_tests() const
const std::string and_conv_tests_name
const std::string alpha_init_name
ScalarMag get_alpha_init() const
const bool and_conv_tests_default
const double g_mag_default
const std::string solverType_name
const int maxIters_default
bool get_alpha_reinit() const
const double g_grad_tol_default
const std::string g_reduct_tol_name
const double g_reduct_tol_default
static RCP< Teuchos::ParameterEntryValidator > solverType_validator_
void setParameterList(RCP< ParameterList > const ¶mList)
const double alpha_init_default
RCP< const ParameterList > getValidParameters() const
NonlinearCG()
Construct with default parameters.
const std::string g_mag_name
const ESolverTypes solverType_default_integral_val
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.
const int minIters_default
NonlinearCGUtils::ESolverTypes get_solverType() const
ScalarMag get_g_grad_tol() const
Fletcher-Reeves Polak-Ribiere Hybrid Method.
ScalarMag get_g_mag() const
const std::string alpha_reinit_name
const std::string solverType_default
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.
const bool alpha_reinit_default
const std::string g_grad_tol_name
ScalarTraits< Scalar >::magnitudeType ScalarMag