This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
unsymmetric eigenvalue problem
- From: Patrick dot Alken at colorado dot edu
- To: gsl-discuss at sourceware dot org
- Date: Wed, 17 May 2006 14:45:10 -0600
- Subject: unsymmetric eigenvalue problem
Hello all,
I have recently needed to compute eigenvalues of unsymmetric
matrices, and so I wrote a code based in GSL to do it since
one does not currently exist. I followed the QR algorithm
outlined in Golub & Van Loan, chapter 7.
I have attached a patch file to this email containing all
my code. Right now only the eigenvalue problem is solved,
the code will not compute eigenvectors. The patch is against
gsl-1.8 so just cd into a fresh gsl-1.8 and
patch -p0 < gsl.unsymm.patch
You will probably need to run automake again in the toplevel
dir to update the Makefile.am's. I have also attached a
test program eig.c to illustrate the use of the solver.
The patch creates 3 new source files:
linalg/hessenberg.c
eigen/balance.c
eigen/unsymm.c
and 1 header file
eigen/balance.h
hessenberg.c contains gsl_linalg_hessenberg() which reduces
a matrix to upper Hessenberg form.
balance.c contains balance_matrix() to do the standard matrix
balancing before finding eigenvalues.
unsymm.c contains the Francis QR algorithm to find the eigenvalues
of a general unsymmetric real matrix:
gsl_eigen_unsymm_alloc()
gsl_eigen_unsymm_free()
gsl_eigen_unsymm() --- Solves A x = \lambda x
I am happy to contribute this to GSL if the developers want to
use it. Also, if there is interest, I would be willing to look
at coding up an eigenvector solver too for the unsymmetric
eigenvalue program - I only needed eigenvalues for my current
project so I haven't looked at the eigenvector problem.
Also if there is interest I'd be interested in looking at the
generalized eigenvalue problem since I feel GSL is incomplete
without an Ax = sBx solver.
Please test and give feedback.
Thanks,
Patrick Alken
/*
* tester program for GSL unsymmetric eigenvalue solver
*
* To compile: gcc -o eig eig.c -lgsl -lgslcblas -lm
*/
#include <stdio.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_vector.h>
#define N 4
int
main(void)
{
#if N == 4
double a_data[] = { 1.0, 2.0, 3.0, 4.0,
3.4, 2.2, 1.0, 9.7,
6.7, 5.0, 33.0, 8.9,
78.7, 0.0, 4.5, 3.3 };
#endif
#if N == 2
double a_data[] = { 3.0, -2.0,
4.0, -1.0 };
#endif
#if N == 3
double a_data[] = { 3.0, -2.0, 0.0,
4.0, -1.0, 0.0,
1.0, 0.0, 1.0 };
#endif
gsl_matrix_view m;
gsl_eigen_unsymm_workspace *w;
gsl_vector_complex *ev;
int i,j;
m = gsl_matrix_view_array(a_data, N, N);
printf("original matrix\n");
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
printf("%e ", gsl_matrix_get(&m.matrix, i, j));
printf("\n");
}
w = gsl_eigen_unsymm_alloc(N);
ev = gsl_vector_complex_alloc(N);
gsl_eigen_unsymm(&m.matrix, ev, w);
printf("Schur decomposed matrix\n");
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
printf("%e ", gsl_matrix_get(&m.matrix, i, j));
printf("\n");
}
for (i = 0; i < N; ++i)
{
gsl_complex a = gsl_vector_complex_get(ev, i);
if (GSL_IMAG(a) == 0.0)
{
printf("ev[%d] = %f\n",
i,
GSL_REAL(a));
}
else
{
printf("ev[%d] = %f + %fi\n",
i,
GSL_REAL(a),
GSL_IMAG(a));
}
}
gsl_eigen_unsymm_free(w);
gsl_vector_complex_free(ev);
return 0;
}
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
--- /home/palken/tmp2/gsl-1.8/eigen/balance.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.c 2006-05-17 11:57:52.000000000 -0600
@@ -0,0 +1,117 @@
+/* eigen/balance.c
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * This program 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) any later version.
+ *
+ * This program 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <stdlib.h>
+#include <math.h>
+
+#include "balance.h"
+
+/*
+balance_matrix()
+ Balance a given matrix by applying a diagonal matrix
+similarity transformation so that the rows and columns
+of the new matrix have norms which are the same order of
+magnitude. This is necessary for the unsymmetric eigenvalue
+problem since the calculation can become numerically unstable
+for unbalanced matrices.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
+
+and
+
+Numerical Recipes in C, Section 11.5
+*/
+
+void
+balance_matrix(gsl_matrix * A)
+{
+ double row_norm,
+ col_norm;
+ int not_converged;
+ const size_t N = A->size1;
+
+ not_converged = 1;
+
+ while (not_converged)
+ {
+ size_t i, j;
+ double g, f, s;
+
+ not_converged = 0;
+
+ for (i = 0; i < N; ++i)
+ {
+ row_norm = 0.0;
+ col_norm = 0.0;
+
+ for (j = 0; j < N; ++j)
+ {
+ if (j != i)
+ {
+ col_norm += fabs(gsl_matrix_get(A, j, i));
+ row_norm += fabs(gsl_matrix_get(A, i, j));
+ }
+ }
+
+ if ((col_norm == 0.0) || (row_norm == 0.0))
+ {
+ continue;
+ }
+
+ g = row_norm / FLOAT_RADIX;
+ f = 1.0;
+ s = col_norm + row_norm;
+
+ /*
+ * find the integer power of the machine radix which
+ * comes closest to balancing the matrix
+ */
+ while (col_norm < g)
+ {
+ f *= FLOAT_RADIX;
+ col_norm *= FLOAT_RADIX_SQ;
+ }
+
+ g = row_norm * FLOAT_RADIX;
+
+ while (col_norm > g)
+ {
+ f /= FLOAT_RADIX;
+ col_norm /= FLOAT_RADIX_SQ;
+ }
+
+ if ((row_norm + col_norm) < 0.95 * s * f)
+ {
+ not_converged = 1;
+
+ g = 1.0 / f;
+
+ /* apply similarity transformation */
+ for (j = 0; j < N; ++j)
+ {
+ gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) * g);
+ }
+ for (j = 0; j < N; ++j)
+ {
+ gsl_matrix_set(A, j, i, gsl_matrix_get(A, j, i) * f);
+ }
+ }
+ }
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
--- /home/palken/tmp2/gsl-1.8/eigen/balance.h 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.h 2006-05-17 10:49:00.000000000 -0600
@@ -0,0 +1,32 @@
+/* eigen/balance.h
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * This program 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) any later version.
+ *
+ * This program 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __GSL_BALANCE_H__
+#define __GSL_BALANCE_H__
+
+#include <gsl/gsl_matrix.h>
+
+/* floating point radix */
+#define FLOAT_RADIX 2.0
+
+#define FLOAT_RADIX_SQ (FLOAT_RADIX * FLOAT_RADIX)
+
+void balance_matrix(gsl_matrix * A);
+
+#endif /* __GSL_BALANCE_H__ */
diff -urN /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
--- /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h 2005-06-26 07:25:34.000000000 -0600
+++ ./eigen/gsl_eigen.h 2006-05-17 12:02:26.000000000 -0600
@@ -59,6 +59,16 @@
typedef struct {
size_t size;
+ gsl_vector *v; /* temporary vector */
+} gsl_eigen_unsymm_workspace;
+
+gsl_eigen_unsymm_workspace * gsl_eigen_unsymm_alloc (const size_t n);
+void gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w);
+int gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_eigen_unsymm_workspace * w);
+
+typedef struct {
+ size_t size;
double * d;
double * sd;
double * tau;
diff -urN /home/palken/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
--- /home/palken/tmp2/gsl-1.8/eigen/Makefile.am 2004-09-11 07:45:45.000000000 -0600
+++ ./eigen/Makefile.am 2006-05-17 12:00:33.000000000 -0600
@@ -3,7 +3,7 @@
check_PROGRAMS = test
pkginclude_HEADERS = gsl_eigen.h
-libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c herm.c hermv.c sort.c
+libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c balance.c
INCLUDES= -I$(top_builddir)
noinst_HEADERS = qrstep.c
diff -urN /home/palken/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
--- /home/palken/tmp2/gsl-1.8/eigen/unsymm.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/unsymm.c 2006-05-17 11:57:57.000000000 -0600
@@ -0,0 +1,378 @@
+/* eigen/unsymm.c
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * This program 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) any later version.
+ *
+ * This program 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <stdlib.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_vector_complex.h>
+#include <gsl/gsl_matrix.h>
+
+#include "balance.h"
+
+/* maximum iterations before failure is reported */
+#define UNSYMM_MAX_ITERATIONS 100
+
+/*
+ * This module computes the eigenvalues of a real unsymmetric
+ * matrix, using the QR decomposition.
+ *
+ * See Golub & Van Loan, "Matrix Computations" (3rd ed),
+ * algorithm 7.5.2
+ */
+
+inline static void zero_subdiag_small_elements(gsl_matrix * A);
+static inline int francis_qrstep(gsl_matrix * H,
+ gsl_eigen_unsymm_workspace * w);
+static void get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+ gsl_complex * e2);
+
+/*
+gsl_eigen_unsymm_alloc()
+
+Allocate a workspace for solving the unsymmetric
+eigenvalue/eigenvector problem
+
+Inputs: n - size of matrix
+
+Return: pointer to workspace
+*/
+
+gsl_eigen_unsymm_workspace *
+gsl_eigen_unsymm_alloc(const size_t n)
+{
+ gsl_eigen_unsymm_workspace *w;
+
+ if (n == 0)
+ {
+ GSL_ERROR_NULL ("matrix dimension must be positive integer",
+ GSL_EINVAL);
+ }
+
+ w = ((gsl_eigen_unsymm_workspace *)
+ malloc (sizeof (gsl_eigen_unsymm_workspace)));
+
+ if (w == 0)
+ {
+ GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
+ }
+
+ w->size = n;
+ w->v = gsl_vector_alloc(3);
+
+ return (w);
+}
+
+void
+gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
+{
+ gsl_vector_free(w->v);
+
+ free(w);
+}
+
+/*
+gsl_eigen_unsymm()
+
+Solve the unsymmetric eigenvalue problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda using algorithm 7.5.2 of
+Golub & Van Loan, "Matrix Computations" (3rd ed)
+
+Inputs: A - matrix
+ eval - where to store eigenvalues
+ w - workspace
+
+Notes: the matrix A is destroyed by this routine
+*/
+
+int
+gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+ gsl_eigen_unsymm_workspace * w)
+{
+ /* check matrix and vector sizes */
+
+ if (A->size1 != A->size2)
+ {
+ GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
+ }
+ else if (eval->size != A->size1)
+ {
+ GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
+ }
+ else
+ {
+ size_t N;
+ size_t q;
+ gsl_complex lambda1, lambda2; /* eigenvalues */
+ gsl_matrix_view m;
+ size_t nev;
+ size_t nit;
+
+ N = A->size1;
+
+ /* special cases */
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(A, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, 0, lambda1);
+ return GSL_SUCCESS;
+ }
+
+ if (N == 2)
+ {
+ /*
+ * The 2x2 case is special since the matrix is already
+ * in upper quasi-triangular form so no Schur decomposition
+ * is necessary
+ */
+ get_2b2_eigenvalues(A, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, 0, lambda1);
+ gsl_vector_complex_set(eval, 1, lambda2);
+ return GSL_SUCCESS;
+ }
+
+ nev = 0;
+
+ /* balance the matrix */
+ balance_matrix(A);
+
+ /* compute the Hessenberg reduction of A */
+ gsl_linalg_hessenberg(A);
+
+ /* set small elements on the subdiagonal to 0 */
+ zero_subdiag_small_elements(A);
+
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+
+ nit = 0;
+ while ((N > 2) && (nit++ < UNSYMM_MAX_ITERATIONS))
+ {
+ francis_qrstep(&m.matrix, w);
+ zero_subdiag_small_elements(&m.matrix);
+
+ q = N - 1;
+ while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
+ {
+ --q;
+ }
+
+ if (q == (N - 1))
+ {
+ /*
+ * the (N, N - 1) element of the matrix is 0 -
+ * m_{NN} is a real eigenvalue
+ */
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, q, q), 0.0);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+
+ --N;
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+ }
+ else if (q == (N - 2))
+ {
+ /*
+ * The bottom right 2x2 block of m is an eigenvalue
+ * system
+ */
+
+ m = gsl_matrix_submatrix(A, q, q, 2, 2);
+ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+
+ N -= 2;
+ m = gsl_matrix_submatrix(A, 0, 0, N, N);
+ }
+ }
+
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ }
+ else if (N == 2)
+ {
+ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+ }
+
+ return GSL_SUCCESS;
+ }
+} /* gsl_eigen_unsymm() */
+
+/********************************************
+ * INTERNAL ROUTINES *
+ ********************************************/
+
+/*
+zero_subdiag_small_elements()
+ Sets to zero all elements on the subdiaganal of a matrix A
+which satisfy
+
+|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
+*/
+
+inline static void
+zero_subdiag_small_elements(gsl_matrix * A)
+{
+ const size_t N = A->size1;
+ size_t i;
+ double dpel = gsl_matrix_get(A, 0, 0);
+
+ for (i = 1; i < N; ++i)
+ {
+ double sel = gsl_matrix_get(A, i, i - 1);
+ double del = gsl_matrix_get(A, i, i);
+
+ if (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel)))
+ {
+ gsl_matrix_set(A, i, i - 1, 0.0);
+ }
+
+ dpel = del;
+ }
+}
+
+/*
+francis_qrstep()
+ Perform a Francis QR step.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed),
+algorithm 7.5.1
+
+Inputs: H - unreduced upper Hessenberg matrix
+ w - workspace
+*/
+
+static inline int
+francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w)
+{
+ const size_t N = H->size1;
+ double s, t, x, y, z;
+ size_t i;
+ gsl_vector_view v;
+ gsl_matrix_view m;
+ double tau_i;
+ size_t q, r;
+
+ /* s = a_1 + a_2 = H_{mm} + H_{nn} */
+ s = gsl_matrix_get(H, N - 2, N - 2) + gsl_matrix_get(H, N - 1, N - 1);
+
+ /* t = a_1 * a_2 = H_{mm} * H_{nn} - H_{mn} * H_{nm} */
+ t = gsl_matrix_get(H, N - 2, N - 2) * gsl_matrix_get(H, N - 1, N - 1) -
+ gsl_matrix_get(H, N - 2, N - 1) * gsl_matrix_get(H, N - 1, N - 2);
+ x = gsl_matrix_get(H, 0, 0) * gsl_matrix_get(H, 0, 0) +
+ gsl_matrix_get(H, 0, 1) * gsl_matrix_get(H, 1, 0) -
+ s*gsl_matrix_get(H, 0, 0) + t;
+ y = gsl_matrix_get(H, 1, 0) *
+ (gsl_matrix_get(H, 0, 0) + gsl_matrix_get(H, 1, 1) - s);
+ z = gsl_matrix_get(H, 1, 0) * gsl_matrix_get(H, 2, 1);
+
+ v = gsl_vector_subvector(w->v, 0, 3);
+
+ for (i = 0; i < N - 2; ++i)
+ {
+ gsl_vector_set(&v.vector, 0, x);
+ gsl_vector_set(&v.vector, 1, y);
+ gsl_vector_set(&v.vector, 2, z);
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ /* q = max(1, i - 1) */
+ q = (1 > ((int)i - 1)) ? 0 : (i - 1);
+
+ /* apply left householder matrix (I - tau_i v v') to H */
+ m = gsl_matrix_submatrix(H, i, q, 3, N - q);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ /* r = min(i + 3, N - 1) */
+ r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
+
+ /* apply right householder matrix (I - tau_i v v') to H */
+ m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+ x = gsl_matrix_get(H, i + 1, i);
+ y = gsl_matrix_get(H, i + 2, i);
+ if (i < (N - 3))
+ {
+ z = gsl_matrix_get(H, i + 3, i);
+ }
+ }
+
+ v = gsl_vector_subvector(w->v, 0, 2);
+ gsl_vector_set(&v.vector, 0, x);
+ gsl_vector_set(&v.vector, 1, y);
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+ return GSL_SUCCESS;
+}
+
+/*
+get_2b2_eigenvalues()
+ Compute the eigenvalues of a 2x2 real matrix
+
+Inputs: A - 2x2 matrix
+ e1 - where to store eigenvalue 1
+ e2 - where to store eigenvalue 2
+*/
+
+static void
+get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+ gsl_complex * e2)
+{
+ double discr; /* discriminant of characteristic poly */
+ double a, b, c, d; /* matrix values */
+
+ a = gsl_matrix_get(A, 0, 0);
+ b = gsl_matrix_get(A, 0, 1);
+ c = gsl_matrix_get(A, 1, 0);
+ d = gsl_matrix_get(A, 1, 1);
+
+ discr = (a + d)*(a + d) - 4.0*(a*d - b*c);
+ if (discr < 0.0)
+ {
+ GSL_SET_REAL(e1, 0.5*(a + d));
+ GSL_SET_REAL(e2, 0.5*(a + d));
+
+ GSL_SET_IMAG(e1, 0.5*sqrt(-discr));
+ GSL_SET_IMAG(e2, -0.5*sqrt(-discr));
+ }
+ else
+ {
+ GSL_SET_REAL(e1, 0.5*(a + d + sqrt(discr)));
+ GSL_SET_REAL(e2, 0.5*(a + d - sqrt(discr)));
+
+ GSL_SET_IMAG(e1, 0.0);
+ GSL_SET_IMAG(e2, 0.0);
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/gsl.unsymm.patch ./gsl.unsymm.patch
--- /home/palken/tmp2/gsl-1.8/gsl.unsymm.patch 1969-12-31 17:00:00.000000000 -0700
+++ ./gsl.unsymm.patch 2006-05-17 12:07:19.000000000 -0600
@@ -0,0 +1,571 @@
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
+--- /home/palken/tmp2/gsl-1.8/eigen/balance.c 1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/balance.c 2006-05-17 11:57:52.000000000 -0600
+@@ -0,0 +1,117 @@
++/* eigen/balance.c
++ *
++ * Copyright (C) 2006 Patrick Alken
++ *
++ * This program 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) any later version.
++ *
++ * This program 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
++ * General Public License for more details.
++ *
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#include <stdlib.h>
++#include <math.h>
++
++#include "balance.h"
++
++/*
++balance_matrix()
++ Balance a given matrix by applying a diagonal matrix
++similarity transformation so that the rows and columns
++of the new matrix have norms which are the same order of
++magnitude. This is necessary for the unsymmetric eigenvalue
++problem since the calculation can become numerically unstable
++for unbalanced matrices.
++
++See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
++
++and
++
++Numerical Recipes in C, Section 11.5
++*/
++
++void
++balance_matrix(gsl_matrix * A)
++{
++ double row_norm,
++ col_norm;
++ int not_converged;
++ const size_t N = A->size1;
++
++ not_converged = 1;
++
++ while (not_converged)
++ {
++ size_t i, j;
++ double g, f, s;
++
++ not_converged = 0;
++
++ for (i = 0; i < N; ++i)
++ {
++ row_norm = 0.0;
++ col_norm = 0.0;
++
++ for (j = 0; j < N; ++j)
++ {
++ if (j != i)
++ {
++ col_norm += fabs(gsl_matrix_get(A, j, i));
++ row_norm += fabs(gsl_matrix_get(A, i, j));
++ }
++ }
++
++ if ((col_norm == 0.0) || (row_norm == 0.0))
++ {
++ continue;
++ }
++
++ g = row_norm / FLOAT_RADIX;
++ f = 1.0;
++ s = col_norm + row_norm;
++
++ /*
++ * find the integer power of the machine radix which
++ * comes closest to balancing the matrix
++ */
++ while (col_norm < g)
++ {
++ f *= FLOAT_RADIX;
++ col_norm *= FLOAT_RADIX_SQ;
++ }
++
++ g = row_norm * FLOAT_RADIX;
++
++ while (col_norm > g)
++ {
++ f /= FLOAT_RADIX;
++ col_norm /= FLOAT_RADIX_SQ;
++ }
++
++ if ((row_norm + col_norm) < 0.95 * s * f)
++ {
++ not_converged = 1;
++
++ g = 1.0 / f;
++
++ /* apply similarity transformation */
++ for (j = 0; j < N; ++j)
++ {
++ gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) * g);
++ }
++ for (j = 0; j < N; ++j)
++ {
++ gsl_matrix_set(A, j, i, gsl_matrix_get(A, j, i) * f);
++ }
++ }
++ }
++ }
++}
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
+--- /home/palken/tmp2/gsl-1.8/eigen/balance.h 1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/balance.h 2006-05-17 10:49:00.000000000 -0600
+@@ -0,0 +1,32 @@
++/* eigen/balance.h
++ *
++ * Copyright (C) 2006 Patrick Alken
++ *
++ * This program 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) any later version.
++ *
++ * This program 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
++ * General Public License for more details.
++ *
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#ifndef __GSL_BALANCE_H__
++#define __GSL_BALANCE_H__
++
++#include <gsl/gsl_matrix.h>
++
++/* floating point radix */
++#define FLOAT_RADIX 2.0
++
++#define FLOAT_RADIX_SQ (FLOAT_RADIX * FLOAT_RADIX)
++
++void balance_matrix(gsl_matrix * A);
++
++#endif /* __GSL_BALANCE_H__ */
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
+--- /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h 2005-06-26 07:25:34.000000000 -0600
++++ ./eigen/gsl_eigen.h 2006-05-17 12:02:26.000000000 -0600
+@@ -59,6 +59,16 @@
+
+ typedef struct {
+ size_t size;
++ gsl_vector *v; /* temporary vector */
++} gsl_eigen_unsymm_workspace;
++
++gsl_eigen_unsymm_workspace * gsl_eigen_unsymm_alloc (const size_t n);
++void gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w);
++int gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
++ gsl_eigen_unsymm_workspace * w);
++
++typedef struct {
++ size_t size;
+ double * d;
+ double * sd;
+ double * tau;
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
+--- /home/palken/tmp2/gsl-1.8/eigen/Makefile.am 2004-09-11 07:45:45.000000000 -0600
++++ ./eigen/Makefile.am 2006-05-17 12:00:33.000000000 -0600
+@@ -3,7 +3,7 @@
+ check_PROGRAMS = test
+
+ pkginclude_HEADERS = gsl_eigen.h
+-libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c herm.c hermv.c sort.c
++libgsleigen_la_SOURCES = jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c balance.c
+ INCLUDES= -I$(top_builddir)
+
+ noinst_HEADERS = qrstep.c
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
+--- /home/palken/tmp2/gsl-1.8/eigen/unsymm.c 1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/unsymm.c 2006-05-17 11:57:57.000000000 -0600
+@@ -0,0 +1,378 @@
++/* eigen/unsymm.c
++ *
++ * Copyright (C) 2006 Patrick Alken
++ *
++ * This program 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) any later version.
++ *
++ * This program 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
++ * General Public License for more details.
++ *
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#include <stdlib.h>
++#include <gsl/gsl_eigen.h>
++#include <gsl/gsl_linalg.h>
++#include <gsl/gsl_math.h>
++#include <gsl/gsl_blas.h>
++#include <gsl/gsl_vector.h>
++#include <gsl/gsl_vector_complex.h>
++#include <gsl/gsl_matrix.h>
++
++#include "balance.h"
++
++/* maximum iterations before failure is reported */
++#define UNSYMM_MAX_ITERATIONS 100
++
++/*
++ * This module computes the eigenvalues of a real unsymmetric
++ * matrix, using the QR decomposition.
++ *
++ * See Golub & Van Loan, "Matrix Computations" (3rd ed),
++ * algorithm 7.5.2
++ */
++
++inline static void zero_subdiag_small_elements(gsl_matrix * A);
++static inline int francis_qrstep(gsl_matrix * H,
++ gsl_eigen_unsymm_workspace * w);
++static void get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
++ gsl_complex * e2);
++
++/*
++gsl_eigen_unsymm_alloc()
++
++Allocate a workspace for solving the unsymmetric
++eigenvalue/eigenvector problem
++
++Inputs: n - size of matrix
++
++Return: pointer to workspace
++*/
++
++gsl_eigen_unsymm_workspace *
++gsl_eigen_unsymm_alloc(const size_t n)
++{
++ gsl_eigen_unsymm_workspace *w;
++
++ if (n == 0)
++ {
++ GSL_ERROR_NULL ("matrix dimension must be positive integer",
++ GSL_EINVAL);
++ }
++
++ w = ((gsl_eigen_unsymm_workspace *)
++ malloc (sizeof (gsl_eigen_unsymm_workspace)));
++
++ if (w == 0)
++ {
++ GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
++ }
++
++ w->size = n;
++ w->v = gsl_vector_alloc(3);
++
++ return (w);
++}
++
++void
++gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
++{
++ gsl_vector_free(w->v);
++
++ free(w);
++}
++
++/*
++gsl_eigen_unsymm()
++
++Solve the unsymmetric eigenvalue problem
++
++A x = \lambda x
++
++for the eigenvalues \lambda using algorithm 7.5.2 of
++Golub & Van Loan, "Matrix Computations" (3rd ed)
++
++Inputs: A - matrix
++ eval - where to store eigenvalues
++ w - workspace
++
++Notes: the matrix A is destroyed by this routine
++*/
++
++int
++gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
++ gsl_eigen_unsymm_workspace * w)
++{
++ /* check matrix and vector sizes */
++
++ if (A->size1 != A->size2)
++ {
++ GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
++ }
++ else if (eval->size != A->size1)
++ {
++ GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
++ }
++ else
++ {
++ size_t N;
++ size_t q;
++ gsl_complex lambda1, lambda2; /* eigenvalues */
++ gsl_matrix_view m;
++ size_t nev;
++ size_t nit;
++
++ N = A->size1;
++
++ /* special cases */
++ if (N == 1)
++ {
++ GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(A, 0, 0), 0.0);
++ gsl_vector_complex_set(eval, 0, lambda1);
++ return GSL_SUCCESS;
++ }
++
++ if (N == 2)
++ {
++ /*
++ * The 2x2 case is special since the matrix is already
++ * in upper quasi-triangular form so no Schur decomposition
++ * is necessary
++ */
++ get_2b2_eigenvalues(A, &lambda1, &lambda2);
++ gsl_vector_complex_set(eval, 0, lambda1);
++ gsl_vector_complex_set(eval, 1, lambda2);
++ return GSL_SUCCESS;
++ }
++
++ nev = 0;
++
++ /* balance the matrix */
++ balance_matrix(A);
++
++ /* compute the Hessenberg reduction of A */
++ gsl_linalg_hessenberg(A);
++
++ /* set small elements on the subdiagonal to 0 */
++ zero_subdiag_small_elements(A);
++
++ m = gsl_matrix_submatrix(A, 0, 0, N, N);
++
++ nit = 0;
++ while ((N > 2) && (nit++ < UNSYMM_MAX_ITERATIONS))
++ {
++ francis_qrstep(&m.matrix, w);
++ zero_subdiag_small_elements(&m.matrix);
++
++ q = N - 1;
++ while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
++ {
++ --q;
++ }
++
++ if (q == (N - 1))
++ {
++ /*
++ * the (N, N - 1) element of the matrix is 0 -
++ * m_{NN} is a real eigenvalue
++ */
++ GSL_SET_COMPLEX(&lambda1,
++ gsl_matrix_get(&m.matrix, q, q), 0.0);
++ gsl_vector_complex_set(eval, nev++, lambda1);
++
++ --N;
++ m = gsl_matrix_submatrix(A, 0, 0, N, N);
++ }
++ else if (q == (N - 2))
++ {
++ /*
++ * The bottom right 2x2 block of m is an eigenvalue
++ * system
++ */
++
++ m = gsl_matrix_submatrix(A, q, q, 2, 2);
++ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
++ gsl_vector_complex_set(eval, nev++, lambda1);
++ gsl_vector_complex_set(eval, nev++, lambda2);
++
++ N -= 2;
++ m = gsl_matrix_submatrix(A, 0, 0, N, N);
++ }
++ }
++
++ if (N == 1)
++ {
++ GSL_SET_COMPLEX(&lambda1,
++ gsl_matrix_get(&m.matrix, 0, 0), 0.0);
++ gsl_vector_complex_set(eval, nev++, lambda1);
++ }
++ else if (N == 2)
++ {
++ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
++ gsl_vector_complex_set(eval, nev++, lambda1);
++ gsl_vector_complex_set(eval, nev++, lambda2);
++ }
++
++ return GSL_SUCCESS;
++ }
++} /* gsl_eigen_unsymm() */
++
++/********************************************
++ * INTERNAL ROUTINES *
++ ********************************************/
++
++/*
++zero_subdiag_small_elements()
++ Sets to zero all elements on the subdiaganal of a matrix A
++which satisfy
++
++|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
++*/
++
++inline static void
++zero_subdiag_small_elements(gsl_matrix * A)
++{
++ const size_t N = A->size1;
++ size_t i;
++ double dpel = gsl_matrix_get(A, 0, 0);
++
++ for (i = 1; i < N; ++i)
++ {
++ double sel = gsl_matrix_get(A, i, i - 1);
++ double del = gsl_matrix_get(A, i, i);
++
++ if (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel)))
++ {
++ gsl_matrix_set(A, i, i - 1, 0.0);
++ }
++
++ dpel = del;
++ }
++}
++
++/*
++francis_qrstep()
++ Perform a Francis QR step.
++
++See Golub & Van Loan, "Matrix Computations" (3rd ed),
++algorithm 7.5.1
++
++Inputs: H - unreduced upper Hessenberg matrix
++ w - workspace
++*/
++
++static inline int
++francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w)
++{
++ const size_t N = H->size1;
++ double s, t, x, y, z;
++ size_t i;
++ gsl_vector_view v;
++ gsl_matrix_view m;
++ double tau_i;
++ size_t q, r;
++
++ /* s = a_1 + a_2 = H_{mm} + H_{nn} */
++ s = gsl_matrix_get(H, N - 2, N - 2) + gsl_matrix_get(H, N - 1, N - 1);
++
++ /* t = a_1 * a_2 = H_{mm} * H_{nn} - H_{mn} * H_{nm} */
++ t = gsl_matrix_get(H, N - 2, N - 2) * gsl_matrix_get(H, N - 1, N - 1) -
++ gsl_matrix_get(H, N - 2, N - 1) * gsl_matrix_get(H, N - 1, N - 2);
++ x = gsl_matrix_get(H, 0, 0) * gsl_matrix_get(H, 0, 0) +
++ gsl_matrix_get(H, 0, 1) * gsl_matrix_get(H, 1, 0) -
++ s*gsl_matrix_get(H, 0, 0) + t;
++ y = gsl_matrix_get(H, 1, 0) *
++ (gsl_matrix_get(H, 0, 0) + gsl_matrix_get(H, 1, 1) - s);
++ z = gsl_matrix_get(H, 1, 0) * gsl_matrix_get(H, 2, 1);
++
++ v = gsl_vector_subvector(w->v, 0, 3);
++
++ for (i = 0; i < N - 2; ++i)
++ {
++ gsl_vector_set(&v.vector, 0, x);
++ gsl_vector_set(&v.vector, 1, y);
++ gsl_vector_set(&v.vector, 2, z);
++ tau_i = gsl_linalg_householder_transform(&v.vector);
++
++ /* q = max(1, i - 1) */
++ q = (1 > ((int)i - 1)) ? 0 : (i - 1);
++
++ /* apply left householder matrix (I - tau_i v v') to H */
++ m = gsl_matrix_submatrix(H, i, q, 3, N - q);
++ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
++
++ /* r = min(i + 3, N - 1) */
++ r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
++
++ /* apply right householder matrix (I - tau_i v v') to H */
++ m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
++ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
++
++ x = gsl_matrix_get(H, i + 1, i);
++ y = gsl_matrix_get(H, i + 2, i);
++ if (i < (N - 3))
++ {
++ z = gsl_matrix_get(H, i + 3, i);
++ }
++ }
++
++ v = gsl_vector_subvector(w->v, 0, 2);
++ gsl_vector_set(&v.vector, 0, x);
++ gsl_vector_set(&v.vector, 1, y);
++ tau_i = gsl_linalg_householder_transform(&v.vector);
++
++ m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
++ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
++
++ m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
++ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
++
++ return GSL_SUCCESS;
++}
++
++/*
++get_2b2_eigenvalues()
++ Compute the eigenvalues of a 2x2 real matrix
++
++Inputs: A - 2x2 matrix
++ e1 - where to store eigenvalue 1
++ e2 - where to store eigenvalue 2
++*/
++
++static void
++get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
++ gsl_complex * e2)
++{
++ double discr; /* discriminant of characteristic poly */
++ double a, b, c, d; /* matrix values */
++
++ a = gsl_matrix_get(A, 0, 0);
++ b = gsl_matrix_get(A, 0, 1);
++ c = gsl_matrix_get(A, 1, 0);
++ d = gsl_matrix_get(A, 1, 1);
++
++ discr = (a + d)*(a + d) - 4.0*(a*d - b*c);
++ if (discr < 0.0)
++ {
++ GSL_SET_REAL(e1, 0.5*(a + d));
++ GSL_SET_REAL(e2, 0.5*(a + d));
++
++ GSL_SET_IMAG(e1, 0.5*sqrt(-discr));
++ GSL_SET_IMAG(e2, -0.5*sqrt(-discr));
++ }
++ else
++ {
++ GSL_SET_REAL(e1, 0.5*(a + d + sqrt(discr)));
++ GSL_SET_REAL(e2, 0.5*(a + d - sqrt(discr)));
++
++ GSL_SET_IMAG(e1, 0.0);
++ GSL_SET_IMAG(e2, 0.0);
++ }
++}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h ./linalg/gsl_linalg.h
--- /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h 2005-11-09 14:13:14.000000000 -0700
+++ ./linalg/gsl_linalg.h 2006-05-17 12:01:43.000000000 -0600
@@ -113,6 +113,10 @@
const gsl_vector_complex * v,
gsl_vector_complex * w);
+/* Hessenberg reduction */
+
+int gsl_linalg_hessenberg (gsl_matrix * A);
+
/* Singular Value Decomposition
* exceptions:
diff -urN /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c ./linalg/hessenberg.c
--- /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c 1969-12-31 17:00:00.000000000 -0700
+++ ./linalg/hessenberg.c 2006-05-16 15:11:12.000000000 -0600
@@ -0,0 +1,76 @@
+/* linalg/hessenberg.c
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * This program 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) any later version.
+ *
+ * This program 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+/*
+gsl_linalg_hessenberg()
+ Compute the Householder reduction to Hessenberg form of a
+square matrix A.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm
+7.4.2
+*/
+
+int
+gsl_linalg_hessenberg(gsl_matrix * A)
+{
+ if (A->size1 != A->size2)
+ {
+ GSL_ERROR ("Hessenberg reduction requires square matrix",
+ GSL_ENOTSQR);
+ }
+ else
+ {
+ const size_t N = A->size1;
+ size_t i, j;
+ gsl_vector_view v;
+ gsl_vector *v_copy;
+ gsl_matrix_view m;
+ double tau_i; /* beta in algorithm 7.4.2 */
+
+ v_copy = gsl_vector_alloc(N);
+
+ for (i = 0; i < N - 2; ++i)
+ {
+ /* make a copy of A(i + 1:n, i) */
+
+ v = gsl_matrix_column(A, i);
+ gsl_vector_memcpy(v_copy, &v.vector);
+ v = gsl_vector_subvector(v_copy, i + 1, N - (i + 1));
+
+ /* compute householder transformation of A(i+1:n,i) */
+ tau_i = gsl_linalg_householder_transform(&v.vector);
+
+ /* apply left householder matrix (I - tau_i v v') to A */
+ m = gsl_matrix_submatrix(A, i + 1, i, N - (i + 1), N - i);
+ gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+ /* apply right householder matrix (I - tau_i v v') to A */
+ m = gsl_matrix_submatrix(A, 0, i + 1, N, N - (i + 1));
+ gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+ }
+
+ gsl_vector_free(v_copy);
+
+ return GSL_SUCCESS;
+ }
+}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/Makefile.am ./linalg/Makefile.am
--- /home/palken/tmp2/gsl-1.8/linalg/Makefile.am 2004-09-13 12:17:04.000000000 -0600
+++ ./linalg/Makefile.am 2006-05-17 12:01:11.000000000 -0600
@@ -4,7 +4,7 @@
INCLUDES= -I$(top_builddir)
-libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
+libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c hessenberg.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
noinst_HEADERS = givens.c apply_givens.c svdstep.c tridiag.h