[gnumeric] EIGEN: change to use GnmMatrix.

[Morten Welinder <mortenw src gnome org>]

commit 7e60229a8b22a2d511281e8ce28eba256727bff4
Author: Morten Welinder <terra gnome org>
Date:   Fri Jan 18 15:21:19 2013 -0500

    EIGEN: change to use GnmMatrix.

 ChangeLog                   |    1 +
 plugins/fn-math/ChangeLog   |    1 +
 plugins/fn-math/functions.c |  131 +++++++++++++++++++++----------------------
 src/mathfunc.c              |   16 ++++-
 src/mathfunc.h              |    3 +-
 5 files changed, 80 insertions(+), 72 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 71b5d38..b915085 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -4,6 +4,7 @@
 	(gnm_matrix_is_empty, gnm_matrix_from_value, gnm_matrix_to_value):
 	New matrix support.
 	(gnm_matrix_multiply): Renamed from mmult and changed to use above.
+	(gnm_matrix_eigen): Change to take GnmMatrix arguments.
 
 2013-01-15  Morten Welinder  <terra gnome org>
 
diff --git a/plugins/fn-math/ChangeLog b/plugins/fn-math/ChangeLog
index f501fcb..9ecb33e 100644
--- a/plugins/fn-math/ChangeLog
+++ b/plugins/fn-math/ChangeLog
@@ -3,6 +3,7 @@
 	* functions.c (gnumeric_minverse, gnumeric_mmult)
 	(gnumeric_leverage, gnumeric_linsolve, gnumeric_mdeterm): Simplify
 	using new matrix support.
+	(compare_gnumeric_eigen_ev): Sort first by absolute value.
 
 2013-01-17  Morten Welinder  <terra gnome org>
 
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index cd0bf1e..bdb76ad 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -3141,18 +3141,6 @@ static GnmFuncHelp const help_eigen[] = {
 };
 
 
-static gnm_float **
-new_matrix (int cols, int rows)
-{
-	gnm_float **res = g_new (gnm_float *, rows);
-	int r;
-
-	for (r = 0; r < rows; r++)
-		res[r] = g_new0 (gnm_float, cols);
-
-	return res;
-}
-
 typedef struct {
 	gnm_float val;
 	int index;
@@ -3161,81 +3149,92 @@ typedef struct {
 static int
 compare_gnumeric_eigen_ev (const void *a, const void *b)
 {
-  const gnumeric_eigen_ev_t *da = a;
-  const gnumeric_eigen_ev_t *db = b;
+	const gnumeric_eigen_ev_t *da = a;
+	const gnumeric_eigen_ev_t *db = b;
+	gnm_float ea = da->val;
+	gnm_float eb = db->val;
+
+	/* Compare first by magnitude (descending).  */
+	if (gnm_abs (ea) > gnm_abs (eb))
+		return -1;
+	else if (gnm_abs (ea) < gnm_abs (eb))
+		return +1;
+
+	/* Then by value (still descending.  */
+	if (ea > eb)
+		return -1;
+	else if (ea < eb)
+		return +1;
+	else
+		return 0;
+}
 
-  if (da->val > db->val)
-	  return -1;
-  else if (da->val == db->val)
-	  return 0;
-  else
-	  return 1;
+static gboolean
+symmetric (GnmMatrix const *m)
+{
+	int c, r;
+
+	if (m->cols != m->rows)
+		return FALSE;
+
+	for (c = 0; c < m->cols; ++c)
+		for (r = c + 1; r < m->rows; ++r)
+			if (m->data[r][c] != m->data[c][r])
+				return FALSE;
+
+	return TRUE;
 }
 
 static GnmValue *
 gnumeric_eigen (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
-	GnmEvalPos const * const ep = ei->pos;
-
-	int	r, rows;
-	int	c, cols;
-	GnmValue *res;
-        GnmValue const *values = argv[0];
-	gnm_float **matrix;
-	gnm_float **eigenvectors;
-	gnm_float *eigenvalues;
-	GnmStdError err;
+	GnmMatrix *A = NULL;
+	GnmMatrix *EIG = NULL;
+	gnm_float *eigenvalues = NULL;
+	GnmValue *res = NULL;
+	int i, c, r;
 	gnumeric_eigen_ev_t *ev_sort;
 
-	if (validate_range_numeric_matrix (ep, values, &rows, &cols, &err))
-		return value_new_error_std (ei->pos, err);
-
-	/* Guarantee shape and non-zero size */
-	if (cols != rows || !rows || !cols)
-		return value_new_error_VALUE (ei->pos);
-
-	matrix = value_to_matrix (values, cols, rows, ep);
+	A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+	if (!A) goto out;
 
-	/* Check for symmetry */
-	for (c = 0; c < cols; ++c)
-		for (r = c + 1; r < rows; ++r)
-			if (matrix[r][c] !=  matrix[c][r]) {
-				free_matrix (matrix, cols, rows);
-				return value_new_error_NUM (ei->pos);
-			}
+	if (A->cols != A->rows || gnm_matrix_is_empty (A) || !symmetric (A)) {
+		res = value_new_error_VALUE (ei->pos);
+		goto out;
+	}
 
-	eigenvectors = new_matrix (rows, cols);
-	eigenvalues = g_new0 (gnm_float, cols);
+	EIG = gnm_matrix_new (A->rows, A->cols);
+	eigenvalues = g_new0 (gnm_float, A->cols);
 
-	if (!gnm_matrix_eigen (matrix, eigenvectors, eigenvalues, cols)) {
-		free_matrix (matrix, cols, rows);
-		free_matrix (eigenvectors, cols, rows);
-		g_free (eigenvalues);
-		return value_new_error_NUM (ei->pos);
+	if (!gnm_matrix_eigen (A, EIG, eigenvalues)) {
+		res = value_new_error_NUM (ei->pos);
+		goto out;
 	}
 
 	/* Sorting eigenvalues */
-	ev_sort = g_new (gnumeric_eigen_ev_t, cols);
-	for (c = 0; c < cols; ++c) {
-		ev_sort[c].val = eigenvalues[c];
-		ev_sort[c].index = c;
+	ev_sort = g_new (gnumeric_eigen_ev_t, A->cols);
+	for (i = 0; i < A->cols; i++) {
+		ev_sort[i].val = eigenvalues[i];
+		ev_sort[i].index = i;
 	}
-	qsort (ev_sort, cols, sizeof (gnumeric_eigen_ev_t), compare_gnumeric_eigen_ev);
+	qsort (ev_sort, A->cols, sizeof (gnumeric_eigen_ev_t), compare_gnumeric_eigen_ev);
 
-	res = value_new_array_non_init (cols, rows + 1);
-	for (c = 0; c < cols; ++c) {
-		res->v_array.vals[c] = g_new (GnmValue *, rows + 1);
+	res = value_new_array_non_init (A->cols, A->rows + 1);
+	for (c = 0; c < A->cols; ++c) {
+		res->v_array.vals[c] = g_new (GnmValue *, A->rows + 1);
 		res->v_array.vals[c][0] = value_new_float (eigenvalues[ev_sort[c].index]);
-		for (r = 0; r < rows; ++r) {
-			gnm_float tmp = eigenvectors[r][ev_sort[c].index];
-			res->v_array.vals[c][r+1] = value_new_float (tmp);
+		for (r = 0; r < A->rows; ++r) {
+			gnm_float tmp = EIG->data[r][ev_sort[c].index];
+			res->v_array.vals[c][r + 1] = value_new_float (tmp);
 		}
 	}
-	free_matrix (matrix, cols, rows);
-	free_matrix (eigenvectors, cols, rows);
-	g_free (eigenvalues);
+
 	g_free (ev_sort);
 
+out:
+	if (A) gnm_matrix_free (A);
+	if (EIG) gnm_matrix_free (EIG);
+	g_free (eigenvalues);
 	return res;
 }
 
diff --git a/src/mathfunc.c b/src/mathfunc.c
index e2d48f7..2829649 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -6631,16 +6631,24 @@ gnm_matrix_eigen_update (guint k, gnm_float t, gnm_float *eigenvalues, gboolean
  * magnitude of off-diagonal elements while preserving eigenvalues.
  */
 gboolean
-gnm_matrix_eigen (gnm_float **matrix, gnm_float **eigenvectors, gnm_float *eigenvalues, int size)
+gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenvalues)
 {
 	guint i, state, usize, *ind;
 	gboolean *changed;
 	guint counter = 0;
+	gnm_float **matrix;
+	gnm_float **eigenvectors;
 
-	if (size < 1)
-		return FALSE;
+	g_return_val_if_fail (m != NULL, FALSE);
+	g_return_val_if_fail (m->rows == m->cols, FALSE);
+	g_return_val_if_fail (EIG != NULL, FALSE);
+	g_return_val_if_fail (EIG->rows == EIG->cols, FALSE);
+	g_return_val_if_fail (EIG->rows == m->rows, FALSE);
+
+	matrix = m->data;
+	eigenvectors = EIG->data;
+	usize = m->rows;
 
-	usize = (guint) size;
 	state = usize;
 
 	ind = g_new (guint, usize);
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 075d926..d737b3f 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -161,8 +161,7 @@ gboolean gnm_matrix_is_empty (GnmMatrix const *m);
 
 void gnm_matrix_multiply (GnmMatrix *C, const GnmMatrix *A, const GnmMatrix *B);
 
-gboolean gnm_matrix_eigen (gnm_float **matrix, gnm_float **eigenvectors,
-			   gnm_float *eigenvalues, int size);
+gboolean gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenvalues);
 /* ------------------------------------------------------------------------- */
 
 gnm_float combin (gnm_float n, gnm_float k);