[gnumeric] Math: clean up most matrix code.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Math: clean up most matrix code.
- Date: Fri, 18 Jan 2013 19:48:06 +0000 (UTC)
commit e67957cd91f3fe6a3a51fddce7ab9903cd214e0a
Author: Morten Welinder <terra gnome org>
Date: Fri Jan 18 14:47:21 2013 -0500
Math: clean up most matrix code.
ChangeLog | 7 ++
NEWS | 1 +
plugins/fn-math/ChangeLog | 6 +
plugins/fn-math/functions.c | 240 ++++++++++++++++---------------------------
src/mathfunc.c | 109 +++++++++++++++++---
src/mathfunc.h | 16 +++-
6 files changed, 213 insertions(+), 166 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index d2331dd..71b5d38 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2013-01-18 Morten Welinder <terra gnome org>
+
+ * src/mathfunc.c (gnm_matrix_new, gnm_matrix_free)
+ (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.
+
2013-01-15 Morten Welinder <terra gnome org>
* src/stf.c (stf_read_workbook_auto_csvtab): Fix crash for text
diff --git a/NEWS b/NEWS
index dfdfdf7..fbdfbc5 100644
--- a/NEWS
+++ b/NEWS
@@ -26,6 +26,7 @@ Morten:
* Restrict size of MUNIT to avoid thrashing the system. [#625544]
* Fix text import crash.
* Add LEVERAGE function for regression tool. [#691913]
+ * Clean up matrix code.
--------------------------------------------------------------------------
Gnumeric 1.12.0
diff --git a/plugins/fn-math/ChangeLog b/plugins/fn-math/ChangeLog
index 5c3df7f..f501fcb 100644
--- a/plugins/fn-math/ChangeLog
+++ b/plugins/fn-math/ChangeLog
@@ -1,3 +1,9 @@
+2013-01-18 Morten Welinder <terra gnome org>
+
+ * functions.c (gnumeric_minverse, gnumeric_mmult)
+ (gnumeric_leverage, gnumeric_linsolve, gnumeric_mdeterm): Simplify
+ using new matrix support.
+
2013-01-17 Morten Welinder <terra gnome org>
* functions.c (gnumeric_leverage): New function.
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index 6f03d82..174f330 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -2664,22 +2664,6 @@ value_to_matrix (GnmValue const *v, int cols, int rows, GnmEvalPos const *ep)
return res;
}
-static gnm_float **
-value_to_tmatrix (GnmValue const *v, int cols, int rows, GnmEvalPos const *ep)
-{
- gnm_float **res = g_new (gnm_float *, cols);
- int r, c;
-
- for (c = 0; c < cols; c++) {
- res[c] = g_new (gnm_float, rows);
- for (r = 0; r < rows; r++)
- res[c][r] =
- value_get_as_float (value_area_get_x_y (v, c, r, ep));
- }
-
- return res;
-}
-
static void
free_matrix (gnm_float **mat, G_GNUC_UNUSED int cols, int rows)
{
@@ -2695,38 +2679,24 @@ free_matrix (gnm_float **mat, G_GNUC_UNUSED int cols, int rows)
static GnmValue *
gnumeric_minverse (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
{
- GnmEvalPos const * const ep = ei->pos;
+ GnmMatrix *A = NULL;
+ GnmValue *res = NULL;
- int r, rows;
- int c, cols;
- GnmValue *res;
- GnmValue const *values = argv[0];
- gnm_float **matrix;
- GnmStdError err;
+ A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+ if (!A) goto out;
- 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);
- if (!gnm_matrix_invert (matrix, rows)) {
- free_matrix (matrix, cols, rows);
- return value_new_error_NUM (ei->pos);
+ if (A->cols != A->rows || gnm_matrix_is_empty (A)) {
+ res = value_new_error_VALUE (ei->pos);
+ goto out;
}
- res = value_new_array_non_init (cols, rows);
- for (c = 0; c < cols; ++c) {
- res->v_array.vals[c] = g_new (GnmValue *, rows);
- for (r = 0; r < rows; ++r) {
- gnm_float tmp = matrix[r][c];
- res->v_array.vals[c][r] = value_new_float (tmp);
- }
- }
- free_matrix (matrix, cols, rows);
+ if (gnm_matrix_invert (A->data, A->rows))
+ res = gnm_matrix_to_value (A);
+ else
+ res = value_new_error_NUM (ei->pos);
+out:
+ if (A) gnm_matrix_free (A);
return res;
}
@@ -2870,55 +2840,30 @@ static GnmFuncHelp const help_mmult[] = {
static GnmValue *
gnumeric_mmult (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
{
- GnmEvalPos const * const ep = ei->pos;
- int r, rows_a, rows_b;
- int c, cols_a, cols_b;
- GnmValue *res;
- GnmValue const *values_a = argv[0];
- GnmValue const *values_b = argv[1];
- gnm_float *A, *B, *product;
- GnmStdError err;
-
- if (validate_range_numeric_matrix (ep, values_a, &rows_a, &cols_a, &err) ||
- validate_range_numeric_matrix (ep, values_b, &rows_b, &cols_b, &err))
- return value_new_error_std (ei->pos, err);
-
- /* Guarantee shape and non-zero size */
- if (cols_a != rows_b || !rows_a || !rows_b || !cols_a || !cols_b)
- return value_new_error_VALUE (ei->pos);
-
- res = value_new_array_non_init (cols_b, rows_a);
+ GnmMatrix *A = NULL;
+ GnmMatrix *B = NULL;
+ GnmMatrix *C = NULL;
+ GnmValue *res = NULL;
- A = g_new (gnm_float, cols_a * rows_a);
- B = g_new (gnm_float, cols_b * rows_b);
- product = g_new (gnm_float, rows_a * cols_b);
+ A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+ if (!A) goto out;
- for (c = 0; c < cols_a; c++)
- for (r = 0; r < rows_a; r++) {
- GnmValue const * a =
- value_area_get_x_y (values_a, c, r, ep);
- A[r + c * rows_a] = value_get_as_float (a);
- }
-
- for (c = 0; c < cols_b; c++)
- for (r = 0; r < rows_b; r++) {
- GnmValue const * b =
- value_area_get_x_y (values_b, c, r, ep);
- B[r + c * rows_b] = value_get_as_float (b);
- }
+ B = gnm_matrix_from_value (argv[1], &res, ei->pos);
+ if (!B) goto out;
- mmult (A, B, cols_a, rows_a, cols_b, product);
-
- for (c = 0; c < cols_b; c++) {
- res->v_array.vals[c] = g_new (GnmValue *, rows_a);
- for (r = 0; r < rows_a; r++)
- res->v_array.vals[c][r] =
- value_new_float (product[r + c * rows_a]);
+ if (A->cols != B->rows || gnm_matrix_is_empty (A) || gnm_matrix_is_empty (B)) {
+ res = value_new_error_VALUE (ei->pos);
+ goto out;
}
- g_free (A);
- g_free (B);
- g_free (product);
+ C = gnm_matrix_new (A->rows, B->cols);
+ gnm_matrix_multiply (C, A, B);
+ res = gnm_matrix_to_value (C);
+
+out:
+ if (A) gnm_matrix_free (A);
+ if (B) gnm_matrix_free (B);
+ if (C) gnm_matrix_free (C);
return res;
}
@@ -2937,36 +2882,33 @@ static GnmFuncHelp const help_leverage[] = {
static GnmValue *
gnumeric_leverage (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
{
- GnmEvalPos const * const ep = ei->pos;
- int rows, cols;
+ GnmMatrix *A = NULL;
+ GnmValue *res = NULL;
GORegressionResult regres;
- gnm_float **A;
- GnmStdError err;
- GnmValue const *mat = argv[0];
gnm_float *x;
- GnmValue *res;
- if (validate_range_numeric_matrix (ep, mat, &rows, &cols, &err))
- return value_new_error_std (ei->pos, err);
+ A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+ if (!A) goto out;
- /* Guarantee shape and non-zero size */
- if (!cols || !rows)
- return value_new_error_VALUE (ei->pos);
+ if (gnm_matrix_is_empty (A)) {
+ res = value_new_error_VALUE (ei->pos);
+ goto out;
+ }
+
+ x = g_new (gnm_float, A->rows);
- A = value_to_matrix (mat, cols, rows, ep);
- x = g_new (gnm_float, rows);
- regres = gnm_linear_regression_leverage (A, x, rows, cols);
- free_matrix (A, cols, rows);
+ regres = gnm_linear_regression_leverage (A->data, x, A->rows, A->cols);
if (regres != GO_REG_ok && regres != GO_REG_near_singular_good) {
- res = value_new_error_VALUE (ei->pos);
+ res = value_new_error_NUM (ei->pos);
} else {
+ int x_rows = A->rows, x_cols = 1;
int c, r;
- res = value_new_array_non_init (1, rows);
- for (c = 0; c < 1; c++) {
- res->v_array.vals[c] = g_new (GnmValue *, rows);
- for (r = 0; r < rows; r++)
+ res = value_new_array_non_init (x_cols, x_rows);
+ for (c = 0; c < x_cols; c++) {
+ res->v_array.vals[c] = g_new (GnmValue *, x_rows);
+ for (r = 0; r < x_rows; r++)
res->v_array.vals[c][r] =
value_new_float (x[r]);
}
@@ -2974,6 +2916,8 @@ gnumeric_leverage (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
g_free (x);
+out:
+ if (A) gnm_matrix_free (A);
return res;
}
@@ -2994,43 +2938,40 @@ static GnmFuncHelp const help_linsolve[] = {
static GnmValue *
gnumeric_linsolve (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
{
- GnmEvalPos const * const ep = ei->pos;
- int mrows, mcols, crows, ccols;
- GORegressionResult regres;
- gnm_float **A;
- gnm_float **b;
- GnmStdError err;
- GnmValue const *mat = argv[0];
- GnmValue const *col = argv[1];
+ GnmMatrix *A = NULL;
+ GnmMatrix *B = NULL;
gnm_float *x;
- GnmValue *res;
+ GnmValue *res = NULL;
+ int i;
+ GORegressionResult regres;
- if (validate_range_numeric_matrix (ep, mat, &mrows, &mcols, &err))
- return value_new_error_std (ei->pos, err);
+ A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+ if (!A) goto out;
- if (validate_range_numeric_matrix (ep, col, &crows, &ccols, &err))
- return value_new_error_std (ei->pos, err);
+ B = gnm_matrix_from_value (argv[1], &res, ei->pos);
+ if (!B) goto out;
- /* Guarantee shape and non-zero size */
- if (mrows != mcols || mrows != crows || ccols != 1 || !mcols)
- return value_new_error_VALUE (ei->pos);
+ if (A->cols != A->rows || gnm_matrix_is_empty (A) ||
+ B->cols != 1 || B->rows != A->rows || gnm_matrix_is_empty (B)) {
+ res = value_new_error_VALUE (ei->pos);
+ goto out;
+ }
+
+ x = g_new (gnm_float, B->rows);
+ for (i = 0; i < B->rows; i++)
+ x[i] = B->data[i][0];
- A = value_to_matrix (mat, mcols, mrows, ep);
- b = value_to_tmatrix (col, ccols, crows, ep);
- x = g_new (gnm_float, crows);
- regres = gnm_linear_solve (A, b[0], crows, x);
- free_matrix (A, mcols, mrows);
- free_matrix (b, crows, ccols); /* tmatrix */
+ regres = gnm_linear_solve (A->data, x, A->rows, x);
if (regres != GO_REG_ok && regres != GO_REG_near_singular_good) {
- res = value_new_error_VALUE (ei->pos);
+ res = value_new_error_NUM (ei->pos);
} else {
int c, r;
- res = value_new_array_non_init (ccols, crows);
- for (c = 0; c < ccols; c++) {
- res->v_array.vals[c] = g_new (GnmValue *, crows);
- for (r = 0; r < crows; r++)
+ res = value_new_array_non_init (B->cols, B->rows);
+ for (c = 0; c < B->cols; c++) {
+ res->v_array.vals[c] = g_new (GnmValue *, B->rows);
+ for (r = 0; r < B->rows; r++)
res->v_array.vals[c][r] =
value_new_float (x[r]);
}
@@ -3038,6 +2979,9 @@ gnumeric_linsolve (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
g_free (x);
+out:
+ if (A) gnm_matrix_free (A);
+ if (B) gnm_matrix_free (B);
return res;
}
@@ -3055,26 +2999,22 @@ static GnmFuncHelp const help_mdeterm[] = {
static GnmValue *
gnumeric_mdeterm (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
{
- GnmEvalPos const * const ep = ei->pos;
-
- int rows, cols;
- gnm_float res;
- gnm_float **matrix;
- GnmStdError err;
- GnmValue const *values = argv[0];
+ GnmMatrix *A = NULL;
+ GnmValue *res = NULL;
- if (validate_range_numeric_matrix (ep, values, &rows, &cols, &err))
- return value_new_error_std (ei->pos, err);
+ A = gnm_matrix_from_value (argv[0], &res, ei->pos);
+ if (!A) goto out;
- /* Guarantee shape and non-zero size */
- if (cols != rows || !rows || !cols)
- return value_new_error_VALUE (ei->pos);
+ if (A->cols != A->rows || gnm_matrix_is_empty (A)) {
+ res = value_new_error_VALUE (ei->pos);
+ goto out;
+ }
- matrix = value_to_matrix (values, cols, rows, ep);
- res = gnm_matrix_determinant (matrix, rows);
- free_matrix (matrix, cols, rows);
+ res = value_new_float (gnm_matrix_determinant (A->data, A->rows));
- return value_new_float (res);
+out:
+ if (A) gnm_matrix_free (A);
+ return res;
}
/***************************************************************************/
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 75ff768..e2d48f7 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -47,6 +47,7 @@
#include <string.h>
#include <goffice/goffice.h>
#include <glib/gstdio.h>
+#include <value.h>
#if defined (HAVE_IEEEFP_H) || defined (HAVE_IEEE754_H)
/* Make sure we have this symbol defined, since the existance of either
@@ -6466,29 +6467,111 @@ pow1pm1 (gnm_float x, gnm_float y)
---------------------------------------------------------------------
*/
-/* Calculates the product of two matrixes.
- */
+/* Note the order: y then x. */
+GnmMatrix *
+gnm_matrix_new (int rows, int cols)
+{
+ GnmMatrix *m = g_new (GnmMatrix, 1);
+ int r;
+
+ m->rows = rows;
+ m->cols = cols;
+ m->data = g_new (gnm_float *, rows);
+ for (r = 0; r < rows; r++)
+ m->data[r] = g_new (gnm_float, cols);
+
+ return m;
+}
+
void
-mmult (gnm_float *A, gnm_float *B, int cols_a, int rows_a, int cols_b,
- gnm_float *product)
+gnm_matrix_free (GnmMatrix *m)
+{
+ int r;
+
+ for (r = 0; r < m->rows; r++)
+ g_free (m->data[r]);
+ g_free (m->data);
+ g_free (m);
+}
+
+gboolean
+gnm_matrix_is_empty (GnmMatrix const *m)
+{
+ return m == NULL || m->rows <= 0 || m->cols <= 0;
+}
+
+GnmMatrix *
+gnm_matrix_from_value (GnmValue const *v, GnmValue **perr, GnmEvalPos const *ep)
{
- int c, r, i;
- void *state = gnm_accumulator_start ();
- GnmAccumulator *acc = gnm_accumulator_new ();
+ int cols, rows;
+ int c, r;
+ GnmMatrix *m = NULL;
+
+ *perr = NULL;
+ cols = value_area_get_width (v, ep);
+ rows = value_area_get_height (v, ep);
+ m = gnm_matrix_new (rows, cols);
+ for (r = 0; r < rows; r++) {
+ for (c = 0; c < cols; c++) {
+ GnmValue const *v1 = value_area_fetch_x_y (v, c, r, ep);
+ if (VALUE_IS_ERROR (v1)) {
+ *perr = value_dup (v1);
+ gnm_matrix_free (m);
+ return NULL;
+ }
+
+ m->data[r][c] = value_get_as_float (v1);
+ }
+ }
+ return m;
+}
+
+GnmValue *
+gnm_matrix_to_value (GnmMatrix const *m)
+{
+ GnmValue *res = value_new_array_non_init (m->cols, m->rows);
+ int c, r;
- for (c = 0; c < cols_b; ++c) {
- for (r = 0; r < rows_a; ++r) {
+ for (c = 0; c < m->cols; c++) {
+ res->v_array.vals[c] = g_new (GnmValue *, m->rows);
+ for (r = 0; r < m->rows; r++)
+ res->v_array.vals[c][r] = value_new_float (m->data[r][c]);
+ }
+ return res;
+}
+
+/* C = A * B */
+void
+gnm_matrix_multiply (GnmMatrix *C, const GnmMatrix *A, const GnmMatrix *B)
+{
+ void *state;
+ GnmAccumulator *acc;
+ int c, r, i;
+
+ g_return_if_fail (C != NULL);
+ g_return_if_fail (A != NULL);
+ g_return_if_fail (B != NULL);
+ g_return_if_fail (C->rows == A->rows);
+ g_return_if_fail (C->cols == B->cols);
+ g_return_if_fail (A->cols == B->rows);
+
+ state = gnm_accumulator_start ();
+ acc = gnm_accumulator_new ();
+
+ for (r = 0; r < C->rows; r++) {
+ for (c = 0; c < C->cols; c++) {
go_accumulator_clear (acc);
- for (i = 0; i < cols_a; ++i) {
+ for (i = 0; i < A->cols; ++i) {
GnmQuad p;
gnm_quad_mul12 (&p,
- A[r + i * rows_a],
- B[i + c * cols_a]);
+ A->data[r][i],
+ B->data[i][c]);
gnm_accumulator_add_quad (acc, &p);
}
- product[r + c * rows_a] = gnm_accumulator_value (acc);
+ C->data[r][c] = gnm_accumulator_value (acc);
}
}
+
gnm_accumulator_free (acc);
gnm_accumulator_end (state);
}
diff --git a/src/mathfunc.h b/src/mathfunc.h
index d306a7c..075d926 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -146,10 +146,20 @@ gnm_float discpfuncinverter (gnm_float p, const gnm_float shape[],
GnmPFunc pfunc);
/* ------------------------------------------------------------------------- */
-
/* Matrix functions. */
-void mmult (gnm_float *A, gnm_float *B, int cols_a, int rows_a, int cols_b,
- gnm_float *product);
+
+typedef struct {
+ gnm_float **data; /* [y][x] */
+ int cols, rows;
+} GnmMatrix;
+
+GnmMatrix *gnm_matrix_new (int rows, int cols); /* Note the order: y then x. */
+void gnm_matrix_free (GnmMatrix *m);
+GnmMatrix *gnm_matrix_from_value (GnmValue const *v, GnmValue **perr, GnmEvalPos const *ep);
+GnmValue *gnm_matrix_to_value (GnmMatrix const *m);
+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);
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]