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;