[gnumeric] nlsolve: add occasional Newton steps when we have analytic gradient

[Morten Welinder <mortenw src gnome org>]

commit f8f818dc966c3f412f8414506be81ac668ac4182
Author: Morten Welinder <terra gnome org>
Date:   Sun Sep 11 18:43:09 2016 -0400

    nlsolve: add occasional Newton steps when we have analytic gradient

 plugins/nlsolve/gnm-nlsolve.c |   29 ++++++++++++++++++++++++++---
 1 files changed, 26 insertions(+), 3 deletions(-)
---
diff --git a/plugins/nlsolve/gnm-nlsolve.c b/plugins/nlsolve/gnm-nlsolve.c
index cbcbf40..15d2232 100644
--- a/plugins/nlsolve/gnm-nlsolve.c
+++ b/plugins/nlsolve/gnm-nlsolve.c
@@ -244,10 +244,27 @@ newton_improve (GnmNlsolve *nl, gnm_float *xs, gnm_float *y, gnm_float ymax)
        if (ok) {
                int i;
                gnm_float y2, *xs2 = g_new (gnm_float, n);
-               gnm_float f, best_f = -1;
+               gnm_float best_f = 42;
+               // We try these step sizes.  We really should not need
+               // negative, but if H isn't positive definite it might
+               // work.
+               static const gnm_float fs[] = {
+                       1.0, 0.5, 1.0 / 16,
+                       -1.0, -1.0 / 16,
+               };
+               unsigned ui;
+
+               if (nl->debug) {
+                       int i;
+                       for (i = 0; i < n; i++)
+                               print_vector (NULL, H[i], n);
+                       print_vector ("d", d, n);
+                       print_vector ("g", g, n);
+               }
 
                ok = FALSE;
-               for (f = 1; f > 1e-4; f /= 2) {
+               for (ui = 0 ; ui < G_N_ELEMENTS (fs); ui++) {
+                       gnm_float f = fs[ui];
                        int i;
                        for (i = 0; i < n; i++)
                                xs2[i] = xs[i] - f * d[i];
@@ -266,7 +283,7 @@ newton_improve (GnmNlsolve *nl, gnm_float *xs, gnm_float *y, gnm_float ymax)
                        }
                }
 
-               if (best_f > 0) {
+               if (best_f != 42) {
                        for (i = 0; i < n; i++)
                                xs[i] = xs[i] - best_f * d[i];
                        *y = ymax;
@@ -349,6 +366,12 @@ rosenbrock_iter (GnmNlsolve *nl)
                }
        }
 
+       if (isol->iterations % 100 == 0 &&
+           gnm_solver_has_analytic_gradient (sol)) {
+               if (newton_improve (nl, isol->xk, &isol->yk, isol->yk))
+                       return TRUE;
+       }
+
        if (isol->iterations % 20 == 0) {
                for (i = 0; i < n; i++)
                        for (j = 0; j < n; j++)