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" 79 template <
class Scalar>
84 typedef Teuchos::ScalarTraits<Scalar>
ST;
88 typedef Teuchos::ScalarTraits<ScalarMag>
SMT;
108 ,
const int defaultMaxNewtonIterations = 1000
109 ,
const bool useDampenedLineSearch =
true 110 ,
const Scalar armijoConstant = 1e-4
111 ,
const int maxLineSearchIterations = 20
115 static RCP<const Teuchos::ParameterList>
142 RCP<const ModelEvaluator<Scalar> >
getModel()
const;
154 RCP<LinearOpWithSolveBase<Scalar> >
get_nonconst_W(
const bool forceUpToDate);
156 RCP<const LinearOpWithSolveBase<Scalar> >
get_W()
const;
164 RCP<Teuchos::ParameterList> paramList_;
165 RCP<const ModelEvaluator<Scalar> > model_;
166 RCP<LinearOpWithSolveBase<Scalar> > J_;
167 RCP<VectorBase<Scalar> > current_x_;
175 template <
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)
191 template <
class Scalar>
192 RCP<const Teuchos::ParameterList>
195 static RCP<const Teuchos::ParameterList> validSolveCriteriaExtraParameters;
196 if(!validSolveCriteriaExtraParameters.get()) {
197 RCP<Teuchos::ParameterList>
198 paramList = Teuchos::rcp(
new Teuchos::ParameterList);
199 paramList->set(
"Max Iters",
int(1000));
200 validSolveCriteriaExtraParameters = paramList;
202 return validSolveCriteriaExtraParameters;
207 template<
class Scalar>
209 RCP<Teuchos::ParameterList>
const& paramList
213 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
214 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
215 paramList_ = paramList;
216 TEUCHOS_TEST_FOR_EXCEPT(
"ToDo: Implement!");
217 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
219 paramList_->validateParameters(*getValidParameters(),0);
220 #endif // TEUCHOS_DEBUG 223 template<
class Scalar>
224 RCP<Teuchos::ParameterList>
230 template<
class Scalar>
231 RCP<Teuchos::ParameterList>
234 RCP<Teuchos::ParameterList> _paramList = paramList_;
235 paramList_ = Teuchos::null;
239 template<
class Scalar>
240 RCP<const Teuchos::ParameterList>
246 template<
class Scalar>
247 RCP<const Teuchos::ParameterList>
250 using Teuchos::setDoubleParameter;
using Teuchos::setIntParameter;
251 static RCP<const Teuchos::ParameterList> validPL;
252 if (is_null(validPL)) {
253 RCP<Teuchos::ParameterList>
254 pl = Teuchos::parameterList();
255 TEUCHOS_TEST_FOR_EXCEPT(
"ToDo: Implement!");
256 Teuchos::setupVerboseObjectSublist(&*pl);
264 template <
class Scalar>
269 TEUCHOS_TEST_FOR_EXCEPT(model.get()==NULL);
272 current_x_ = Teuchos::null;
273 J_is_current_ =
false;
276 template <
class Scalar>
277 RCP<const ModelEvaluator<Scalar> >
283 template <
class Scalar>
297 TEUCHOS_TEST_FOR_EXCEPT(0==x_inout);
299 "DampenedNewtonNonlinearSolver<Scalar>::solve(...)",
300 *x_inout->
space(), *model_->get_x_space() );
304 const RCP<Teuchos::FancyOStream> out = this->getOStream();
305 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
306 const bool showNewtonIters = (verbLevel==Teuchos::VERB_LOW);
307 const bool showLineSearchIters = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM));
308 const bool showNewtonDetails = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
309 const bool dumpAll = (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME));
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();
319 RCP<VectorBase<Scalar> > f = createMember(model_->get_f_space());
320 RCP<VectorBase<Scalar> > x = Teuchos::rcp(x_inout,
false);
321 RCP<VectorBase<Scalar> > dx = createMember(model_->get_x_space());
322 RCP<VectorBase<Scalar> > x_new = createMember(model_->get_x_space());
323 RCP<VectorBase<Scalar> > ee = createMember(model_->get_x_space());
324 V_S(ee.ptr(),ST::zero());
328 int maxIters = this->defaultMaxNewtonIterations();
330 TEUCHOS_TEST_FOR_EXCEPTION(
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() );
398 J_->solve(
NOTRANS,*f,dx.ptr());
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 ) {
412 TEUCHOS_OSTAB_DIFF(lineSearchIter);
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;
423 if( Teuchos::ScalarTraits<Scalar>::isnaninf(phi_new) ) {
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;
488 template <
class Scalar>
489 RCP<const VectorBase<Scalar> >
495 template <
class Scalar>
498 return J_is_current_;
501 template <
class Scalar>
502 RCP<LinearOpWithSolveBase<Scalar> >
506 TEUCHOS_TEST_FOR_EXCEPT(forceUpToDate);
511 template <
class Scalar>
512 RCP<const LinearOpWithSolveBase<Scalar> >
518 template <
class Scalar>
521 J_is_current_ = W_is_current;
528 #endif // THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP bool useDefault() const
Return if this is a default solve measure (default constructed).
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
Return a smart pointer to the vector space that this vector belongs to.
RCP< Teuchos::ParameterList > unsetParameterList()
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
SolveStatus< Scalar > solve(VectorBase< Scalar > *x, const SolveCriteria< Scalar > *solveCriteria, VectorBase< Scalar > *delta)
#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...
ScalarMag requestedTol
The requested solve tolerance (what the client would like to see). Only significant if !this->solveMe...
RCP< const Teuchos::ParameterList > getParameterList() const
Simple dampended Newton solver using a Armijo line search :-)
Teuchos::ScalarTraits< Scalar > ST
std::string message
A simple one-line message (i.e. no newlines) returned from the solver.
Use the non-transposed operator.
Base class for all nonlinear equation solvers.
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
DampenedNewtonNonlinearSolver(const ScalarMag defaultTol=1e-2, const int defaultMaxNewtonIterations=1000, const bool useDampenedLineSearch=true, const Scalar armijoConstant=1e-4, const int maxLineSearchIterations=20)
ST::magnitudeType ScalarMag
RCP< const LinearOpWithSolveBase< Scalar > > get_W() const
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< const VectorBase< Scalar > > get_current_x() const
Teuchos::ScalarTraits< ScalarMag > SMT
Abstract interface for finite-dimensional dense vectors.
RCP< Teuchos::ParameterList > getNonconstParameterList()
Simple struct for the return status from a solve.
Norm of the right-hand side (i.e. ||b||)
void set_W_is_current(bool W_is_current)
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...
bool is_W_current() const
RCP< ParameterList > extraParameters
Any extra control parameters (e.g. max iterations).
void setModel(const RCP< const ModelEvaluator< Scalar > > &model)
STANDARD_MEMBER_COMPOSITION_MEMBERS(ScalarMag, defaultTol)
The default solution tolerance.
static RCP< const Teuchos::ParameterList > getValidSolveCriteriaExtraParameters()
SolveMeasureType solveMeasureType
The type of solve tolerance requested as given in this->requestedTol.
The requested solution criteria has likely been achieved.
Exception type thrown on an catastrophic solve failure.
The requested solution criteria has likely not been achieved.
Simple struct that defines the requested solution criteria for a solve.
RCP< LinearOpWithSolveBase< Scalar > > get_nonconst_W(const bool forceUpToDate)
virtual RCP< VectorBase< Scalar > > clone_v() const =0
Returns a cloned copy of *this vector.
RCP< const ModelEvaluator< Scalar > > getModel() const
Norm of the current residual vector (i.e. ||A*x-b||)