GNOME.org
 

[gnumeric] Add EIGEN function to calculate eigenvalues and eigenvectors of real symmetric matrices

commit 8f5f51a31ebe08e220282f459f3c5be294380939
Author: Andreas J. Guelzow <aguelzow pyrshep ca>
Date:   Mon Dec 21 23:07:00 2009 -0700

    Add EIGEN function to calculate eigenvalues and eigenvectors of real symmetric matrices
    
    2009-12-21  Andreas J. Guelzow <aguelzow pyrshep ca>
    
    	* src/mathfunc.h (gnm_matrix_eigen): new
    	* src/mathfunc.c (gnm_matrix_eigen): new
    	(gnm_matrix_eigen_max_index): new
    	(gnm_matrix_eigen_update): new
    
    2009-121  Andreas J. Guelzow <aguelzow pyrshep ca>
    
    	* plugin.xml.in: add EIGEN
    	* functions.c (help_eigen): new
    	(new_matrix): new
    	)compare_doubles): new
    	(gnumeric_eigen): new
    	(math_functions): add EIGEN

 ChangeLog                     |    7 +++
 NEWS                          |    2 +
 plugins/fn-math/ChangeLog     |    9 +++
 plugins/fn-math/functions.c   |  114 ++++++++++++++++++++++++++++++++++++++++
 plugins/fn-math/plugin.xml.in |    1 +
 src/mathfunc.c                |  115 +++++++++++++++++++++++++++++++++++++++++
 src/mathfunc.h                |    2 +
 7 files changed, 250 insertions(+), 0 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 21c894a..bbd0bbc 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2009-12-21  Andreas J. Guelzow <aguelzow pyrshep ca>
+
+	* src/mathfunc.h (gnm_matrix_eigen): new
+	* src/mathfunc.c (gnm_matrix_eigen): new
+	(gnm_matrix_eigen_max_index): new
+	(gnm_matrix_eigen_update): new
+	
 2009-12-20  Andreas J. Guelzow <aguelzow pyrshep ca>
 
 	* src/sheet-style.h (sheet_style_set_list): change arguments
diff --git a/NEWS b/NEWS
index ec8c5df..c9c5ba3 100644
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,8 @@ Gnumeric 1.9.18
 
 Andreas:
 	* Add paste special flip horizontally and vertically [#393367]
+	* Add EIGEN function to calculate eigenvalues and eigenvectors
+	  of real symmetric matrices
 
 Jean
 	* Fix import export of line type in scatter plots. [#605043]
diff --git a/plugins/fn-math/ChangeLog b/plugins/fn-math/ChangeLog
index a052eed..8521e76 100644
--- a/plugins/fn-math/ChangeLog
+++ b/plugins/fn-math/ChangeLog
@@ -1,3 +1,12 @@
+2009-121  Andreas J. Guelzow <aguelzow pyrshep ca>
+
+	* plugin.xml.in: add EIGEN
+	* functions.c (help_eigen): new
+	(new_matrix): new
+	)compare_doubles): new
+	(gnumeric_eigen): new
+	(math_functions): add EIGEN
+
 2009-12-15  Morten Welinder <terra gnome org>
 
 	* Release 1.9.17
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index ba2fcc1..936dbed 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -3020,6 +3020,117 @@ done:
 
 /***************************************************************************/
 
+static GnmFuncHelp const help_eigen[] = {
+        { GNM_FUNC_HELP_NAME, F_("EIGEN:eigenvalues and eigenvectors of the symmetric @{matrix}")},
+        { GNM_FUNC_HELP_ARG, F_("matrix:a symmetric matrix")},
+	{ GNM_FUNC_HELP_NOTE, F_("If @{matrix} is not symmetric, EIGEN returns #NUM!") },
+	{ GNM_FUNC_HELP_NOTE, F_("If @{matrix} does not contain an equal number of columns and rows, EIGEN returns #VALUE!") },
+        { GNM_FUNC_HELP_END}
+};
+
+
+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;
+}
+
+int
+compare_doubles (const void *a, const void *b)
+{
+  const gnm_float *da = (const double *) a;
+  const gnm_float *db = (const double *) b;
+
+  if (*da > *db)
+	  return -1;
+  else if (*da == *db)
+	  return 0;
+  else
+	  return 1;
+}
+
+
+typedef struct {
+	gnm_float val;
+	int index;
+} gnumeric_eigen_ev_t;
+
+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;
+	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);
+
+	/* 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);
+			}
+	
+	eigenvectors = new_matrix (rows, cols);
+	eigenvalues = g_new0 (gnm_float, 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);
+	}
+
+	/* 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;
+	}
+	qsort (ev_sort, cols, sizeof (gnumeric_eigen_ev_t), compare_doubles);
+
+	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->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);
+		}
+	}
+	free_matrix (matrix, cols, rows);
+	free_matrix (eigenvectors, cols, rows);
+	g_free (eigenvalues);
+	g_free (ev_sort);
+
+	return res;
+}
+
+
+/***************************************************************************/
+
 GnmFuncDescriptor const math_functions[] = {
 	{ "abs",     "f",     help_abs,
 	  gnumeric_abs, NULL, NULL, NULL, NULL,
@@ -3283,6 +3394,9 @@ GnmFuncDescriptor const math_functions[] = {
 	  gnumeric_munit, NULL, NULL, NULL, NULL,
 	  GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_COMPLETE,
 	  GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
+	{ "eigen","A",      help_eigen,
+	  gnumeric_eigen, NULL, NULL, NULL, NULL,
+	  GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE},
 #if 0
 	{ "logmdeterm", "A|si",
 	  help_logmdeterm, gnumeric_logmdeterm, NULL, NULL, NULL },
diff --git a/plugins/fn-math/plugin.xml.in b/plugins/fn-math/plugin.xml.in
index a79e8ed..d740914 100644
--- a/plugins/fn-math/plugin.xml.in
+++ b/plugins/fn-math/plugin.xml.in
@@ -38,6 +38,7 @@
 				<function name="csc"/>
 				<function name="csch"/>
 				<function name="degrees"/>
+				<function name="eigen"/>
 				<function name="even"/>
 				<function name="exp"/>
 				<function name="expm1"/>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 42a3608..086961d 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7795,6 +7795,121 @@ mmult (gnm_float *A, gnm_float *B, int cols_a, int rows_a, int cols_b,
 	}
 }
 
+/***************************************************************************/
+
+static int
+gnm_matrix_eigen_max_index (gnm_float *row, int row_n, int size) 
+{
+	int i, res = row_n + 1;
+	gnm_float max = gnm_abs (row[res]);
+	for (i = res + 1; i < size; i++)
+		if (gnm_abs (row[i]) > max) {
+			res = i;
+			max = gnm_abs (row[i]);
+		}
+	return res;
+}
+
+static void
+gnm_matrix_eigen_rotate (gnm_float **matrix, int k, int l, int i, int j, gnm_float c, gnm_float s)
+{
+	gnm_float x = c * matrix[k][l] - s * matrix[i][j];
+	gnm_float y = s * matrix[k][l] + c * matrix[i][j];
+
+	matrix[k][l] = x;
+	matrix[i][j] = y;
+}
+
+static void
+gnm_matrix_eigen_update (int k, gnm_float t, gnm_float *eigenvalues, gboolean *changed, int *state)
+{
+	gnm_float y = eigenvalues[k];
+	eigenvalues[k] += t;
+	if (changed[k] && y == eigenvalues[k]) {
+		changed[k] = FALSE;
+		(*state)--;
+	} else if ((!changed[k]) && (y != eigenvalues[k])) {
+		changed[k] = TRUE;
+		(*state)++;
+	}
+}
+
+/* Calculates the eigenvalues and eigenvectors of a real symmetric matrix.
+ */
+gboolean    
+gnm_matrix_eigen (gnm_float **matrix, gnm_float **eigenvectors, gnm_float *eigenvalues, int size) 
+{
+	int i, state = size, *ind;
+	gboolean *changed;
+	int counter = 0;
+
+	ind = g_new (int, size);
+	changed =  g_new (gboolean, size);
+	
+	for (i = 0; i < size; i++) {
+		int j;
+		for (j = 0; j < size; j++)
+			eigenvectors[j][i] = 0.;
+		eigenvectors[i][i] = 1.;
+		eigenvalues[i] = matrix[i][i];
+		ind[i] = gnm_matrix_eigen_max_index (matrix[i], i, size);
+		changed[i] = TRUE;
+	}
+
+	while (state != 0) {
+		int k, l, m = 0;
+		gnm_float c, s, y, pivot, t;
+		
+		counter++;
+		if (counter > 400000) {
+			g_free (ind);
+			g_free (changed);
+			g_print ("gnm_matrix_eigen exceeded iterations\n"); 
+			return FALSE;
+		}
+		for (k = 1; k < (size-1); k++)
+			if (gnm_abs (matrix[k][ind[k]]) > gnm_abs (matrix[m][ind[m]]))
+				m = k;
+		l = ind[m];
+		pivot = matrix[m][l];
+		/* pivot is (m,l) */
+		
+		y = (eigenvalues[l] - eigenvalues[m])/2;
+		t = gnm_abs (y) + gnm_sqrt (pivot*pivot+y*y);
+		s = gnm_sqrt (pivot*pivot+t*t);
+		c = t/s;
+		s = pivot/s;
+		t = pivot * pivot /t;
+		if (y < 0) {
+			s = -s;
+			t = -t;
+		}
+		matrix[m][l] = 0.;
+		gnm_matrix_eigen_update (m, -t, eigenvalues, changed, &state);
+		gnm_matrix_eigen_update (l, t, eigenvalues, changed, &state);
+		for (i = 0; i < m; i++)
+			gnm_matrix_eigen_rotate (matrix, i, m, i, l, c, s);
+		for (i = m + 1; i < l; i++)
+			gnm_matrix_eigen_rotate (matrix, m, i, i, l, c, s);
+		for (i = l + 1; i < size; i++)
+			gnm_matrix_eigen_rotate (matrix, m, i, l, i, c, s);
+		for (i = 0; i < size; i++) {
+			gnm_float x = c * eigenvectors[i][m] - s * eigenvectors[i][l];
+			gnm_float y = s * eigenvectors[i][m] + c * eigenvectors[i][l];
+
+			eigenvectors[i][m] = x;
+			eigenvectors[i][l] = y;
+		}
+		ind[m] = gnm_matrix_eigen_max_index (matrix[m], m, size);
+		ind[l] = gnm_matrix_eigen_max_index (matrix[l], l, size);
+	}
+
+	g_free (ind);
+	g_free (changed);
+
+	return TRUE;
+}
+
 /* ------------------------------------------------------------------------- */
 
 #ifdef NEED_FAKE_ERFGNUM
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 38796be..3d335f2 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -184,6 +184,8 @@ gnm_float discpfuncinverter (gnm_float p, const gnm_float shape[],
 void    mmult (gnm_float *A, gnm_float *B, int cols_a, int rows_a, int cols_b,
 	       gnm_float *product);
 
+gboolean gnm_matrix_eigen (gnm_float **matrix, gnm_float **eigenvectors,
+			   gnm_float *eigenvalues, int size);
 /* ------------------------------------------------------------------------- */
 
 gnm_float combin (gnm_float n, gnm_float k);
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]