[gnumeric] Solver: return Hessian as a GnmMatrix.

[Morten Welinder <mortenw src gnome org>]

commit 1c54df66562719fb9623234e7302586dd3000b10
Author: Morten Welinder <terra gnome org>
Date:   Sun Oct 2 16:03:31 2016 -0400

    Solver: return Hessian as a GnmMatrix.

 src/tools/gnm-solver.c |   35 ++++++++++++++++++++++-------------
 src/tools/gnm-solver.h |    2 +-
 2 files changed, 23 insertions(+), 14 deletions(-)
---
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index c162b88..8303be4 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -9,6 +9,7 @@
 #include "rangefunc.h"
 #include "ranges.h"
 #include "gutils.h"
+#include "mathfunc.h"
 #include "gnm-solver.h"
 #include "workbook-view.h"
 #include "workbook-control.h"
@@ -2143,29 +2144,37 @@ gnm_solver_has_analytic_hessian (GnmSolver *sol)
  * will be n+(n-1)+...+2+1 elements long containing the triangular
  * Hessian.  Use symmetry to obtain the full Hessian.
  */
-gnm_float *
+GnmMatrix *
 gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs)
 {
-       int i, hlen;
-       gnm_float *H;
+       int i, j, k;
+       GnmMatrix *H;
        GnmEvalPos ep;
+       int const n = sol->input_cells->len;
 
        if (!gnm_solver_has_analytic_hessian (sol))
                return NULL;
 
        gnm_solver_set_vars (sol, xs);
 
-       hlen = sol->hessian->len;
-       H = g_new (gnm_float, hlen);
+       H = gnm_matrix_new (n, n);
        eval_pos_init_cell (&ep, sol->target);
-       for (i = 0; i < hlen; i++) {
-               GnmExprTop const *te = g_ptr_array_index (sol->hessian, i);
-               GnmValue *v = gnm_expr_top_eval
-                       (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
-               H[i] = VALUE_IS_NUMBER (v) ? value_get_as_float (v) : gnm_nan;
-               if (sol->flip_sign)
-                       H[i] = 0 - H[i];
-               value_release (v);
+       for (i = k = 0; i < n; i++) {
+               for (j = i; j < n; j++, k++) {
+                       GnmExprTop const *te =
+                               g_ptr_array_index (sol->hessian, k);
+                       GnmValue *v = gnm_expr_top_eval
+                               (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
+                       gnm_float x = VALUE_IS_NUMBER (v)
+                               ? value_get_as_float (v)
+                               : gnm_nan;
+                       if (sol->flip_sign)
+                               x = 0 - x;
+                       value_release (v);
+
+                       H->data[i][j] = x;
+                       H->data[j][i] = x;
+               }
        }
 
        return H;
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index 63eea58..7f33146 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -309,7 +309,7 @@ gboolean gnm_solver_has_analytic_gradient (GnmSolver *sol);
 gnm_float *gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs);
 
 gboolean gnm_solver_has_analytic_hessian (GnmSolver *sol);
-gnm_float *gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs);
+GnmMatrix *gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs);
 
 gnm_float gnm_solver_line_search (GnmSolver *sol,
                                  gnm_float const *x0, gnm_float const *dir,