Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DampenedNewtonNonlinearSolver.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
43#define THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
44
45#include "Thyra_NonlinearSolverBase.hpp"
46#include "Thyra_ModelEvaluatorHelpers.hpp"
47#include "Thyra_TestingTools.hpp"
50#include "Teuchos_VerboseObject.hpp"
51#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
52#include "Teuchos_StandardParameterEntryValidators.hpp"
53#include "Teuchos_as.hpp"
54
55
56namespace Thyra {
57
58
79template <class Scalar>
81public:
82
86 typedef typename ST::magnitudeType ScalarMag;
89
92
94 STANDARD_MEMBER_COMPOSITION_MEMBERS( int, defaultMaxNewtonIterations );
95
97 STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, useDampenedLineSearch );
98
100 STANDARD_MEMBER_COMPOSITION_MEMBERS( Scalar, armijoConstant );
101
103 STANDARD_MEMBER_COMPOSITION_MEMBERS( int, maxLineSearchIterations );
104
107 const ScalarMag defaultTol = 1e-2
108 ,const int defaultMaxNewtonIterations = 1000
109 ,const bool useDampenedLineSearch = true
110 ,const Scalar armijoConstant = 1e-4
111 ,const int maxLineSearchIterations = 20
112 );
113
117
120
122 void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
131
133
136
138 void setModel(
139 const RCP<const ModelEvaluator<Scalar> > &model
140 );
146 const SolveCriteria<Scalar> *solveCriteria,
147 VectorBase<Scalar> *delta
148 );
152 bool is_W_current() const;
154 RCP<LinearOpWithSolveBase<Scalar> > get_nonconst_W(const bool forceUpToDate);
158 void set_W_is_current(bool W_is_current);
159
161
162private:
163
167 RCP<VectorBase<Scalar> > current_x_;
168 bool J_is_current_;
169
170};
171
172// ////////////////////////
173// Defintions
174
175template <class Scalar>
177 const ScalarMag my_defaultTol
178 ,const int my_defaultMaxNewtonIterations
179 ,const bool my_useDampenedLineSearch
180 ,const Scalar my_armijoConstant
181 ,const int my_maxLineSearchIterations
182 )
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)
189{}
190
191template <class Scalar>
194{
195 static RCP<const Teuchos::ParameterList> validSolveCriteriaExtraParameters;
196 if(!validSolveCriteriaExtraParameters.get()) {
198 paramList = Teuchos::rcp(new Teuchos::ParameterList);
199 paramList->set("Max Iters",int(1000));
200 validSolveCriteriaExtraParameters = paramList;
201 }
202 return validSolveCriteriaExtraParameters;
203}
204
205// Overridden from Teuchos::ParameterListAcceptor
206
207template<class Scalar>
209 RCP<Teuchos::ParameterList> const& paramList
210 )
211{
212 using Teuchos::get;
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);
218#ifdef TEUCHOS_DEBUG
219 paramList_->validateParameters(*getValidParameters(),0);
220#endif // TEUCHOS_DEBUG
221}
222
223template<class Scalar>
229
230template<class Scalar>
233{
234 RCP<Teuchos::ParameterList> _paramList = paramList_;
235 paramList_ = Teuchos::null;
236 return _paramList;
237}
238
239template<class Scalar>
242{
243 return paramList_;
244}
245
246template<class Scalar>
249{
250 using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
252 if (is_null(validPL)) {
254 pl = Teuchos::parameterList();
255 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement!");
256 Teuchos::setupVerboseObjectSublist(&*pl);
257 validPL = pl;
258 }
259 return validPL;
260}
261
262// Overridden from NonlinearSolverBase
263
264template <class Scalar>
266 const RCP<const ModelEvaluator<Scalar> > &model
267 )
268{
269 TEUCHOS_TEST_FOR_EXCEPT(model.get()==NULL);
270 model_ = model;
271 J_ = Teuchos::null;
272 current_x_ = Teuchos::null;
273 J_is_current_ = false;
274}
275
276template <class Scalar>
279{
280 return model_;
281}
282
283template <class Scalar>
286 VectorBase<Scalar> *x_inout
287 ,const SolveCriteria<Scalar> *solveCriteria
288 ,VectorBase<Scalar> *delta
289 )
290{
291
292 using std::endl;
293 using Teuchos::as;
294
295 // Validate input
296#ifdef TEUCHOS_DEBUG
297 TEUCHOS_TEST_FOR_EXCEPT(0==x_inout);
299 "DampenedNewtonNonlinearSolver<Scalar>::solve(...)",
300 *x_inout->space(), *model_->get_x_space() );
301#endif
302
303 // Get the output stream and verbosity level
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";
315 }
316
317 // Initialize storage for algorithm
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());
325
326 // Get convergence criteria
327 ScalarMag tol = this->defaultTol();
328 int maxIters = this->defaultMaxNewtonIterations();
329 if(solveCriteria && !solveCriteria->solveMeasureType.useDefault()) {
332 ,"DampenedNewtonNonlinearSolver<Scalar>::solve(...): Error, can only support resudual-based"
333 " convergence criteria!");
334 tol = solveCriteria->requestedTol;
335 if(solveCriteria->extraParameters.get()) {
336 solveCriteria->extraParameters->validateParameters(*getValidSolveCriteriaExtraParameters());
337 maxIters = solveCriteria->extraParameters->get("Max Iters",int(maxIters));
338 }
339 }
340
341 if(out.get() && showNewtonDetails)
342 *out << "\nCompute the initial starting point ...\n";
343
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_;
350 }
351
352 // Peform the Newton iterations
353 int newtonIter, num_residual_evals = 1;
354 SolveStatus<Scalar> solveStatus;
356
357 for( newtonIter = 1; newtonIter <= maxIters; ++newtonIter ) {
358
359 if(out.get() && showNewtonDetails) *out << "\n*** newtonIter = " << newtonIter << endl;
360
361 // Check convergence
362 if(out.get() && showNewtonDetails) *out << "\nChecking for convergence ... : ";
363 const Scalar phi = scalarProd(*f,*f), sqrt_phi = ST::squareroot(phi); // merit function: phi(f) = <f,f>
364 solveStatus.achievedTol = sqrt_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;
371 if(isConverged) {
372 if(x_inout != x.get()) assign( ptr(x_inout), *x ); // Assign the solution if we have to
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";
378 }
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();
383 break;
384 }
385 if(out.get() && showNewtonDetails) *out << "\nWe have to keep going :-(\n";
386
387 // Compute the Jacobian if we have not already
388 if(newtonIter > 1) {
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_;
392 }
393
394 // Compute the newton step: dx = -inv(J)*f
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() ); // Initial guess for the linear solve
398 J_->solve(NOTRANS,*f,dx.ptr()); // Solve: J*dx = f
399 Vt_S( dx.ptr(), Scalar(-ST::one()) ); // dx *= -1.0
400 Vp_V( ee.ptr(), *dx); // ee += dx
401 if(out.get() && showNewtonDetails) *out << "\n||dx||inf = " << norm_inf(*dx) << endl;
402 if(out.get() && dumpAll) *out << "\ndy =\n" << *dx;
403
404 // Perform backtracking armijo line search
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; // D(phi(x+alpha*dx))/D(alpha) at alpha=0.0 => -2.0*<f,c>: where dx = -inv(J)*f
408 Scalar alpha = 1.0; // Try a full step initially since it will eventually be accepted near solution
409 int lineSearchIter;
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;
414 // x_new = x + alpha*dx
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;
418 // Compute the residual at the updated point
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";
425 alpha *= 0.1;
426 continue;
427 }
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";
440 break;
441 }
442 if(out.get() && showNewtonDetails) *out << "\nBacktracking (alpha = 0.5*alpha) ...\n";
443 alpha *= 0.5;
444 }
445
446 // Check for line search failure
447 if( lineSearchIter > maxLineSearchIterations() ) {
448 std::ostringstream oss;
449 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;
454 goto exit;
455 }
456
457 // Take the Newton step
458 std::swap<RCP<VectorBase<Scalar> > >( x_new, x ); // Now x is current point!
459
460 }
461
462exit:
463
464 if(out.get() && showNewtonIters) *out
465 << "\n[Final] newton_iters="<<newtonIter<<", num_residual_evals="<<num_residual_evals<<"\n";
466
467 if(newtonIter > maxIters) {
468 std::ostringstream oss;
469 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;
474 }
475
476 if(x_inout != x.get()) assign( ptr(x_inout), *x ); // Assign the final point
477 if(delta != NULL) assign( ptr(delta), *ee );
478 current_x_ = x_inout->clone_v(); // Remember the final point
479 J_is_current_ = newtonIter==1; // J is only current with x if initial point was converged!
480
481 if(out.get() && showNewtonDetails) *out
482 << "\n*** Ending dampended Newton solve." << endl;
483
484 return solveStatus;
485
486}
487
488template <class Scalar>
491{
492 return current_x_;
493}
494
495template <class Scalar>
497{
498 return J_is_current_;
499}
500
501template <class Scalar>
504{
505 if (forceUpToDate) {
506 TEUCHOS_TEST_FOR_EXCEPT(forceUpToDate);
507 }
508 return J_;
509}
510
511template <class Scalar>
514{
515 return J_;
516}
517
518template <class Scalar>
520{
521 J_is_current_ = W_is_current;
522}
523
524
525} // namespace Thyra
526
527
528#endif // THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
Ptr< T > ptr() const
T * get() const
Exception type thrown on an catastrophic solve failure.
Simple dampended Newton solver using a Armijo line search :-)
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
STANDARD_MEMBER_COMPOSITION_MEMBERS(int, maxLineSearchIterations)
Set the maximum number of backtracking line search iterations to take.
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.
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
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)
T_To & dyn_cast(T_From &from)
#define TEUCHOS_OSTAB_DIFF(DIFF)
#define TEUCHOS_OSTAB
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...