27 const Problem& P,
const Preconditioner&
T,
28 const Field& u_old,
Float Tu_old, Field& delta_u, Real slope, Real norm_delta_u_max,
31 const Float alpha = 1e-4;
32 const Float eps_mach = std::numeric_limits<Float>::epsilon();
33 Float norm_delta_u = P.space_norm(delta_u);
34 if (norm_delta_u > norm_delta_u_max) {
35 Float c = norm_delta_u_max/norm_delta_u;
36 if (p_derr) *p_derr <<
"# damped-Newton/backtrack: warning: delta_u bounded by factor " << c << std::endl << std::flush;
40 Float norm_u = P.space_norm(
u);
41 if (norm_u < norm_delta_u) norm_u = norm_delta_u;
42 Float lambda_min = eps_mach*norm_u/norm_delta_u;
47 lambda = std::max (Real(0.0), std::min (Real(1.0),
lambda));
48 Float lambda_prev = 0;
50 Float Tu_prev_old = 0;
51 for (
size_t k = 0;
true; k++) {
59 }
else if (Tu <= Tu_old + alpha*
lambda*slope) {
63 lambda_next = - 0.5*slope/(Tu - Tu_old - slope);
67 Float z_prev = Tu_prev - Tu_prev_old - lambda_prev*slope;
72 lambda_next = - slope/(2*b);
74 Float Delta = sqr(b) - 3*a*slope;
76 if (p_derr) *p_derr <<
"# damped-Newton/backtrack: warning: machine precision reached" << std::endl << std::flush;
79 lambda_next = (-b + sqrt(Delta))/(3*a);
81 lambda_next = std::min (
lambda/2, lambda_next);
int newton_backtrack(const Problem &P, const Preconditioner &T, const Field &u_old, Float Tu_old, Field &delta_u, Real slope, Real norm_delta_u_max, Field &u, Field &Fu, Real &Tu, Real &lambda, odiststream *p_derr=0)