42#ifndef THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
43#define THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
45#include "Thyra_NonlinearSolverBase.hpp"
46#include "Thyra_ModelEvaluatorHelpers.hpp"
47#include "Thyra_TestingTools.hpp"
48#include "Teuchos_StandardMemberCompositionMacros.hpp"
49#include "Teuchos_StandardCompositionMacros.hpp"
50#include "Teuchos_VerboseObject.hpp"
51#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
52#include "Teuchos_StandardParameterEntryValidators.hpp"
53#include "Teuchos_as.hpp"
79template <
class Scalar>
108 ,
const int defaultMaxNewtonIterations = 1000
109 ,
const bool useDampenedLineSearch =
true
110 ,
const Scalar armijoConstant = 1e-4
111 ,
const int maxLineSearchIterations = 20
175template <
class Scalar>
178 ,
const int my_defaultMaxNewtonIterations
179 ,
const bool my_useDampenedLineSearch
180 ,
const Scalar my_armijoConstant
181 ,
const int my_maxLineSearchIterations
183 :defaultTol_(my_defaultTol)
184 ,defaultMaxNewtonIterations_(my_defaultMaxNewtonIterations)
185 ,useDampenedLineSearch_(my_useDampenedLineSearch)
186 ,armijoConstant_(my_armijoConstant)
187 ,maxLineSearchIterations_(my_maxLineSearchIterations)
188 ,J_is_current_(false)
191template <
class Scalar>
196 if(!validSolveCriteriaExtraParameters.
get()) {
199 paramList->set(
"Max Iters",
int(1000));
200 validSolveCriteriaExtraParameters = paramList;
202 return validSolveCriteriaExtraParameters;
207template<
class Scalar>
214 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
215 paramList_ = paramList;
217 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
219 paramList_->validateParameters(*getValidParameters(),0);
223template<
class Scalar>
230template<
class Scalar>
235 paramList_ = Teuchos::null;
239template<
class Scalar>
246template<
class Scalar>
250 using Teuchos::setDoubleParameter;
using Teuchos::setIntParameter;
252 if (is_null(validPL)) {
254 pl = Teuchos::parameterList();
256 Teuchos::setupVerboseObjectSublist(&*pl);
264template <
class Scalar>
272 current_x_ = Teuchos::null;
273 J_is_current_ =
false;
276template <
class Scalar>
283template <
class Scalar>
299 "DampenedNewtonNonlinearSolver<Scalar>::solve(...)",
300 *x_inout->
space(), *model_->get_x_space() );
311 if(out.
get() && showNewtonIters) {
312 *out <<
"\nBeginning dampended Newton solve of model = " << model_->description() <<
"\n\n";
313 if (!useDampenedLineSearch())
314 *out <<
"\nDoing undampened newton ...\n\n";
318 if(!J_.get()) J_ = model_->create_W();
324 V_S(ee.
ptr(),ST::zero());
328 int maxIters = this->defaultMaxNewtonIterations();
332 ,
"DampenedNewtonNonlinearSolver<Scalar>::solve(...): Error, can only support resudual-based"
333 " convergence criteria!");
336 solveCriteria->
extraParameters->validateParameters(*getValidSolveCriteriaExtraParameters());
337 maxIters = solveCriteria->
extraParameters->get(
"Max Iters",
int(maxIters));
341 if(out.
get() && showNewtonDetails)
342 *out <<
"\nCompute the initial starting point ...\n";
344 eval_f_W( *model_, *x, &*f, &*J_ );
345 if(out.
get() && dumpAll) {
346 *out <<
"\nInitial starting point:\n";
347 *out <<
"\nx =\n" << *x;
348 *out <<
"\nf =\n" << *f;
349 *out <<
"\nJ =\n" << *J_;
353 int newtonIter, num_residual_evals = 1;
357 for( newtonIter = 1; newtonIter <= maxIters; ++newtonIter ) {
359 if(out.
get() && showNewtonDetails) *out <<
"\n*** newtonIter = " << newtonIter << endl;
362 if(out.
get() && showNewtonDetails) *out <<
"\nChecking for convergence ... : ";
363 const Scalar phi = scalarProd(*f,*f), sqrt_phi = ST::squareroot(phi);
365 const bool isConverged = sqrt_phi <= tol;
366 if(out.
get() && showNewtonDetails) *out
367 <<
"sqrt(phi) = sqrt(<f,f>) = ||f|| = " << sqrt_phi << ( isConverged ?
" <= " :
" > " ) <<
"tol = " << tol << endl;
368 if(out.
get() && showNewtonIters) *out
369 <<
"newton_iter="<<newtonIter<<
": Check convergence: ||f|| = "
370 << sqrt_phi << ( isConverged ?
" <= " :
" > " ) <<
"tol = " << tol << ( isConverged ?
", Converged!!!" :
"" ) << endl;
372 if(x_inout != x.
get()) assign( ptr(x_inout), *x );
373 if(out.
get() && showNewtonDetails) {
374 *out <<
"\nWe have converged :-)\n"
375 <<
"\n||x||inf = " << norm_inf(*x) << endl;
376 if(dumpAll) *out <<
"\nx =\n" << *x;
377 *out <<
"\nExiting SimpleNewtonSolver::solve(...)\n";
379 std::ostringstream oss;
380 oss <<
"Converged! ||f|| = " << sqrt_phi <<
", num_newton_iters="<<newtonIter<<
", num_residual_evals="<<num_residual_evals<<
".";
382 solveStatus.
message = oss.str();
385 if(out.
get() && showNewtonDetails) *out <<
"\nWe have to keep going :-(\n";
389 if(out.
get() && showNewtonDetails) *out <<
"\nComputing the Jacobian J_ at current point ...\n";
390 eval_f_W<Scalar>( *model_, *x, NULL, &*J_ );
391 if(out.
get() && dumpAll) *out <<
"\nJ =\n" << *J_;
395 if(out.
get() && showNewtonDetails) *out <<
"\nComputing the Newton step: dx = - inv(J)*f ...\n";
396 if(out.
get() && showNewtonIters) *out <<
"newton_iter="<<newtonIter<<
": Computing Newton step ...\n";
397 assign( dx.
ptr(), ST::zero() );
399 Vt_S( dx.
ptr(), Scalar(-ST::one()) );
400 Vp_V( ee.
ptr(), *dx);
401 if(out.
get() && showNewtonDetails) *out <<
"\n||dx||inf = " << norm_inf(*dx) << endl;
402 if(out.
get() && dumpAll) *out <<
"\ndy =\n" << *dx;
405 if(out.
get() && showNewtonDetails) *out <<
"\nStarting backtracking line search iterations ...\n";
406 if(out.
get() && showNewtonIters) *out <<
"newton_iter="<<newtonIter<<
": Starting backtracking line search ...\n";
407 const Scalar Dphi = -2.0*phi;
410 ++num_residual_evals;
411 for( lineSearchIter = 1; lineSearchIter <= maxLineSearchIterations(); ++lineSearchIter, ++num_residual_evals ) {
413 if(out.
get() && showNewtonDetails) *out <<
"\n*** lineSearchIter = " << lineSearchIter << endl;
415 assign( x_new.
ptr(), *x ); Vp_StV( x_new.
ptr(), alpha, *dx );
416 if(out.
get() && showNewtonDetails) *out <<
"\n||x_new||inf = " << norm_inf(*x_new) << endl;
417 if(out.
get() && dumpAll) *out <<
"\nx_new =\n" << *x_new;
419 eval_f(*model_,*x_new,&*f);
420 if(out.
get() && dumpAll) *out <<
"\nf_new =\n" << *f;
421 const Scalar phi_new = scalarProd(*f,*f), phi_frac = phi + alpha * armijoConstant() * Dphi;
422 if(out.
get() && showNewtonDetails) *out <<
"\nphi_new = <f_new,f_new> = " << phi_new << endl;
424 if(out.
get() && showNewtonDetails) *out <<
"\nphi_new is not a valid number, backtracking (alpha = 0.1*alpha) ...\n";
428 const bool acceptPoint = (phi_new <= phi_frac);
429 if(out.
get() && showNewtonDetails) *out
430 <<
"\nphi_new = " << phi_new << ( acceptPoint ?
" <= " :
" > " )
431 <<
"phi + alpha * eta * Dphi = " << phi <<
" + " << alpha <<
" * " << armijoConstant() <<
" * " << Dphi
432 <<
" = " << phi_frac << endl;
433 if(out.
get() && (showLineSearchIters || (showNewtonIters && acceptPoint))) *out
434 <<
"newton_iter="<<newtonIter<<
", ls_iter="<<lineSearchIter<<
" : "
435 <<
"phi(alpha="<<alpha<<
") = "<<phi_new<<(acceptPoint ?
" <=" :
" >")<<
" armijo_cord = " << phi_frac << endl;
436 if (out.
get() && showNewtonDetails && !useDampenedLineSearch())
437 *out <<
"\nUndamped newton, always accpeting the point!\n";
438 if( acceptPoint || !useDampenedLineSearch() ) {
439 if(out.
get() && showNewtonDetails) *out <<
"\nAccepting the current step with step length alpha = " << alpha <<
"!\n";
442 if(out.
get() && showNewtonDetails) *out <<
"\nBacktracking (alpha = 0.5*alpha) ...\n";
447 if( lineSearchIter > maxLineSearchIterations() ) {
448 std::ostringstream oss;
450 <<
"lineSearchIter = " << lineSearchIter <<
" > maxLineSearchIterations = " << maxLineSearchIterations()
451 <<
": Linear search failure! Algorithm terminated!";
452 solveStatus.
message = oss.str();
453 if(out.
get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
458 std::swap<RCP<VectorBase<Scalar> > >( x_new, x );
464 if(out.
get() && showNewtonIters) *out
465 <<
"\n[Final] newton_iters="<<newtonIter<<
", num_residual_evals="<<num_residual_evals<<
"\n";
467 if(newtonIter > maxIters) {
468 std::ostringstream oss;
470 <<
"newton_iter = " << newtonIter <<
" > maxIters = " << maxIters
471 <<
": Newton algorithm terminated!";
472 solveStatus.
message = oss.str();
473 if( out.
get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
476 if(x_inout != x.
get()) assign( ptr(x_inout), *x );
477 if(delta != NULL) assign( ptr(delta), *ee );
478 current_x_ = x_inout->
clone_v();
479 J_is_current_ = newtonIter==1;
481 if(out.
get() && showNewtonDetails) *out
482 <<
"\n*** Ending dampended Newton solve." << endl;
488template <
class Scalar>
495template <
class Scalar>
498 return J_is_current_;
501template <
class Scalar>
511template <
class Scalar>
518template <
class Scalar>
521 J_is_current_ = W_is_current;
Exception type thrown on an catastrophic solve failure.
Simple dampended Newton solver using a Armijo line search :-)
void set_W_is_current(bool W_is_current)
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
RCP< Teuchos::ParameterList > unsetParameterList()
STANDARD_MEMBER_COMPOSITION_MEMBERS(int, maxLineSearchIterations)
Set the maximum number of backtracking line search iterations to take.
Teuchos::ScalarTraits< ScalarMag > SMT
DampenedNewtonNonlinearSolver(const ScalarMag defaultTol=1e-2, const int defaultMaxNewtonIterations=1000, const bool useDampenedLineSearch=true, const Scalar armijoConstant=1e-4, const int maxLineSearchIterations=20)
STANDARD_MEMBER_COMPOSITION_MEMBERS(bool, useDampenedLineSearch)
The default maximum number of iterations.
Teuchos::ScalarTraits< Scalar > ST
RCP< const VectorBase< Scalar > > get_current_x() const
RCP< const Teuchos::ParameterList > getValidParameters() const
void setModel(const RCP< const ModelEvaluator< Scalar > > &model)
SolveStatus< Scalar > solve(VectorBase< Scalar > *x, const SolveCriteria< Scalar > *solveCriteria, VectorBase< Scalar > *delta)
RCP< LinearOpWithSolveBase< Scalar > > get_nonconst_W(const bool forceUpToDate)
STANDARD_MEMBER_COMPOSITION_MEMBERS(Scalar, armijoConstant)
Set the armijo constant for the line search.
RCP< const ModelEvaluator< Scalar > > getModel() const
static RCP< const Teuchos::ParameterList > getValidSolveCriteriaExtraParameters()
STANDARD_MEMBER_COMPOSITION_MEMBERS(int, defaultMaxNewtonIterations)
The default maximum number of iterations.
STANDARD_MEMBER_COMPOSITION_MEMBERS(ScalarMag, defaultTol)
The default solution tolerance.
RCP< const LinearOpWithSolveBase< Scalar > > get_W() const
RCP< const Teuchos::ParameterList > getParameterList() const
ST::magnitudeType ScalarMag
RCP< Teuchos::ParameterList > getNonconstParameterList()
bool is_W_current() const
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Base class for all nonlinear equation solvers.
Abstract interface for finite-dimensional dense vectors.
virtual RCP< VectorBase< Scalar > > clone_v() const =0
Returns a seprate cloned copy of *this vector with the same values but different storage.
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
Return a smart pointer to the vector space that this vector belongs to.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
@ SOLVE_STATUS_UNCONVERGED
The requested solution criteria has likely not been achieved.
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
@ SOLVE_MEASURE_NORM_RESIDUAL
Norm of the current residual vector (i.e. ||A*x-b||)
@ SOLVE_MEASURE_NORM_RHS
Norm of the right-hand side (i.e. ||b||)
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible.
@ NOTRANS
Use the non-transposed operator.
TypeTo as(const TypeFrom &t)
#define TEUCHOS_OSTAB_DIFF(DIFF)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Simple struct that defines the requested solution criteria for a solve.
ScalarMag requestedTol
The requested solve tolerance (what the client would like to see). Only significant if !...
RCP< ParameterList > extraParameters
Any extra control parameters (e.g. max iterations).
SolveMeasureType solveMeasureType
The type of solve tolerance requested as given in this->requestedTol.
bool useDefault() const
Return if this is a default solve measure (default constructed).
Simple struct for the return status from a solve.
std::string message
A simple one-line message (i.e. no newlines) returned from the solver.
ESolveStatus solveStatus
The return status of the solve.
ScalarMag achievedTol
The maximum final tolerance actually achieved by the (block) linear solve. A value of unknownToleranc...