[goffice] Math: extract and formalize matrix stuff.



commit 5168578d9ff0229961f7327669acc0550b5c4401
Author: Morten Welinder <terra gnome org>
Date:   Sun May 5 17:01:40 2013 -0400

    Math: extract and formalize matrix stuff.
    
    This is probably not the final word.  In particular, we need a better
    interface for all the non-quad stuff.

 goffice/Makefile.am          |    2 +
 goffice/math/go-matrix.c     |  727 +++++++++++++++++++++++++++++++++++++
 goffice/math/go-matrix.h     |   95 +++++
 goffice/math/go-regression.c |  820 +++++++-----------------------------------
 goffice/math/goffice-math.h  |    5 +
 5 files changed, 958 insertions(+), 691 deletions(-)
---
diff --git a/goffice/Makefile.am b/goffice/Makefile.am
index 98bc860..8acde66 100644
--- a/goffice/Makefile.am
+++ b/goffice/Makefile.am
@@ -340,6 +340,7 @@ math_SOURCES =                                      \
        math/go-cspline.c                       \
        math/go-complex.c                       \
        math/go-fft.c                           \
+       math/go-matrix.c                        \
        math/go-matrix3x3.c                     \
        math/go-quad.c                          \
        math/go-R.c                             \
@@ -355,6 +356,7 @@ math_HEADERS =                                      \
        math/go-cspline.h                       \
        math/go-complex.h                       \
        math/go-fft.h                           \
+       math/go-matrix.h                        \
        math/go-matrix3x3.h                     \
        math/go-quad.h                          \
        math/go-R.h                             \
diff --git a/goffice/math/go-matrix.c b/goffice/math/go-matrix.c
new file mode 100644
index 0000000..7b4769c
--- /dev/null
+++ b/goffice/math/go-matrix.c
@@ -0,0 +1,727 @@
+/*
+ * go-matrix.c: Matrix routines.
+ *
+ * Authors:
+ *   Morten Welinder <terra gnome org>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) version 3.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
+ * USA.
+ *
+ */
+
+#include <goffice/goffice-config.h>
+#include <goffice/goffice.h>
+#include <math.h>
+
+
+#ifndef DOUBLE
+
+#define QUAD SUFFIX(GOQuad)
+#define QQR SUFFIX(GOQuadQR)
+#define QMATRIX SUFFIX(GOQuadMatrix)
+
+#define DOUBLE double
+#define SUFFIX(_n) _n
+
+struct GOQuadQR_ {
+       QMATRIX *V;
+       QMATRIX *R;
+       QUAD det;
+};
+
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+#include "go-matrix.c"
+#undef DOUBLE
+#undef SUFFIX
+#define DOUBLE long double
+#define SUFFIX(_n) _n ## l
+
+struct GOQuadQRl_ {
+       QMATRIX *V;
+       QMATRIX *R;
+       QUAD det;
+};
+
+#endif
+
+#endif
+
+
+/**
+ * go_quad_matrix_new: (skip)
+ * @m: number of rows
+ * @n: number of columns
+ *
+ * Returns: a new zero matrix.
+ **/
+/**
+ * go_quad_matrix_newl: (skip)
+ **/
+QMATRIX *
+SUFFIX(go_quad_matrix_new) (int m, int n)
+{
+       QMATRIX *res;
+       int i;
+
+       g_return_val_if_fail (m >= 1, NULL);
+       g_return_val_if_fail (n >= 1, NULL);
+
+       res = g_new (QMATRIX, 1);
+       res->m = m;
+       res->n = n;
+       res->data = g_new (QUAD *, m);
+
+       for (i = 0; i < m; i++)
+               res->data[i] = g_new0 (QUAD, n);
+
+       return res;
+}
+
+void
+SUFFIX(go_quad_matrix_free) (QMATRIX *A)
+{
+       int i;
+
+       for (i = 0; i < A->m; i++)
+               g_free (A->data[i]);
+       g_free (A->data);
+       g_free (A);
+}
+
+
+/**
+ * go_quad_matrix_dup: (skip)
+ **/
+/**
+ * go_quad_matrix_dupl: (skip)
+ **/
+QMATRIX *
+SUFFIX(go_quad_matrix_dup) (const QMATRIX *A)
+{
+       QMATRIX *res;
+
+       g_return_val_if_fail (A != NULL, NULL);
+       res = SUFFIX(go_quad_matrix_new) (A->m, A->n);
+       SUFFIX(go_quad_matrix_copy) (res, A);
+       return res;
+}
+
+/**
+ * go_quad_matrix_copy:
+ * @A: Destination matrix.
+ * @B: Source matrix.
+ *
+ * Copies B to A.
+ */
+void
+SUFFIX(go_quad_matrix_copy) (QMATRIX *A, const QMATRIX *B)
+{
+       int i, j;
+
+       g_return_if_fail (A != NULL);
+       g_return_if_fail (B != NULL);
+       g_return_if_fail (A->m == B->m && A->n == B->n);
+
+       if (A == B)
+               return;
+
+       for (i = 0; i < A->m; i++)
+               for (j = 0; j < A->n; j++)
+                       A->data[i][j] = B->data[i][j];
+}
+
+/**
+ * go_quad_matrix_transpose:
+ * @A: Destination matrix.
+ * @B: Source matrix.
+ *
+ * Transposes B into A.
+ */
+void
+SUFFIX(go_quad_matrix_transpose) (QMATRIX *A, const QMATRIX *B)
+{
+       int i, j;
+
+       g_return_if_fail (A != NULL);
+       g_return_if_fail (B != NULL);
+       g_return_if_fail (A->m == B->n && A->n == B->m);
+
+       if (A == B)
+               return;
+
+       for (i = 0; i < A->m; i++)
+               for (j = 0; j < A->n; j++)
+                       A->data[i][j] = B->data[j][i];
+}
+
+
+/**
+ * go_quad_matrix_multiply:
+ * @C: Destination matrix.
+ * @A: Source matrix.
+ * @B: Source matrix.
+ *
+ * Multiplies A*B and stores the result in C.
+ */
+void
+SUFFIX(go_quad_matrix_multiply) (QMATRIX *C,
+                                const QMATRIX *A,
+                                const QMATRIX *B)
+{
+       int i, j, k;
+
+       g_return_if_fail (C != NULL);
+       g_return_if_fail (A != NULL);
+       g_return_if_fail (B != NULL);
+       g_return_if_fail (C->m == A->m && A->n == B->m && B->n == C->n);
+       g_return_if_fail (C != A && C != B);
+
+       for (i = 0; i < C->m; i++) {
+               for (j = 0; j < C->n; j++) {
+                       QUAD p, acc = SUFFIX(go_quad_zero);
+                       for (k = 0; k < A->n; k++) {
+                               SUFFIX(go_quad_mul) (&p,
+                                                    &A->data[i][k],
+                                                    &B->data[k][j]);
+                               SUFFIX(go_quad_add) (&acc, &acc, &p);
+                       }
+                       C->data[i][j] = acc;
+               }
+       }
+}
+
+/**
+ * go_quad_matrix_inverse: (skip)
+ * @A: Source matrix.
+ * @threshold: condition number threshold.
+ *
+ * Returns: The inverse matrix of A.  If any eigenvalues divided by the largest
+ * eigenvalue is less than or equal to the given threshold, %NULL is returning
+ * indicating a matrix that cannot be inverted.  (Note: this doesn't actually
+ * use the eigenvalues of A, but of A after an orthogonal transformation.)
+ */
+/**
+ * go_quad_matrix_inversel: (skip)
+ */
+QMATRIX *
+SUFFIX(go_quad_matrix_inverse) (const QMATRIX *A, DOUBLE threshold)
+{
+       QQR *qr;
+       int i, k, n;
+       QMATRIX *Z;
+       const QMATRIX *R;
+       gboolean ok;
+       QUAD *x, *QTk;
+       DOUBLE emin, emax;
+
+       g_return_val_if_fail (A != NULL, NULL);
+       g_return_val_if_fail (A->m == A->n, NULL);
+       g_return_val_if_fail (threshold >= 0, NULL);
+
+       qr = SUFFIX(go_quad_qr_new) (A);
+       if (!qr)
+               return NULL;
+
+       n = A->n;
+       Z = SUFFIX(go_quad_matrix_new) (n, n);
+       x = g_new (QUAD, n);
+       QTk = g_new (QUAD, n);
+
+       R = SUFFIX(go_quad_qr_r) (qr);
+       SUFFIX(go_quad_matrix_eigen_range) (R, &emin, &emax);
+       ok = (emin > emax * threshold);
+
+       for (k = 0; ok && k < n; k++) {
+               /* Compute Q^T's k-th column.  */
+               for (i = 0; i < n; i++)
+                       SUFFIX(go_quad_init)(&QTk[i], i == k);
+               SUFFIX(go_quad_qr_multiply_qt) (qr, QTk);
+
+               /* Solve R x = Q^T e_k */
+               if (SUFFIX(go_quad_matrix_back_solve) (R, x, QTk)) {
+                       ok = FALSE;
+                       break;
+               }
+
+               for (i = 0; i < n; i++)
+                       Z->data[i][k] = x[i];
+       }
+
+       SUFFIX(go_quad_qr_free) (qr);
+       g_free (QTk);
+       g_free (x);
+
+       if (!ok) {
+               SUFFIX(go_quad_matrix_free) (Z);
+               return NULL;
+       }
+
+       return Z;
+}
+
+void
+SUFFIX(go_quad_matrix_determinant) (const QMATRIX *A, QUAD *res)
+{
+       QQR *qr;
+
+       g_return_if_fail (A != NULL);
+       g_return_if_fail (A->m == A->n);
+       g_return_if_fail (res != NULL);
+
+       if (A->m == 1) {
+               *res = A->data[0][0];
+               return;
+       }
+
+       if (A->m == 2) {
+               QUAD a, b;
+               SUFFIX(go_quad_mul)(&a, &A->data[0][0], &A->data[1][1]);
+               SUFFIX(go_quad_mul)(&b, &A->data[1][0], &A->data[0][1]);
+               SUFFIX(go_quad_sub)(res, &a, &b);
+               return;
+       }
+
+       qr = SUFFIX(go_quad_qr_new) (A);
+       if (!qr) {
+               /* Hmm... */
+               SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
+               return;
+       }
+
+       SUFFIX(go_quad_qr_determinant) (qr, res);
+       SUFFIX(go_quad_qr_free) (qr);
+}
+
+/**
+ * go_quad_matrix_pseudo_inverse: (skip)
+ * @A: An arbitrary matrix.
+ * @threshold: condition number threshold.
+ */
+/**
+ * go_quad_matrix_pseudo_inversel: (skip)
+ */
+QMATRIX *
+SUFFIX(go_quad_matrix_pseudo_inverse) (const QMATRIX *A, DOUBLE threshold)
+{
+       int i, j, m, n;
+       QQR *qr;
+       const QMATRIX *R;
+       QMATRIX *RT;
+       QMATRIX *RTR;
+       QMATRIX *B0;
+       QMATRIX *Bi;
+       QMATRIX *B;
+       DOUBLE emax;
+       gboolean full_rank;
+       QUAD delta;
+       int steps;
+       QUAD *x;
+
+       g_return_val_if_fail (A != NULL, NULL);
+       g_return_val_if_fail (threshold >= 0, NULL);
+
+       m = A->m;
+       n = A->n;
+       B = SUFFIX(go_quad_matrix_new) (n, m);
+
+       if (m < n) {
+               /*
+                * The main code assumes m >= n.  Luckily, taking the
+                * pseudo-inverse commutes with transposition.
+                */
+
+               QMATRIX *AT = SUFFIX(go_quad_matrix_new) (n, m);
+               QMATRIX *BT;
+
+               SUFFIX(go_quad_matrix_transpose) (AT, A);
+               BT = SUFFIX(go_quad_matrix_pseudo_inverse) (AT, threshold);
+               SUFFIX(go_quad_matrix_transpose) (B, BT);
+
+               SUFFIX (go_quad_matrix_free) (AT);
+               SUFFIX (go_quad_matrix_free) (BT);
+               return B;
+       }
+
+       qr = SUFFIX(go_quad_qr_new) (A);
+       if (!qr)
+               goto out;
+       R = SUFFIX(go_quad_qr_r) (qr);
+
+       SUFFIX(go_quad_matrix_eigen_range) (R, NULL, &emax);
+       if (emax == 0)
+               goto out;
+
+       full_rank = TRUE;
+       for (i = 0; i < n; i++) {
+               DOUBLE abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value)(&R->data[i][i]));
+               if (abs_e <= emax * threshold) {
+                       full_rank = FALSE;
+                       R->data[i][i] = SUFFIX(go_quad_zero);
+               }
+       }
+
+       SUFFIX(go_quad_init) (&delta, full_rank ? 0 : emax * threshold);
+
+       /*
+        * Starting point for the iteration:
+        *
+        *    Bi := (RT R + delta I)^-1 * RT
+        *
+        * (RT R) is positive semi-definite and (delta I) is positive definite,
+        * so the sum is positive definite and invertible.
+        */
+       RT = SUFFIX(go_quad_matrix_new) (n, n);
+       SUFFIX(go_quad_matrix_transpose) (RT, R);
+       RTR = SUFFIX(go_quad_matrix_new) (n, n);
+       SUFFIX(go_quad_matrix_multiply) (RTR, RT, R);
+       for (i = 0; i < n; i++)
+               SUFFIX(go_quad_add) (&RTR->data[i][i],
+                                    &RTR->data[i][i],
+                                    &delta);
+       B0 = SUFFIX(go_quad_matrix_inverse) (RTR, 0.0);
+       Bi = SUFFIX(go_quad_matrix_new) (n, n);
+       SUFFIX(go_quad_matrix_multiply) (Bi, B0, RT);
+       SUFFIX(go_quad_matrix_free) (B0);
+       SUFFIX(go_quad_matrix_free) (RTR);
+       SUFFIX(go_quad_matrix_free) (RT);
+
+       /* B_{i+1} = (2 I - B_i R) B_i */
+       for (steps = 0; steps < 10; steps++) {
+               QMATRIX *W = SUFFIX(go_quad_matrix_new) (n, n);
+               QMATRIX *Bip1 = SUFFIX(go_quad_matrix_new) (n, n);
+               QUAD two;
+
+               SUFFIX(go_quad_init)(&two, 2);
+
+               SUFFIX(go_quad_matrix_multiply) (W, Bi, R);
+               for (i = 0; i < n; i++)
+                       for (j = 0; j < n; j++)
+                               SUFFIX(go_quad_sub) (&W->data[i][j],
+                                                    i == j ? &two : &SUFFIX(go_quad_zero),
+                                                    &W->data[i][j]);
+               SUFFIX(go_quad_matrix_multiply) (Bip1, W, Bi);
+               SUFFIX(go_quad_matrix_copy) (Bi, Bip1);
+
+               SUFFIX(go_quad_matrix_free) (Bip1);
+               SUFFIX(go_quad_matrix_free) (W);
+       }
+
+       /* B := (Bi|O) Q^T */
+       x = g_new (QUAD, m);
+       for (j = 0; j < m; j++) {
+               /* Compute Q^T e_j.  */
+               for (i = 0; i < m; i++)
+                       SUFFIX(go_quad_init)(&x[i], i == j ? 1 : 0);
+               SUFFIX(go_quad_qr_multiply_qt) (qr, x);
+
+               for (i = 0; i < n; i++) {
+                       int k;
+
+                       B->data[i][j] = SUFFIX(go_quad_zero);
+                       for (k = 0; k < n /* Only n */; k++) {
+                               QUAD p;
+                               SUFFIX(go_quad_mul) (&p, &Bi->data[i][k], &x[k]);
+                               SUFFIX(go_quad_add) (&B->data[i][j], &B->data[i][j], &p);
+                       }
+               }
+       }
+       g_free (x);
+       SUFFIX(go_quad_matrix_free) (Bi);
+
+out:
+       SUFFIX(go_quad_qr_free) (qr);
+
+       return B;
+}
+
+/**
+ * go_quad_matrix_fwd_solve:
+ * @R: An upper triangular matrix.
+ * @x: (out): Result vector.
+ * @b: Input vector.
+ *
+ * Returns: %TRUE on error.
+ *
+ * This function solves the triangular system RT*x=b.
+ */
+gboolean
+SUFFIX(go_quad_matrix_fwd_solve) (const QMATRIX *R, QUAD *x, const QUAD *b)
+{
+       int i, j, n;
+
+       g_return_val_if_fail (R != NULL, TRUE);
+       g_return_val_if_fail (R->m == R-> n, TRUE);
+       g_return_val_if_fail (x != NULL, TRUE);
+       g_return_val_if_fail (b != NULL, TRUE);
+
+       n = R->m;
+
+       for (i = 0; i < n; i++) {
+               QUAD d = b[i];
+               for (j = 0; j < i; j++) {
+                       QUAD p;
+                       SUFFIX(go_quad_mul) (&p, &R->data[j][i], &x[j]);
+                       SUFFIX(go_quad_sub) (&d, &d, &p);
+               }
+               if (SUFFIX(go_quad_value)(&R->data[i][i]) == 0) {
+                       while (i < n) x[i++] = SUFFIX(go_quad_zero);
+                       return TRUE;
+               }
+               SUFFIX(go_quad_div) (&x[i], &d, &R->data[i][i]);
+       }
+
+       return FALSE;
+}
+
+/**
+ * go_quad_matrix_back_solve:
+ * @R: An upper triangular matrix.
+ * @x: (out): Result vector.
+ * @b: Input vector.
+ *
+ * Returns: %TRUE on error.
+ *
+ * This function solves the triangular system R*x=b.
+ */
+gboolean
+SUFFIX(go_quad_matrix_back_solve) (const QMATRIX *R, QUAD *x, const QUAD *b)
+{
+       int i, j, n;
+
+       g_return_val_if_fail (R != NULL, TRUE);
+       g_return_val_if_fail (R->m == R-> n, TRUE);
+       g_return_val_if_fail (x != NULL, TRUE);
+       g_return_val_if_fail (b != NULL, TRUE);
+
+       n = R->m;
+
+       for (i = n - 1; i >= 0; i--) {
+               QUAD d = b[i];
+               for (j = i + 1; j < n; j++) {
+                       QUAD p;
+                       SUFFIX(go_quad_mul) (&p, &R->data[i][j], &x[j]);
+                       SUFFIX(go_quad_sub) (&d, &d, &p);
+               }
+               if (SUFFIX(go_quad_value)(&R->data[i][i]) == 0) {
+                       while (i >= 0)
+                               x[i--] = SUFFIX(go_quad_zero);
+                       return TRUE;
+               }
+               SUFFIX(go_quad_div) (&x[i], &d, &R->data[i][i]);
+       }
+
+       return FALSE;
+}
+
+/**
+ * go_quad_matrix_eigen_range:
+ * @A: Triangular matrix.
+ * @emin: (out): Smallest absolute eigen value.
+ * @emax: (out): Largest absolute eigen value.
+ */
+void
+SUFFIX(go_quad_matrix_eigen_range) (const QMATRIX *A,
+                                   DOUBLE *emin, DOUBLE *emax)
+{
+       int i;
+       DOUBLE abs_e;
+
+       g_return_if_fail (A != NULL);
+       g_return_if_fail (A->m == A->n);
+
+       abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value) (&A->data[0][0]));
+       if (emin) *emin = abs_e;
+       if (emax) *emax = abs_e;
+       for (i = 1; i < A->m; i++) {
+               abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value) (&A->data[i][i]));
+               if (emin) *emin = MIN (abs_e, *emin);
+               if (emax) *emax = MAX (abs_e, *emax);
+       }
+}
+
+
+/* -------------------------------------------------------------------------- */
+
+/**
+ * go_quad_qr_new: (skip)
+ * @A: Source matrix.
+ *
+ * QR decomposition of a matrix using Householder matrices.
+ *
+ * A (input) is an m-times-n matrix.  A[0...m-1][0..n-1]
+ * If qAT is TRUE, this parameter is transposed.
+ *
+ * V is a pre-allocated output m-times-n matrix.  V will contrain
+ * n vectors of different lengths: n, n-1, ..., 1.  These are the
+ * Householder vectors (or null for the degenerate case).  The
+ * matrix Q of size m-times-m is implied from V.
+ *
+ * R is a matrix of size n-times-n.  (To get the m-times-n version
+ * of R, simply add m-n null rows.)
+ */
+
+/**
+ * go_quad_qr_newl: (skip)
+ */
+QQR *
+SUFFIX(go_quad_qr_new) (const QMATRIX *A)
+{
+       QQR *qr;
+       int qdet = 1;
+       QMATRIX *R;
+       QMATRIX *V;
+       int i, j, k, m, n;
+       QUAD *tmp;
+
+       g_return_val_if_fail (A != NULL, NULL);
+       g_return_val_if_fail (A->m >= A->n, NULL);
+
+       m = A->m;
+       n = A->n;
+
+       qr = g_new (QQR, 1);
+       V = qr->V = SUFFIX(go_quad_matrix_new) (m, n);
+       qr->R = SUFFIX(go_quad_matrix_new) (n, n);
+
+       /* Temporary m-by-n version of R.  */
+       R = SUFFIX(go_quad_matrix_dup) (A);
+
+       tmp = g_new (QUAD, n);
+
+       for (k = 0; k < n; k++) {
+               QUAD L, L2 = SUFFIX(go_quad_zero), L2p = L2, s;
+
+               for (i = m - 1; i >= k; i--) {
+                       V->data[i][k] = R->data[i][k];
+                       SUFFIX(go_quad_mul)(&s, &V->data[i][k], &V->data[i][k]);
+                       L2p = L2;
+                       SUFFIX(go_quad_add)(&L2, &L2, &s);
+               }
+               SUFFIX(go_quad_sqrt)(&L, &L2);
+
+               (SUFFIX(go_quad_value)(&V->data[k][k]) < 0
+                ? SUFFIX(go_quad_sub)
+                : SUFFIX(go_quad_add)) (&V->data[k][k], &V->data[k][k], &L);
+
+               /* Normalize v[k] to length 1.  */
+               SUFFIX(go_quad_mul)(&s, &V->data[k][k], &V->data[k][k]);
+               SUFFIX(go_quad_add)(&L2p, &L2p, &s);
+               SUFFIX(go_quad_sqrt)(&L, &L2p);
+               if (SUFFIX(go_quad_value)(&L) == 0) {
+                       /* This will be an identity so no determinant sign */
+                       continue;
+               }
+               for (i = k; i < m; i++)
+                       SUFFIX(go_quad_div)(&V->data[i][k], &V->data[i][k], &L);
+
+               /* Householder matrices have determinant -1.  */
+               qdet = -qdet;
+
+               /* Calculate tmp = v[k]^t * R[k:m,k:n] */
+               for (j = k; j < n; j++) {
+                       tmp[j] = SUFFIX(go_quad_zero);
+                       for (i = k ; i < m; i++) {
+                               QUAD p;
+                               SUFFIX(go_quad_mul) (&p, &V->data[i][k], &R->data[i][j]);
+                               SUFFIX(go_quad_add) (&tmp[j], &tmp[j], &p);
+                       }
+               }
+
+               /* R[k:m,k:n] -= v[k] * tmp */
+               for (j = k; j < n; j++) {
+                       for (i = k; i < m; i++) {
+                               QUAD p;
+                               SUFFIX(go_quad_mul) (&p, &V->data[i][k], &tmp[j]);
+                               SUFFIX(go_quad_add) (&p, &p, &p);
+                               SUFFIX(go_quad_sub) (&R->data[i][j], &R->data[i][j], &p);
+                       }
+               }
+
+               /* Explicitly zero what should become zero.  */
+               for (i = k + 1; i < m; i++)
+                       R->data[i][k] = SUFFIX(go_quad_zero);
+       }
+
+       g_free (tmp);
+
+       for (i = 0; i < n /* Only n */; i++)
+               for (j = 0; j < n; j++)
+                       qr->R->data[i][j] = R->data[i][j];
+
+       SUFFIX(go_quad_init) (&qr->det, qdet);
+       for (i = 0; i < n; i++)
+               SUFFIX(go_quad_mul) (&qr->det, &qr->det, &R->data[i][i]);
+
+       SUFFIX(go_quad_matrix_free) (R);
+
+       return qr;
+}
+
+void
+SUFFIX(go_quad_qr_free) (QQR *qr)
+{
+       g_return_if_fail (qr != NULL);
+
+       SUFFIX(go_quad_matrix_free) (qr->V);
+       SUFFIX(go_quad_matrix_free) (qr->R);
+       g_free (qr);
+}
+
+void
+SUFFIX(go_quad_qr_determinant) (const QQR *qr, QUAD *det)
+{
+       g_return_if_fail (qr != NULL);
+       g_return_if_fail (det != NULL);
+
+       *det = qr->det;
+}
+
+const QMATRIX *
+SUFFIX(go_quad_qr_r) (const QQR *qr)
+{
+       g_return_val_if_fail (qr != NULL, NULL);
+
+       return qr->R;
+}
+
+/**
+ * go_quad_qr_multiply_qt:
+ * @qr: A QR decomposition.
+ * @x: (inout): a vector.
+ *
+ * Replaces @x by Q^t * x
+ */
+void
+SUFFIX(go_quad_qr_multiply_qt) (const QQR *qr, QUAD *x)
+{
+       int i, k;
+       QMATRIX *V = qr->V;
+
+       for (k = 0; k < V->n; k++) {
+               QUAD s = SUFFIX(go_quad_zero);
+               for (i = k; i < V->m; i++) {
+                       QUAD p;
+                       SUFFIX(go_quad_mul) (&p, &x[i], &V->data[i][k]);
+                       SUFFIX(go_quad_add) (&s, &s, &p);
+               }
+               SUFFIX(go_quad_add) (&s, &s, &s);
+               for (i = k; i < V->m; i++) {
+                       QUAD p;
+                       SUFFIX(go_quad_mul) (&p, &s, &V->data[i][k]);
+                       SUFFIX(go_quad_sub) (&x[i], &x[i], &p);
+               }
+       }
+}
diff --git a/goffice/math/go-matrix.h b/goffice/math/go-matrix.h
new file mode 100644
index 0000000..b27aa27
--- /dev/null
+++ b/goffice/math/go-matrix.h
@@ -0,0 +1,95 @@
+#ifndef GO_MATRIX_H__
+#define GO_MATRIX_H__
+
+#include <glib.h>
+
+G_BEGIN_DECLS
+
+struct GOQuadMatrix_ {
+       GOQuad **data;   /* [m][n] */
+       int m; /* down */
+       int n; /* right */
+};
+
+GOQuadMatrix *go_quad_matrix_new (int m, int n);
+GOQuadMatrix *go_quad_matrix_dup (const GOQuadMatrix *A);
+void go_quad_matrix_free (GOQuadMatrix *A);
+
+void go_quad_matrix_copy (GOQuadMatrix *A, const GOQuadMatrix *B);
+void go_quad_matrix_transpose (GOQuadMatrix *A, const GOQuadMatrix *B);
+
+void go_quad_matrix_multiply (GOQuadMatrix *C,
+                             const GOQuadMatrix *A,
+                             const GOQuadMatrix *B);
+
+GOQuadMatrix *go_quad_matrix_inverse (const GOQuadMatrix *A, double threshold);
+
+void go_quad_matrix_determinant (const GOQuadMatrix *A, GOQuad *res);
+
+GOQuadMatrix *go_quad_matrix_pseudo_inverse (const GOQuadMatrix *A,
+                                            double threshold);
+
+gboolean go_quad_matrix_back_solve (const GOQuadMatrix *R, GOQuad *x,
+                                   const GOQuad *b);
+gboolean go_quad_matrix_fwd_solve (const GOQuadMatrix *R, GOQuad *x,
+                                  const GOQuad *b);
+
+void go_quad_matrix_eigen_range (const GOQuadMatrix *A,
+                                double *emin, double *emax);
+
+/* ---------------------------------------- */
+
+GOQuadQR *go_quad_qr_new (const GOQuadMatrix *A);
+void go_quad_qr_free (GOQuadQR *qr);
+void go_quad_qr_determinant (const GOQuadQR *qr, GOQuad *det);
+const GOQuadMatrix *go_quad_qr_r (const GOQuadQR *qr);
+void go_quad_qr_multiply_qt (const GOQuadQR *qr, GOQuad *x);
+
+/* ---------------------------------------- */
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+struct GOQuadMatrixl_ {
+       GOQuadl **data;   /* [m][n] */
+       int m; /* down */
+       int n; /* right */
+};
+
+GOQuadMatrixl *go_quad_matrix_newl (int m, int n);
+GOQuadMatrixl *go_quad_matrix_dupl (const GOQuadMatrixl *A);
+void go_quad_matrix_freel (GOQuadMatrixl *A);
+
+void go_quad_matrix_copyl (GOQuadMatrixl *A, const GOQuadMatrixl *B);
+void go_quad_matrix_transposel (GOQuadMatrixl *A, const GOQuadMatrixl *B);
+
+void go_quad_matrix_multiplyl (GOQuadMatrixl *C,
+                              const GOQuadMatrixl *A,
+                              const GOQuadMatrixl *B);
+
+GOQuadMatrixl *go_quad_matrix_inversel (const GOQuadMatrixl *A, long double threshold);
+
+void go_quad_matrix_determinantl (const GOQuadMatrixl *A, GOQuadl *res);
+
+GOQuadMatrixl *go_quad_matrix_pseudo_inversel (const GOQuadMatrixl *A,
+                                              long double threshold);
+
+gboolean go_quad_matrix_back_solvel (const GOQuadMatrixl *R, GOQuadl *x,
+                                    const GOQuadl *b);
+gboolean go_quad_matrix_fwd_solvel (const GOQuadMatrixl *R, GOQuadl *x,
+                                   const GOQuadl *b);
+
+void go_quad_matrix_eigen_rangel (const GOQuadMatrixl *A,
+                                 long double *emin, long double *emax);
+
+/* ---------------------------------------- */
+
+GOQuadQRl *go_quad_qr_newl (const GOQuadMatrixl *A);
+void go_quad_qr_freel (GOQuadQRl *qr);
+void go_quad_qr_determinantl (const GOQuadQRl *qr, GOQuadl *det);
+const GOQuadMatrixl *go_quad_qr_rl (const GOQuadQRl *qr);
+void go_quad_qr_multiply_qtl (const GOQuadQRl *qr, GOQuadl *x);
+
+#endif
+
+G_END_DECLS
+
+#endif
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index dec0e36..f207fe7 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -90,9 +90,12 @@
 
 #ifndef DOUBLE
 
+#define DEFAULT_THRESHOLD (256 * DOUBLE_EPS)
+
+
 #define DEFINE_COMMON
 #define DOUBLE double
-#define DOUBLE_MANT_DIG DBL_MANT_DIG
+#define DOUBLE_EPS DBL_EPSILON
 #define SUFFIX(_n) _n
 #define FORMAT_f "f"
 #define FORMAT_g "g"
@@ -116,7 +119,7 @@ general_linear_regressionl (long double *const *const xssT, int n,
 #undef BOUNCE
 #undef DEFINE_COMMON
 #undef DOUBLE
-#undef DOUBLE_MANT_DIG
+#undef DOUBLE_EPS
 #undef SUFFIX
 #undef FORMAT_f
 #undef FORMAT_g
@@ -124,7 +127,7 @@ general_linear_regressionl (long double *const *const xssT, int n,
 #include <sunmath.h>
 #endif
 #define DOUBLE long double
-#define DOUBLE_MANT_DIG LDBL_MANT_DIG
+#define DOUBLE_EPS LDBL_EPSILON
 #define SUFFIX(_n) _n ## l
 #define FORMAT_f "Lf"
 #define FORMAT_g "Lg"
@@ -264,6 +267,31 @@ go_regression_statl_get_type (void)
  *
  */
 
+static SUFFIX(GOQuadMatrix)
+*
+SUFFIX(quad_matrix_from_matrix) (CONSTMATRIX A, int m, int n)
+{
+       int i, j;
+       SUFFIX(GOQuadMatrix) *qA = SUFFIX(go_quad_matrix_new) (m, n);
+
+       for (i = 0; i < m; i++)
+               for (j = 0; j < n; j++)
+                       SUFFIX(go_quad_init) (&qA->data[i][j], A[i][j]);
+
+       return qA;
+}
+
+static void
+SUFFIX(copy_quad_matrix_to_matrix) (MATRIX A, const SUFFIX(GOQuadMatrix) *qA)
+{
+       int i, j;
+
+       for (i = 0; i < qA->m; i++)
+               for (j = 0; j < qA->n; j++)
+                       A[i][j] = SUFFIX(go_quad_value) (&qA->data[i][j]);
+
+}
+
 /* ------------------------------------------------------------------------- */
 #ifdef BOUNCE
 /* Suppport routines for bouncing double to long double.  */
@@ -312,159 +340,7 @@ copy_stats (go_regression_stat_t *s1,
 
 /* ------------------------------------------------------------------------- */
 
-/*
- * QR decomposition of a matrix using Householder matrices.
- *
- * A (input) is an m-times-n matrix.  A[0...m-1][0..n-1]
- * If qAT is TRUE, this parameter is transposed.
- *
- * V is a pre-allocated output m-times-n matrix.  V will contrain
- * n vectors of different lengths: n, n-1, ..., 1.  These are the
- * Householder vectors (or null for the degenerate case).  The
- * matrix Q of size m-times-m is implied from V.
- *
- * Rfinal is a pre-allocated output matrix of size n-times-n.  (To
- * get the m-times-n version of R, simply add m-n null rows.)
- *
- * qdet is an optional output location for det(Q).
- *
- * Note: the output is V and R, not Q and R.
- */
-static gboolean
-SUFFIX(QRH) (CONSTMATRIX A, gboolean qAT,
-            QMATRIX V, QMATRIX Rfinal, int m, int n,
-            int *qdet)
-{
-       QMATRIX R = NULL;
-       int i, j, k;
-       QUAD *tmp = g_new (QUAD, n);
-       gboolean ok = TRUE;
-       gboolean debug = FALSE;
-
-       if (m < n || n < 1) {
-               ok = FALSE;
-               goto out;
-       }
-
-       ALLOC_MATRIX2 (R, m, n, QUAD);
-       for (i = 0; i < m; i++) {
-               for (j = 0; j < n; j++) {
-                       DOUBLE Aij = qAT ? A[j][i] : A[i][j];
-                       if (!SUFFIX(go_finite)(Aij)) {
-                               ok = FALSE;
-                               goto out;
-                       }
-                       SUFFIX(go_quad_init) (&R[i][j], Aij);
-                       SUFFIX(go_quad_init) (&V[i][j], -42);
-               }
-       }
-
-       if (qdet)
-               *qdet = 1;
-
-       for (k = 0; k < n; k++) {
-               QUAD L, L2, L2p, s;
-
-               SUFFIX(go_quad_init) (&L2, 0);
-               L2p = L2;
-               for (i = m - 1; i >= k; i--) {
-                       V[i][k] = R[i][k];
-                       SUFFIX(go_quad_mul)(&s, &V[i][k], &V[i][k]);
-                       L2p = L2;
-                       SUFFIX(go_quad_add)(&L2, &L2, &s);
-               }
-               SUFFIX(go_quad_sqrt)(&L, &L2);
-
-               (SUFFIX(go_quad_value)(&V[k][k]) < 0
-                ? SUFFIX(go_quad_sub)
-                : SUFFIX(go_quad_add)) (&V[k][k], &V[k][k], &L);
-
-               /* Normalize v[k] to length 1.  */
-               SUFFIX(go_quad_mul)(&s, &V[k][k], &V[k][k]);
-               SUFFIX(go_quad_add)(&L2p, &L2p, &s);
-               SUFFIX(go_quad_sqrt)(&L, &L2p);
-               if (SUFFIX(go_quad_value)(&L) == 0) {
-                       if (debug)
-                               g_printerr ("L=0 in round %d:\n", k);
-                       /* This will be an identity so no determinant sign */
-                       continue;
-               }
-               for (i = k; i < m; i++)
-                       SUFFIX(go_quad_div)(&V[i][k], &V[i][k], &L);
-
-               /* Householder matrices have determinant -1.  */
-               if (qdet)
-                       *qdet = -*qdet;
-
-               /* Calculate tmp = v[k]^t * R[k:m,k:n] */
-               for (j = k; j < n; j++) {
-                       SUFFIX(go_quad_init) (&tmp[j], 0);
-                       for (i = k ; i < m; i++) {
-                               QUAD p;
-                               SUFFIX(go_quad_mul) (&p, &V[i][k], &R[i][j]);
-                               SUFFIX(go_quad_add) (&tmp[j], &tmp[j], &p);
-                       }
-               }
-
-               /* R[k:m,k:n] -= v[k] * tmp */
-               for (j = k; j < n; j++) {
-                       for (i = k; i < m; i++) {
-                               QUAD p;
-                               SUFFIX(go_quad_mul) (&p, &V[i][k], &tmp[j]);
-                               SUFFIX(go_quad_add) (&p, &p, &p);
-                               SUFFIX(go_quad_sub) (&R[i][j], &R[i][j], &p);
-                       }
-               }
-
-               /* Explicitly zero what should become zero.  */
-               for (i = k + 1; i < m; i++)
-                       SUFFIX(go_quad_init) (&R[i][k], 0);
-
-               if (debug) {
-                       g_printerr ("After round %d:\n", k);
-                       g_printerr ("R:\n");
-                       PRINT_QMATRIX(R, m, n);
-                       g_printerr ("\n");
-                       g_printerr ("V:\n");
-                       PRINT_QMATRIX(V, m, n);
-                       g_printerr ("\n\n");
-               }
-       }
-
-       for (i = 0; i < n /* Only n */; i++)
-               for (j = 0; j < n; j++)
-                       Rfinal[i][j] = R[i][j];
-
-out:
-       if (R)
-               FREE_MATRIX (R, m, n);
-       g_free (tmp);
-
-       return ok;
-}
-
-static void
-SUFFIX(multiply_Qt) (QUAD *x, CONSTQMATRIX V, int m, int n)
-{
-       int i, k;
-
-       for (k = 0; k < n; k++) {
-               QUAD s;
-               SUFFIX(go_quad_init)(&s, 0);
-               for (i = k; i < m; i++) {
-                       QUAD p;
-                       SUFFIX(go_quad_mul) (&p, &x[i], &V[i][k]);
-                       SUFFIX(go_quad_add) (&s, &s, &p);
-               }
-               SUFFIX(go_quad_add) (&s, &s, &s);
-               for (i = k; i < m; i++) {
-                       QUAD p;
-                       SUFFIX(go_quad_mul) (&p, &s, &V[i][k]);
-                       SUFFIX(go_quad_sub) (&x[i], &x[i], &p);
-               }
-       }
-}
-
+#if 0
 static void
 SUFFIX(print_Q) (CONSTQMATRIX V, int m, int n)
 {
@@ -489,104 +365,17 @@ SUFFIX(print_Q) (CONSTQMATRIX V, int m, int n)
        FREE_MATRIX (Q, m, m);
        g_free (x);
 }
-
-/*
- * Solve Rx=b.
- *
- * Return TRUE in case of error.
- */
-static gboolean
-SUFFIX(back_solve) (CONSTQMATRIX R, QUAD *x, const QUAD *b, int n)
-{
-       int i, j;
-
-       for (i = n - 1; i >= 0; i--) {
-               QUAD d = b[i];
-               for (j = i + 1; j < n; j++) {
-                       QUAD p;
-                       SUFFIX(go_quad_mul) (&p, &R[i][j], &x[j]);
-                       SUFFIX(go_quad_sub) (&d, &d, &p);
-               }
-               if (SUFFIX(go_quad_value)(&R[i][i]) == 0) {
-                       while (i >= 0) SUFFIX(go_quad_init) (&x[i--], 0);
-                       return TRUE;
-               }
-               SUFFIX(go_quad_div) (&x[i], &d, &R[i][i]);
-       }
-
-       return FALSE;
-}
-
-/*
- * Solve RTx=b.
- *
- * Return TRUE in case of error.
- */
-static gboolean
-SUFFIX(fwd_solve) (CONSTQMATRIX R, QUAD *x, const QUAD *b, int n)
-{
-       int i, j;
-
-       for (i = 0; i < n; i++) {
-               QUAD d = b[i];
-               for (j = 0; j < i; j++) {
-                       QUAD p;
-                       SUFFIX(go_quad_mul) (&p, &R[j][i], &x[j]);
-                       SUFFIX(go_quad_sub) (&d, &d, &p);
-               }
-               if (SUFFIX(go_quad_value)(&R[i][i]) == 0) {
-                       while (i < n) SUFFIX(go_quad_init) (&x[i], 0);
-                       return TRUE;
-               }
-               SUFFIX(go_quad_div) (&x[i], &d, &R[i][i]);
-       }
-
-       return FALSE;
-}
-
-static GORegressionResult
-SUFFIX(regres_from_condition) (CONSTQMATRIX R, int n)
-{
-       DOUBLE emin = go_pinf;
-       DOUBLE emax = go_ninf;
-       DOUBLE c, lc;
-       int i;
-
-       for (i = 0; i < n; i++) {
-               DOUBLE e = SUFFIX(fabs)(SUFFIX(go_quad_value)(&R[i][i]));
-               if (e < emin) emin = e;
-               if (e > emax) emax = e;
-       }
-
-       if (emin == 0)
-               return GO_REG_singular;
-
-       /* Condition number */
-       c = emax / emin;
-
-       /* Number of bits destroyed by matrix.  Since the calculations are
-          done with QUADs, we can afford to lose a lot.  */
-       lc = SUFFIX(log)(c) / SUFFIX(log)(FLT_RADIX);
-
-#if 1
-       if (lc > 0.95* DOUBLE_MANT_DIG)
-               return GO_REG_near_singular_bad;
-       if (lc > 0.75 * DOUBLE_MANT_DIG)
-               return GO_REG_near_singular_good;
 #endif
 
-       return GO_REG_ok;
-}
-
 #ifndef BOUNCE
 
 static void
 SUFFIX(calc_residual) (CONSTMATRIX AT, const DOUBLE *b, int m, int n,
-                      const QUAD *y, QUAD *residual, QUAD *N2)
+                      const QUAD *y, QUAD *N2)
 {
        int i, j;
 
-       SUFFIX(go_quad_init) (N2, 0);
+       *N2 = SUFFIX(go_quad_zero);
 
        for (i = 0; i < m; i++) {
                QUAD d;
@@ -600,85 +389,11 @@ SUFFIX(calc_residual) (CONSTMATRIX AT, const DOUBLE *b, int m, int n,
                        SUFFIX(go_quad_sub) (&d, &d, &e);
                }
 
-               if (residual)
-                       residual[i] = d;
-
                SUFFIX(go_quad_mul) (&d, &d, &d);
                SUFFIX(go_quad_add) (N2, N2, &d);
        }
 }
 
-static void
-SUFFIX(refine) (CONSTMATRIX AT, const DOUBLE *b, int m, int n,
-               CONSTQMATRIX V, CONSTQMATRIX R, QUAD *result)
-{
-       QUAD *residual = g_new (QUAD, m);
-       QUAD *newresidual = g_new (QUAD, m);
-       QUAD *QTres = g_new (QUAD, m);
-       QUAD *delta = g_new (QUAD, n);
-       QUAD *newresult = g_new (QUAD, n);
-       int pass;
-       QUAD best_N2;
-
-       SUFFIX(calc_residual) (AT, b, m, n, result, residual, &best_N2);
-#ifdef DEBUG_REFINEMENT
-       g_printerr ("-: Residual norm = %20.15" FORMAT_g "\n",
-                   SUFFIX(go_quad_value) (&best_N2));
-#endif
-
-       for (pass = 1; pass < 10; pass++) {
-               int i;
-               QUAD dres, N2;
-
-               /* Compute Q^T residual.  */
-               for (i = 0; i < m; i++)
-                       QTres[i] = residual[i];
-               SUFFIX(multiply_Qt)(QTres, V, m, n);
-
-               /* Solve R delta = Q^T residual */
-               if (SUFFIX(back_solve) (R, delta, QTres, n))
-                       break;
-
-               /* newresult = result + delta */
-               for (i = 0; i < n; i++) {
-                       SUFFIX(go_quad_add) (&newresult[i],
-                                            &delta[i],
-                                            &result[i]);
-#ifdef DEBUG_REFINEMENT
-                       g_printerr ("d[%2d] = %20.15" FORMAT_f
-                                   "  %20.15" FORMAT_f
-                                   "  %20.15" FORMAT_f "\n",
-                                   i,
-                                   SUFFIX(go_quad_value) (&delta[i]),
-                                   SUFFIX(go_quad_value) (&result[i]),
-                                   SUFFIX(go_quad_value) (&newresult[i]));
-#endif
-               }
-
-               SUFFIX(calc_residual) (AT, b, m, n, newresult, newresidual, &N2);
-               SUFFIX (go_quad_sub) (&dres, &N2, &best_N2);
-
-#ifdef DEBUG_REFINEMENT
-               g_printerr ("%d: Residual norm = %20.15" FORMAT_g
-                           " (%.5" FORMAT_g ")\n",
-                           pass,
-                           SUFFIX(go_quad_value) (&N2),
-                           SUFFIX(go_quad_value) (&dres));
-#endif
-               if (SUFFIX(go_quad_value) (&dres) >= 0)
-                       break;
-
-               best_N2 = N2;
-               COPY_VECTOR (result, newresult, n);
-               COPY_VECTOR (residual, newresidual, m);
-       }
-
-       g_free (residual);
-       g_free (newresidual);
-       g_free (QTres);
-       g_free (newresult);
-       g_free (delta);
-}
 
 #endif
 
@@ -692,11 +407,13 @@ SUFFIX(refine) (CONSTMATRIX AT, const DOUBLE *b, int m, int n,
 GORegressionResult
 SUFFIX(go_linear_solve_multiple) (CONSTMATRIX A, MATRIX B, int n, int bn)
 {
-       QMATRIX V;
-       QMATRIX R;
        void *state;
        GORegressionResult regres;
-       gboolean ok;
+       SUFFIX(GOQuadMatrix) *qA;
+       const SUFFIX(GOQuadMatrix) *R;
+       SUFFIX(GOQuadQR) *qr;
+       QUAD *QTb, *x;
+       int i, j;
 
        if (n < 1 || bn < 1)
                return GO_REG_not_enough_data;
@@ -715,35 +432,37 @@ SUFFIX(go_linear_solve_multiple) (CONSTMATRIX A, MATRIX B, int n, int bn)
 
        state = SUFFIX(go_quad_start) ();
 
-       ALLOC_MATRIX2 (V, n, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
-       ok = SUFFIX(QRH)(A, FALSE, V, R, n, n, NULL);
-       if (ok) {
-               QUAD *QTb = g_new (QUAD, n);
-               QUAD *x = g_new (QUAD, n);
-               int i, j;
+       qA = SUFFIX(quad_matrix_from_matrix) (A, n, n);
+       qr = SUFFIX(go_quad_qr_new) (qA);
+       if (!qr) {
+               regres = GO_REG_invalid_data;
+               goto done;
+       }
+       R = SUFFIX(go_quad_qr_r) (qr);
 
-               regres = SUFFIX(regres_from_condition)(R, n);
+       QTb = g_new (QUAD, n);
+       x = g_new (QUAD, n);
 
-               for (j = 0; j < bn; j++) {
-                       /* Compute Q^T b.  */
-                       for (i = 0; i < n; i++)
-                               SUFFIX(go_quad_init)(&QTb[i], B[i][j]);
-                       SUFFIX(multiply_Qt)(QTb, V, n, n);
+       for (j = 0; j < bn; j++) {
+               /* Compute Q^T b.  */
+               for (i = 0; i < n; i++)
+                       SUFFIX(go_quad_init)(&QTb[i], B[i][j]);
+               SUFFIX(go_quad_qr_multiply_qt)(qr, QTb);
 
-                       /* Solve R x = Q^T b  */
-                       if (SUFFIX(back_solve) (R, x, QTb, n))
-                               regres = GO_REG_singular;
+               /* Solve R x = Q^T b  */
+               if (SUFFIX(go_quad_matrix_back_solve) (R, x, QTb))
+                       regres = GO_REG_singular;
 
-                       /* Round for output.  */
-                       for (i = 0; i < n; i++)
-                               B[i][j] = SUFFIX(go_quad_value) (&x[i]);
-               }
+               /* Round for output.  */
+               for (i = 0; i < n; i++)
+                       B[i][j] = SUFFIX(go_quad_value) (&x[i]);
+       }
 
-               g_free (x);
-               g_free (QTb);
-       } else
-               regres = GO_REG_invalid_data;
+       SUFFIX(go_quad_qr_free) (qr);
+       g_free (x);
+       g_free (QTb);
+done:
+       SUFFIX(go_quad_matrix_free) (qA);
 
        SUFFIX(go_quad_end) (state);
 
@@ -778,99 +497,42 @@ SUFFIX(go_linear_solve) (CONSTMATRIX A, const DOUBLE *b, int n, DOUBLE *res)
 gboolean
 SUFFIX(go_matrix_invert) (MATRIX A, int n)
 {
-       QMATRIX V;
-       QMATRIX R;
-       gboolean has_result;
-       void *state = SUFFIX(go_quad_start) ();
-
-       ALLOC_MATRIX2 (V, n, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
-
-       has_result = SUFFIX(QRH)(A, FALSE, V, R, n, n, NULL);
-
-       if (has_result) {
-               int i, k;
-               QUAD *x = g_new (QUAD, n);
-               QUAD *QTk = g_new (QUAD, n);
-               GORegressionResult regres =
-                       SUFFIX(regres_from_condition) (R, n);
-
-               if (regres != GO_REG_ok && regres != GO_REG_near_singular_good)
-                       has_result = FALSE;
-
-               for (k = 0; has_result && k < n; k++) {
-                       /* Compute Q^T's k-th column.  */
-                       for (i = 0; i < n; i++)
-                               SUFFIX(go_quad_init)(&QTk[i], i == k);
-                       SUFFIX(multiply_Qt)(QTk, V, n, n);
-
-                       /* Solve R x = Q^T e_k */
-                       if (SUFFIX(back_solve) (R, x, QTk, n))
-                               has_result = FALSE;
-
-                       /* Round for output.  */
-                       for (i = 0; i < n; i++)
-                               A[i][k] = SUFFIX(go_quad_value) (&x[i]);
-               }
+       SUFFIX(GOQuadMatrix) *qA, *qZ;
+       DOUBLE threshold = DEFAULT_THRESHOLD;
+       void *state;
+       gboolean ok;
 
-               g_free (QTk);
-               g_free (x);
+       state = SUFFIX(go_quad_start) ();
+       qA = SUFFIX(quad_matrix_from_matrix) (A, n, n);
+       qZ = SUFFIX(go_quad_matrix_inverse) (qA, threshold);
+       ok = (qZ != NULL);
+       SUFFIX(go_quad_matrix_free) (qA);
+       if (qZ) {
+               SUFFIX(copy_quad_matrix_to_matrix) (A, qZ);
+               SUFFIX(go_quad_matrix_free) (qZ);
        }
-
-       FREE_MATRIX (V, n, n);
-       FREE_MATRIX (R, n, n);
-
        SUFFIX(go_quad_end) (state);
-
-       return has_result;
+       return ok;
 }
 
 DOUBLE
 SUFFIX(go_matrix_determinant) (CONSTMATRIX A, int n)
 {
-       gboolean ok;
-       QMATRIX R;
-       QMATRIX V;
-       int qdet;
        DOUBLE res;
+       QUAD qres;
        void *state;
+       SUFFIX(GOQuadMatrix) *qA;
 
        if (n < 1)
                return 0;
 
-       /* Special case.  */
-       if (n == 1)
-               return A[0][0];
-
        state = SUFFIX(go_quad_start) ();
 
-       /* Special case.  */
-       if (n == 2) {
-               QUAD a, b, det;
-               SUFFIX(go_quad_mul12)(&a, A[0][0], A[1][1]);
-               SUFFIX(go_quad_mul12)(&b, A[1][0], A[0][1]);
-               SUFFIX(go_quad_sub)(&det, &a, &b);
-               res = SUFFIX(go_quad_value)(&det);
-               goto done;
-       }
+       qA = SUFFIX(quad_matrix_from_matrix) (A, n, n);
+       SUFFIX(go_quad_matrix_determinant) (qA, &qres);
+       SUFFIX(go_quad_matrix_free) (qA);
+       res = SUFFIX(go_quad_value) (&qres);
 
-       ALLOC_MATRIX2 (V, n, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
-       ok = SUFFIX(QRH)(A, FALSE, V, R, n, n, &qdet);
-       if (ok) {
-               int i;
-               QUAD det;
-               SUFFIX(go_quad_init)(&det, qdet);
-               for (i = 0; i < n; i++)
-                       SUFFIX(go_quad_mul)(&det, &det, &R[i][i]);
-               res = SUFFIX(go_quad_value)(&det);
-       } else
-               res = go_nan;
-
-       FREE_MATRIX (V, n, n);
-       FREE_MATRIX (R, n, n);
-
-done:
        SUFFIX(go_quad_end) (state);
 
        return res;
@@ -910,12 +572,12 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
        g_free (ys2);
        FREE_MATRIX (xssT2, n, m);
 #else
-       QMATRIX V;
-       QMATRIX R;
+       SUFFIX(GOQuadMatrix) *xss;
+       SUFFIX(GOQuadQR) *qr;
        QUAD *qresult;
        QUAD *QTy;
        QUAD *inv;
-       int i, k;
+       int i, j, k;
        gboolean has_result;
        void *state;
        gboolean has_stat;
@@ -935,29 +597,32 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 
        state = SUFFIX(go_quad_start) ();
 
-       ALLOC_MATRIX2 (V, m, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
+       xss = SUFFIX(go_quad_matrix_new) (m, n);
+       for (i = 0; i < m; i++)
+               for (j = 0; j < n; j++)
+                       SUFFIX(go_quad_init) (&xss->data[i][j], xssT[j][i]);
+       qr = SUFFIX(go_quad_qr_new) (xss);
+       SUFFIX(go_quad_matrix_free) (xss);
+
        qresult = g_new0 (QUAD, n);
        QTy = g_new (QUAD, m);
        inv = g_new (QUAD, n);
 
-       has_result = SUFFIX(QRH) (xssT, TRUE, V, R, m, n, NULL);
+       has_result = (qr != NULL);
 
        if (has_result) {
+               const SUFFIX(GOQuadMatrix) *R = SUFFIX(go_quad_qr_r) (qr);
                regerr = GO_REG_ok;
 
                /* Compute Q^T ys.  */
                for (i = 0; i < m; i++)
                        SUFFIX(go_quad_init)(&QTy[i], ys[i]);
-               SUFFIX(multiply_Qt)(QTy, V, m, n);
+               SUFFIX(go_quad_qr_multiply_qt)(qr, QTy);
 
                /* Solve R res = Q^T ys */
-               if (SUFFIX(back_solve) (R, qresult, QTy, n))
+               if (SUFFIX(go_quad_matrix_back_solve) (R, qresult, QTy))
                        has_result = FALSE;
 
-               if (n > 2)
-                       SUFFIX(refine) (xssT, ys, m, n, V, R, qresult);
-
                /* Round to plain precision.  */
                for (i = 0; i < n; i++) {
                        result[i] = SUFFIX(go_quad_value) (&qresult[i]);
@@ -982,7 +647,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
                        g_assert (err == 0);
                }
 
-               SUFFIX(calc_residual) (xssT, ys, m, n, qresult, NULL, &N2);
+               SUFFIX(calc_residual) (xssT, ys, m, n, qresult, &N2);
                stat_->ss_resid = SUFFIX(go_quad_value) (&N2);
 
                stat_->sqr_r = (stat_->ss_total == 0)
@@ -1000,7 +665,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
                stat_->adj_sqr_r = 1 - stat_->ss_resid * (m - 1) /
                        ((m - n) * stat_->ss_total);
                if (m == n)
-                       SUFFIX(go_quad_init) (&N2, 0);
+                       N2 = SUFFIX(go_quad_zero);
                else {
                        QUAD d;
                        SUFFIX(go_quad_init) (&d, m - n);
@@ -1017,13 +682,13 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
                                SUFFIX(go_quad_init) (&inv[i], i == k ? 1 : 0);
 
                        /* Solve R^T inv = e_k */
-                       if (SUFFIX(fwd_solve) (R, inv, inv, n)) {
+                       if (SUFFIX(go_quad_matrix_fwd_solve) (R, inv, inv)) {
                                regerr = GO_REG_singular;
                                break;
                        }
 
                        /* Solve R newinv = inv */
-                       if (SUFFIX(back_solve) (R, inv, inv, n)) {
+                       if (SUFFIX(go_quad_matrix_back_solve) (R, inv, inv)) {
                                regerr = GO_REG_singular;
                                break;
                        }
@@ -1077,8 +742,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
        g_free (QTy);
        g_free (qresult);
        g_free (inv);
-       FREE_MATRIX (V, m, n);
-       FREE_MATRIX (R, n, n);
+       if (qr) SUFFIX(go_quad_qr_free) (qr);
 out:
        if (!has_stat)
                SUFFIX(go_regression_stat_destroy) (stat_);
@@ -2113,49 +1777,49 @@ SUFFIX(go_non_linear_regression) (SUFFIX(GORegressionFunction) f,
 GORegressionResult
 SUFFIX(go_linear_regression_leverage) (MATRIX A, DOUBLE *d, int m, int n)
 {
-       QMATRIX V;
-       QMATRIX R;
-       gboolean has_result;
+       SUFFIX(GOQuadMatrix) *qA;
+       SUFFIX(GOQuadQR) *qr;
        void *state = SUFFIX(go_quad_start) ();
        GORegressionResult regres;
+       DOUBLE threshold = DEFAULT_THRESHOLD;
 
-       ALLOC_MATRIX2 (V, m, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
-
-       has_result = SUFFIX(QRH)(A, FALSE, V, R, m, n, NULL);
-
-       if (has_result) {
+       qA = SUFFIX(quad_matrix_from_matrix) (A, m, n);
+       qr = SUFFIX(go_quad_qr_new) (qA);
+       if (qr) {
                int k;
                QUAD *b = g_new (QUAD, n);
+               const SUFFIX(GOQuadMatrix) *R = SUFFIX(go_quad_qr_r) (qr);
+               DOUBLE emin, emax;
 
-               regres = SUFFIX(regres_from_condition) (R, n);
+               SUFFIX(go_quad_matrix_eigen_range) (R, &emin, &emax);
+               regres = (emin > emax * threshold)
+                       ? GO_REG_ok
+                       : GO_REG_singular;
 
                for (k = 0; k < m; k++) {
                        int i;
-                       QUAD acc;
+                       QUAD acc = SUFFIX(go_quad_zero);
 
                        /* b = AT e_k  */
                        for (i = 0; i < n; i++)
-                               SUFFIX(go_quad_init) (&b[i], A[k][i]);
+                               b[i] = qA->data[k][i];
 
-                       /* Solve R^T newb = b */
-                       if (SUFFIX(fwd_solve) (R, b, b, n)) {
+                       /* Solve R^T b = AT e_k */
+                       if (SUFFIX(go_quad_matrix_fwd_solve) (R, b, b)) {
                                regres = GO_REG_singular;
                                break;
                        }
 
                        /* Solve R newb = b */
-                       if (SUFFIX(back_solve) (R, b, b, n)) {
+                       if (SUFFIX(go_quad_matrix_back_solve) (R, b, b)) {
                                regres = GO_REG_singular;
                                break;
                        }
 
                        /* acc = (Ab)_k */
-                       SUFFIX(go_quad_init) (&acc, 0);
                        for (i = 0; i < n; i++) {
                                QUAD p;
-                               SUFFIX(go_quad_init) (&p, A[k][i]);
-                               SUFFIX(go_quad_mul) (&p, &p, &b[i]);
+                               SUFFIX(go_quad_mul) (&p, &qA->data[k][i], &b[i]);
                                SUFFIX(go_quad_add) (&acc, &acc, &p);
                        }
 
@@ -2163,41 +1827,15 @@ SUFFIX(go_linear_regression_leverage) (MATRIX A, DOUBLE *d, int m, int n)
                }
 
                g_free (b);
+               SUFFIX(go_quad_qr_free) (qr);
        } else
                regres = GO_REG_invalid_data;
 
-       FREE_MATRIX (V, m, n);
-       FREE_MATRIX (R, n, n);
-
        SUFFIX(go_quad_end) (state);
 
        return regres;
 }
 
-static void
-SUFFIX(qmatrix_mul) (QMATRIX C, QMATRIX A, QMATRIX B,
-                    int C_m, int C_n, int A_n)
-{
-       int i, j, k;
-
-       for (i = 0; i < C_n; i++) {
-               for (j = 0; j < C_m; j++) {
-                       QUAD acc;
-
-                       SUFFIX(go_quad_init) (&acc, 0);
-                       for (k = 0; k < A_n; k++) {
-                               QUAD p;
-
-                               SUFFIX(go_quad_mul) (&p, &A[i][k], &B[k][j]);
-                               SUFFIX(go_quad_add) (&acc, &acc, &p);
-                       }
-
-                       C[i][j] = acc;
-               }
-       }
-}
-
-
 /*
  * Compute the pseudo-inverse of a matrix, see
  *     http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse
@@ -2220,216 +1858,16 @@ SUFFIX(go_matrix_pseudo_inverse) (CONSTMATRIX A, int m, int n,
                                  DOUBLE threshold,
                                  MATRIX B)
 {
-       QMATRIX V;
-       QMATRIX R;
-       QMATRIX Bi = NULL;
-       MATRIX RTR = NULL;
-       gboolean has_result;    
+       SUFFIX(GOQuadMatrix) *qA, *qZ;
        void *state;
-       int i, j;
-       DOUBLE emin, emax, delta;
-       QUAD *x;
-       int steps;
-       gboolean full_rank;
-       gboolean debug = FALSE;
-
-       if (debug) {
-               g_printerr ("A:\n");
-               PRINT_MATRIX (A, m, n);
-       }
-
-       if (m < n) {
-               MATRIX AT;
-               MATRIX BT;
-
-               if (debug)
-                       g_printerr ("Pseudo-inverse via transpose\n");
-
-               /*
-                * The code below assumes m >= n.  Luckily, taking the
-                * pseudo-inverse commutes with transposition.
-                */
-
-               ALLOC_MATRIX (AT, n, m);
-               ALLOC_MATRIX (BT, m, n);
-               for (i = 0; i < m; i++)
-                       for (j = 0; j < n; j++)
-                               AT[j][i] = A[i][j];
-               SUFFIX(go_matrix_pseudo_inverse)(AT, n, m, threshold, BT);
-               for (i = 0; i < m; i++)
-                       for (j = 0; j < n; j++)
-                               B[j][i] = BT[i][j];
-               FREE_MATRIX (AT, n, m);
-               FREE_MATRIX (BT, m, n);
-               return;
-       }
 
        state = SUFFIX(go_quad_start) ();
-
-       ALLOC_MATRIX2 (V, m, n, QUAD);
-       ALLOC_MATRIX2 (R, n, n, QUAD);
-
-       has_result = SUFFIX(QRH)(A, FALSE, V, R, m, n, NULL);
-       if (!has_result) {
-               if (debug)
-                       g_printerr ("QRH failed\n");
-               goto null_result;
-       }
-
-       if (debug) {
-               g_printerr ("Q:\n");
-               SUFFIX(print_Q) (V, m, n);
-               g_printerr ("R:\n");
-               PRINT_QMATRIX (R, n, n);
+       qA = SUFFIX(quad_matrix_from_matrix) (A, m, n);
+       qZ = SUFFIX(go_quad_matrix_pseudo_inverse) (qA, threshold);
+       SUFFIX(go_quad_matrix_free) (qA);
+       if (qZ) {
+               SUFFIX(copy_quad_matrix_to_matrix) (B, qZ);
+               SUFFIX(go_quad_matrix_free) (qZ);
        }
-
-       emax = 0;
-       for (i = 0; i < n; i++) {
-               DOUBLE abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value)(&R[i][i]));
-               emax = MAX (emax, abs_e);
-       }
-       if (emax == 0)
-               goto null_result;
-
-       emin = emax;
-       full_rank = TRUE;
-       for (i = 0; i < n; i++) {
-               DOUBLE abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value)(&R[i][i]));
-               if (abs_e <= emax * threshold) {
-                       full_rank = FALSE;
-                       SUFFIX(go_quad_init)(&R[i][i], 0);
-               } else
-                       emin = MIN (emin, abs_e);
-       }
-
-       delta = full_rank ? 0 : emax * threshold;
-
-       /*
-        * Starting point for the iteration:
-        *
-        *    Bi := (RT R + delta I)^-1 * RT
-        *
-        * (RT R) is positive semi-definite and (delta I) is positive definite,
-        * so the sum is positive definite and invertible.
-        */
-       ALLOC_MATRIX (RTR, n, n);
-       ALLOC_MATRIX2 (Bi, n, n, QUAD);
-       for (i = 0; i < n; i++) {
-               for (j = 0; j < n; j++) {
-                       QUAD acc;
-                       int k;
-
-                       SUFFIX(go_quad_init) (&acc, i == j ? delta : 0.0);
-                       for (k = 0; k <= i && k <= j; k++) {
-                               QUAD p;
-
-                               SUFFIX(go_quad_mul) (&p, &R[k][i], &R[k][j]);
-                               SUFFIX(go_quad_add) (&acc, &acc, &p);
-                       }
-
-                       RTR[i][j] = SUFFIX(go_quad_value) (&acc);
-               }
-       }
-       if (debug) {
-               g_printerr ("RTR:\n");
-               PRINT_MATRIX (RTR, n, n);
-       }
-       if (!SUFFIX(go_matrix_invert)(RTR, n))
-               goto null_result;
-       if (debug) {
-               g_printerr ("RTR inverse:\n");
-               PRINT_MATRIX (RTR, n, n);
-       }
-       for (i = 0; i < n; i++) {
-               for (j = 0; j < n; j++) {
-                       QUAD acc;
-                       int k;
-
-                       SUFFIX(go_quad_init) (&acc, 0.0);
-                       for (k = 0; k < n; k++) {
-                               QUAD p, a;
-
-                               SUFFIX(go_quad_init) (&a, RTR[i][k]);
-                               SUFFIX(go_quad_mul) (&p, &a, &R[j][k]);
-                               SUFFIX(go_quad_add) (&acc, &acc, &p);
-                       }
-
-                       Bi[i][j] = acc;
-               }
-       }
-       if (debug) {
-               g_printerr ("B0:\n");
-               PRINT_QMATRIX (Bi, n, n);
-       }
-
-       /* B_{i+1} = (2 I - B_i R) B_i */
-       for (steps = 0; steps < 10; steps++) {
-               QMATRIX W;
-               QMATRIX Bip1;
-               QUAD two;
-
-               SUFFIX(go_quad_init)(&two, 2);
-
-               ALLOC_MATRIX2 (W, n, n, QUAD);
-               ALLOC_MATRIX2 (Bip1, n, n, QUAD);
-
-               SUFFIX(qmatrix_mul) (W, Bi, R, n, n, n);
-               for (i = 0; i < n; i++)
-                       for (j = 0; j < n; j++)
-                               SUFFIX(go_quad_sub) (&W[i][j],
-                                                    i == j ? &two : &SUFFIX(go_quad_zero),
-                                                    &W[i][j]);
-               SUFFIX(qmatrix_mul) (Bip1, W, Bi, n, n, n);
-               COPY_MATRIX (Bi, Bip1, n, n);
-               if (debug) {
-                       g_printerr ("B%d:\n", steps + 1);
-                       PRINT_QMATRIX (Bi, n, n);
-               }
-
-               FREE_MATRIX (W, n, n);
-               FREE_MATRIX (Bip1, n, n);
-       }
-
-       /* B := (Bi|O) Q^T */
-       x = g_new (QUAD, m);
-       for (j = 0; j < m; j++) {
-               QUAD acc;
-
-               /* Compute Q^T e_j.  */
-               for (i = 0; i < m; i++)
-                       SUFFIX(go_quad_init)(&x[i], i == j ? 1 : 0);
-               SUFFIX(multiply_Qt) (x, V, m, n);
-
-               SUFFIX(go_quad_init) (&acc, 0);
-               for (i = 0; i < n; i++) {
-                       QUAD acc;
-                       int k;
-
-                       SUFFIX(go_quad_init) (&acc, 0);
-                       for (k = 0; k < n /* Only n */; k++) {
-                               QUAD p;
-                               SUFFIX(go_quad_mul) (&p, &Bi[i][k], &x[k]);
-                               SUFFIX(go_quad_add) (&acc, &acc, &p);
-                       }
-
-                       B[i][j] = SUFFIX(go_quad_value)(&acc);
-               }
-       }
-       g_free (x);
-
-out:
-       FREE_MATRIX (V, m, n);
-       FREE_MATRIX (R, n, n);
-       if (RTR) FREE_MATRIX (RTR, n, n);
-       if (Bi) FREE_MATRIX (Bi, n, n);
-
        SUFFIX(go_quad_end) (state);
-       return;
-
-null_result:
-       if (debug) g_printerr ("Null result\n");
-       for (i = 0; i < n; i++)
-               for (j = 0; j < m; j++)
-                       B[i][j] = 0;
-       goto out;
 }
diff --git a/goffice/math/goffice-math.h b/goffice/math/goffice-math.h
index 66be62f..b9fd172 100644
--- a/goffice/math/goffice-math.h
+++ b/goffice/math/goffice-math.h
@@ -4,10 +4,14 @@
 /* Forward stuff */
 typedef struct GOAccumulator_ GOAccumulator;
 typedef struct GOQuad_ GOQuad;
+typedef struct GOQuadMatrix_ GOQuadMatrix;
+typedef struct GOQuadQR_ GOQuadQR;
 
 #ifdef GOFFICE_WITH_LONG_DOUBLE
 typedef struct GOAccumulatorl_ GOAccumulatorl;
 typedef struct GOQuadl_ GOQuadl;
+typedef struct GOQuadMatrixl_ GOQuadMatrixl;
+typedef struct GOQuadQRl_ GOQuadQRl;
 #endif
 
 #include <goffice/math/go-accumulator.h>
@@ -16,6 +20,7 @@ typedef struct GOQuadl_ GOQuadl;
 #include <goffice/math/go-distribution.h>
 #include <goffice/math/go-fft.h>
 #include <goffice/math/go-math.h>
+#include <goffice/math/go-matrix.h>
 #include <goffice/math/go-matrix3x3.h>
 #include <goffice/math/go-quad.h>
 #include <goffice/math/go-R.h>



[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]