29 Float x0 = fabs(x[0]);
31 Float wmin = 0, wmax = 1, w0 = 0.5;
33 Float x_star = (x0 < 0.5) ? x0+0.5 : x0-0.5;
40 if (g_star < w_star) {
41 w0 = (w_star + wmin)/2;
44 w0 = (w_star + wmax)/2;
46 w0 = (w_star + w_star2)/2;
52 Float wini_list[] = {w0, wmax, wmin};
57 const Float tol = 1e2*numeric_limits<Float>::epsilon();
59 for (
Float wini : wini_list) {
60 if (
verbose) derr <<
"[harten] wmin="<<wmin<<
", wmax="<<wmax<<
", wini="<<wini<<endl;
61 if (
verbose) derr <<
"[harten] i r dw w, for t0=" <<
t0 <<
", x0="<<x0<< endl;
68 if (
verbose) derr <<
"[harten] " << i <<
" " << fabs(r)
69 <<
" " << dw <<
" " << w << endl;
70 if (fabs(r) <= tol && fabs(dw) <= tol) {
break; }
71 if (++i >=
max_iter) { status = 1;
break; }
72 if (w+dw > 0 && w+dw < wmax) {
74 }
else if (w+dw <= 0) {
76 check_macro (1+w != w,
"Newton iteration: machine precision problem.");
81 if (status == 0)
break;
82 if (
verbose) derr <<
"[harten] failed: restart with another starting point..." << endl;
84 check_macro (status == 0,
"t = " <<
t0 <<
", x = " << x0 <<
" : Newton iteration " << i
85 <<
": precision " << tol <<
" not reached: " << fabs(r));
86 return (x[0] >= 0) ? w : -w;
see the Float page for the full documentation
see the point page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
static const size_t max_iter
Float operator()(const point &x) const
Float f(Float w, Float t, Float x) const
Float df_dw(Float w, Float t, Float x) const