* Re: unsymmetric eigenvalue problem
@ 2006-05-18 21:36 Patrick.Alken
0 siblings, 0 replies; 16+ messages in thread
From: Patrick.Alken @ 2006-05-18 21:36 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 809 bytes --]
Hi again,
There was a nasty underflow problem in the last patch I submitted, which was a
result of the gsl_linalg_householder_transform routine. In some steps of the
Francis QR algorithm it is necessary to compute a householder vector with
small elements, and the householder_transform wasn't set up to do this. I have
attached a new patch which adds some code to gsl_linalg_householder_transform
to basically scale the input vector v -> v / |v|. This seems to correct the
underflow problem.
On another note, I haven't gotten the convergence criteria completely right
yet for the eigenvalue solver, so I need to study Golub & Van Loan a bit more -
they aren't completely clear about the convergence test. But the patch I
submitted should work a lot better than the previous one.
Thanks,
Patrick Alken
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: gsl.unsymm.patch --]
[-- Type: text/x-patch; name="gsl.unsymm.patch", Size: 21797 bytes --]
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-16 16:38:14.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-18 15:04:48.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-18 15:20:58.000000000 -0600
@@ -0,0 +1,400 @@
+/* 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 <math.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 "gsl_eigen.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: On output, A contains the upper block triangular Schur
+ decomposition.
+*/
+
+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 p, q;
+ gsl_complex lambda1, lambda2; /* eigenvalues */
+ gsl_matrix_view m;
+ size_t nev; /* number of eigenvalues found so far */
+ size_t nit; /* number of iterations */
+ size_t i,j,k;
+
+ 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);
+
+ /* Find the largest q such that A_{q,q-1} is 0 */
+ q = N - 1;
+ while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
+ {
+ --q;
+ }
+
+ p = 0;
+ while ((p < (N - 1)) && (gsl_matrix_get(A, p + 1, p) != 0.0))
+ {
+ ++p;
+ }
+
+ 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 (nit > UNSYMM_MAX_ITERATIONS)
+ {
+ GSL_ERROR("maximum iterations exceeded", GSL_EMAXITER);
+ }
+
+ 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);
+ }
+
+ printf("nit = %d\n", nit);
+
+ 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);
+
+ if (tau_i != 0.0)
+ {
+
+ /* 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-16 10:36:55.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/householder.c ./linalg/householder.c
--- /home/palken/tmp2/gsl-1.8/linalg/householder.c 2005-06-26 07:25:35.000000000 -0600
+++ ./linalg/householder.c 2006-05-18 14:56:00.000000000 -0600
@@ -41,10 +41,25 @@
else
{
double alpha, beta, tau ;
+ size_t i;
+ double xnorm;
gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
- double xnorm = gsl_blas_dnrm2 (&x.vector);
+ double vnorm = gsl_blas_dnrm2 (v);
+
+ if (vnorm == 0)
+ {
+ return 0.0;
+ }
+
+ /* v -> v / |v| to prevent overflow/underflow problems */
+ for (i = 0; i < n; ++i)
+ {
+ gsl_vector_set(v, i, gsl_vector_get(v, i) / vnorm);
+ }
+
+ xnorm = gsl_blas_dnrm2 (&x.vector);
if (xnorm == 0)
{
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-18 15:06:29.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
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-26 10:43 Patrick.Alken
@ 2006-05-27 21:37 ` Brian Gough
0 siblings, 0 replies; 16+ messages in thread
From: Brian Gough @ 2006-05-27 21:37 UTC (permalink / raw)
To: Patrick.Alken; +Cc: gsl-discuss
Patrick.Alken@colorado.edu writes:
> Thanks for making the test program - its a very good idea! I ran it on a variety
> of cases and found that there are some convergence issues in my code. Most of
> these can be fixed by increasing UNSYMM_MAX_ITERATIONS in unsymm.c, but on the
> other hand I don't think it should be necessary to have so many QR iterations
> for convergence.
There is a Lapack Working Note (LAWN) archive on netlib, all good
stuff -- there might be some extra information about the
implementation there.
--
Brian Gough
Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
@ 2006-05-26 10:43 Patrick.Alken
2006-05-27 21:37 ` Brian Gough
0 siblings, 1 reply; 16+ messages in thread
From: Patrick.Alken @ 2006-05-26 10:43 UTC (permalink / raw)
To: gsl-discuss
Thanks for making the test program - its a very good idea! I ran it on a variety
of cases and found that there are some convergence issues in my code. Most of
these can be fixed by increasing UNSYMM_MAX_ITERATIONS in unsymm.c, but on the
other hand I don't think it should be necessary to have so many QR iterations
for convergence.
After a cursory look at LAPACK it seems to me based on a comment in dhseqr that
it reports failure after 30*N total unsuccessful iterations. This seems like a
lot to me, my method reports failure after 30 iterations without finding an
eigenvalue (I don't keep track of total iterations).
Another possibility is the actual algorithm. LAPACK uses the multi-shift QR
algorithm published in 1989, whereas I implemented the original double-shift
algorithm from the 1960s which was described in Golub and Van Loan. I might
look into the multi-shift method when I get some time and see how its
convergence compares with the double-shift method. I'll be away for a week or
so, so I won't have a chance to look at this for a little while.
Thanks,
Patrick Alken
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-22 22:26 ` Brian Gough
@ 2006-05-23 18:35 ` Brian Gough
0 siblings, 0 replies; 16+ messages in thread
From: Brian Gough @ 2006-05-23 18:35 UTC (permalink / raw)
To: Patrick Alken, gsl-discuss
You can find my test program here
http://www.briangough.ukfsn.org/unsymm/
It's not totally finished, example output below.
--
Brian Gough
Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/
$ ./a.out -n 3 -i -l -1 -u 1
testing N=3 incrementally on element range [-1,1]
==== CASE 17 ==========================
A = [
0, 1, 0;
-1, -1, -1;
-1, -1, -1;
]
=== eval - lapack dgees ===
-1.813985266369727769e-32 0.000000000000000000e+00
-9.999999999999996669e-01 1.950115636758225308e-08
-9.999999999999996669e-01 -1.950115636758225308e-08
=== eval - alken ===
-6.631056343480628004e-33 0.000000000000000000e+00
-9.999999999999997780e-01 1.633303231307825499e-08
-9.999999999999997780e-01 -1.633303231307825499e-08
max abserr = 3.16812e-09 max relerr = 0.162458
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-20 23:58 ` Patrick Alken
2006-05-21 5:17 ` James Bergstra
2006-05-21 10:13 ` Jochen Küpper
@ 2006-05-22 22:26 ` Brian Gough
2006-05-23 18:35 ` Brian Gough
2 siblings, 1 reply; 16+ messages in thread
From: Brian Gough @ 2006-05-22 22:26 UTC (permalink / raw)
To: Patrick Alken; +Cc: gsl-discuss
Patrick Alken writes:
> I haven't gotten much feedback so far so I hope others will find
> this useful..
I've tested the earlier versions with valgrind -- all ok.
I am writing a program to do a direct comparison with the LAPACK DGEES
routine -- I'm about halfway through that.
I'll use it to do some independent exhaustive tests for small integer
matrices and leave it running on large random matrices to compare the
accuracy.
--
Brian Gough
Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-18 16:24 ` Linas Vepstas
@ 2006-05-22 21:37 ` Brian Gough
0 siblings, 0 replies; 16+ messages in thread
From: Brian Gough @ 2006-05-22 21:37 UTC (permalink / raw)
To: Linas Vepstas; +Cc: Patrick.Alken, gsl-discuss
Linas Vepstas writes:
> > + f *= FLOAT_RADIX;
>
> Is this right? isn't there a "DOUBLE_RADIX" ?
I think the C89 standard defines FLT_RADIX only.
--
Brian Gough
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-20 23:58 ` Patrick Alken
2006-05-21 5:17 ` James Bergstra
@ 2006-05-21 10:13 ` Jochen Küpper
2006-05-22 22:26 ` Brian Gough
2 siblings, 0 replies; 16+ messages in thread
From: Jochen Küpper @ 2006-05-21 10:13 UTC (permalink / raw)
To: gsl-discuss
Patrick Alken <patrick.alken@colorado.edu> writes:
> and compared the results with octave and find perfect agreement.
Good;)
> In terms of speed, my code runs slower than octave's eigenvalue
> solver (which I think is essentially the fortran code from LAPACK)
> but I only notice the slowdown at around N = 500 matrices.
You are suggesting that your code is roughly as fast as LAPACK up to
N=500? Not knowing anything about the internals of octave, but they
might be switching over to a different LAPACK routine/algorithm for
larger matrices, not?
> There is certainly room for improvement here but I believe the
> priority for now is to get an eigenvalue solver into GSL.
Sounds good.
> I haven't gotten much feedback so far so I hope others will find
> this useful..
It definitely closes one hole in GSL, and as such I find it very
important. Personally I only have to deal with symmetric/hermitian
matrices (for now), so there is not much to say about your code from
my point of view.
I hope to see this code in GSL soon, and I hope that /time/ will
improve other parts of GSL as well.
Greetings,
Jochen
--
Einigkeit und Recht und Freiheit http://www.Jochen-Kuepper.de
Liberté, Égalité, Fraternité GnuPG key: CC1B0B4D
(Part 3 you find in my messages before fall 2003.)
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-20 23:58 ` Patrick Alken
@ 2006-05-21 5:17 ` James Bergstra
2006-05-21 10:13 ` Jochen Küpper
2006-05-22 22:26 ` Brian Gough
2 siblings, 0 replies; 16+ messages in thread
From: James Bergstra @ 2006-05-21 5:17 UTC (permalink / raw)
To: Patrick Alken; +Cc: gsl-discuss
On Sat, May 20, 2006 at 05:58:49PM -0600, Patrick Alken wrote:
...
> I haven't gotten much feedback so far so I hope others will find
> this useful..
Yeah, this discussion list isn't a warm fuzzy support group. No comments could
mean anything from impeccable to irrelevent. While personally, I'm still waiting
[passively] for a LAPACK wrapper / tutorial, I think you've identified an
important gap in linalg and I hope that this code finds its way into a future
release.
I was hoping to pass on minor coding tricks since you were so brave as to post
your entire patch... but it looks good to me!
'night
--
James Bergstra
http://www-etud.iro.umontreal.ca/~bergstrj
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-18 23:46 Patrick.Alken
@ 2006-05-20 23:58 ` Patrick Alken
2006-05-21 5:17 ` James Bergstra
` (2 more replies)
0 siblings, 3 replies; 16+ messages in thread
From: Patrick Alken @ 2006-05-20 23:58 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 815 bytes --]
Hello all,
I have implemented some convergence tricks suggested by
Numerical Recipes and have attached the latest patch.
This should be the last patch for a while, since I have now
tested it on randomly generated 1000x1000 matrices and compared
the results with octave and find perfect agreement.
In terms of speed, my code runs slower than octave's eigenvalue
solver (which I think is essentially the fortran code from LAPACK)
but I only notice the slowdown at around N = 500 matrices. There
is certainly room for improvement here but I believe the priority
for now is to get an eigenvalue solver into GSL. I will continue
trying to optimize the code but it will involve studying the
LAPACK source for a while.
I haven't gotten much feedback so far so I hope others will find
this useful..
Patrick Alken
[-- Attachment #2: gsl.unsymm.patch --]
[-- Type: text/plain, Size: 25575 bytes --]
diff -urN /home/cosine/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
--- /home/cosine/tmp2/gsl-1.8/eigen/Makefile.am 2004-09-11 07:45:45.000000000 -0600
+++ ./eigen/Makefile.am 2006-05-20 13:41:10.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/cosine/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
--- /home/cosine/tmp2/gsl-1.8/eigen/balance.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.c 2006-05-20 13:41:10.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/cosine/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
--- /home/cosine/tmp2/gsl-1.8/eigen/balance.h 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.h 2006-05-20 13:41:10.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/cosine/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
--- /home/cosine/tmp2/gsl-1.8/eigen/gsl_eigen.h 2005-06-26 07:25:34.000000000 -0600
+++ ./eigen/gsl_eigen.h 2006-05-20 17:32:45.000000000 -0600
@@ -59,6 +59,17 @@
typedef struct {
size_t size;
+ gsl_vector *v2; /* temporary 2x1 vector */
+ gsl_vector *v3; /* temporary 3x1 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/cosine/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
--- /home/cosine/tmp2/gsl-1.8/eigen/unsymm.c 1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/unsymm.c 2006-05-20 17:44:52.000000000 -0600
@@ -0,0 +1,550 @@
+/* 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 <math.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 30
+
+/*
+ * 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 size_t schur_decomp(gsl_matrix * H, size_t top,
+ size_t bot, gsl_vector_complex * eval,
+ size_t evidx,
+ size_t * nit,
+ gsl_eigen_unsymm_workspace * w);
+inline static void zero_subdiag_small_elements(gsl_matrix * A);
+static inline int francis_qrstep(gsl_matrix * H,
+ gsl_eigen_unsymm_workspace * w,
+ double s, double t);
+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->v2 = gsl_vector_alloc(2);
+ w->v3 = gsl_vector_alloc(3);
+
+ return (w);
+}
+
+void
+gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
+{
+ gsl_vector_free(w->v2);
+
+ gsl_vector_free(w->v3);
+
+ 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: On output, A contains the upper block triangular Schur
+ decomposition.
+*/
+
+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;
+ gsl_complex lambda1, lambda2; /* eigenvalues */
+ 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;
+ }
+
+ /* 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);
+
+ nit = 0;
+
+ /*
+ * compute Schur decomposition of A and store eigenvalues
+ * into eval
+ */
+ schur_decomp(A, 0, N - 1, eval, 0, &nit, w);
+
+ if (nit > UNSYMM_MAX_ITERATIONS)
+ {
+ GSL_ERROR("maximum iterations exceeded", GSL_EMAXITER);
+ }
+
+ return GSL_SUCCESS;
+ }
+} /* gsl_eigen_unsymm() */
+
+/********************************************
+ * INTERNAL ROUTINES *
+ ********************************************/
+
+/*
+schur_decomp()
+ Compute the Schur decomposition of the submatrix of H
+starting from (top, top) to (bot, bot)
+
+Inputs: H - hessenberg matrix
+ top - top index
+ bot - bottom index
+ eval - where to store eigenvalues
+ evidx - index into eval
+ nit - running total of number of QR iterations since
+ we last found an eigenvalue
+ (must be initialized before calling this function)
+ w - workspace
+
+Return: number of eigenvalues found
+*/
+
+inline static size_t
+schur_decomp(gsl_matrix * H, size_t top, size_t bot,
+ gsl_vector_complex * eval, size_t evidx,
+ size_t * nit,
+ gsl_eigen_unsymm_workspace * w)
+{
+ gsl_matrix_view m;
+ size_t N; /* size of matrix */
+ double s, t; /* shifts */
+ size_t nev; /* number of eigenvalues found so far */
+ size_t q;
+ gsl_complex lambda1, /* eigenvalues */
+ lambda2;
+
+ N = bot - top + 1;
+
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(H, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, evidx, lambda1);
+ *nit = 0;
+ return 1;
+ }
+ else if (N == 2)
+ {
+ get_2b2_eigenvalues(H, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, evidx, lambda1);
+ gsl_vector_complex_set(eval, evidx + 1, lambda2);
+ *nit = 0;
+ return 2;
+ }
+
+ m = gsl_matrix_submatrix(H, top, top, N, N);
+
+ nev = 0;
+ while ((N > 2) && ((*nit)++ < UNSYMM_MAX_ITERATIONS))
+ {
+ if ((*nit == 10) || (*nit == 20))
+ {
+ /*
+ * We have gone 10 or 20 iterations without finding
+ * a new eigenvalue, try a new choice of shifts.
+ * See Numerical Recipes in C, eq 11.6.27
+ */
+ t = fabs(gsl_matrix_get(&m.matrix, N - 1, N - 2)) +
+ fabs(gsl_matrix_get(&m.matrix, N - 2, N - 3));
+ s = 1.5 * t;
+ t *= t;
+ }
+ else
+ {
+ /* s = a_1 + a_2 = m_{mm} + m_{nn} where m = n - 1 */
+ s = gsl_matrix_get(&m.matrix, N - 2, N - 2) +
+ gsl_matrix_get(&m.matrix, N - 1, N - 1);
+
+ /* t = a_1 * a_2 = m_{mm} * m_{nn} - m_{mn} * m_{nm} */
+ t = (gsl_matrix_get(&m.matrix, N - 2, N - 2) *
+ gsl_matrix_get(&m.matrix, N - 1, N - 1)) -
+ (gsl_matrix_get(&m.matrix, N - 2, N - 1) *
+ gsl_matrix_get(&m.matrix, N - 1, N - 2));
+ }
+
+ francis_qrstep(&m.matrix, w, s, t);
+ zero_subdiag_small_elements(&m.matrix);
+
+ /* Find the largest q such that A_{q,q-1} is 0 */
+ q = N - 1;
+ while ((q > 0) && (gsl_matrix_get(&m.matrix, 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, evidx + nev++, lambda1);
+ *nit = 0;
+
+ --N;
+ m = gsl_matrix_submatrix(&m.matrix, 0, 0, N, N);
+ }
+ else if (q == (N - 2))
+ {
+ gsl_matrix_view v;
+
+ /*
+ * The bottom right 2x2 block of m is an eigenvalue
+ * system
+ */
+
+ v = gsl_matrix_submatrix(&m.matrix, q, q, 2, 2);
+ get_2b2_eigenvalues(&v.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda1);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda2);
+ *nit = 0;
+
+ N -= 2;
+ m = gsl_matrix_submatrix(&m.matrix, 0, 0, N, N);
+ }
+ else if (q == 1)
+ {
+ /* the first matrix element is an eigenvalue */
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda1);
+ *nit = 0;
+
+ --N;
+ m = gsl_matrix_submatrix(&m.matrix, 1, 1, N, N);
+ }
+ else if (q == 2)
+ {
+ gsl_matrix_view v;
+
+ /* the upper left 2x2 block is an eigenvalue system */
+
+ v = gsl_matrix_submatrix(&m.matrix, 0, 0, 2, 2);
+ get_2b2_eigenvalues(&v.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda1);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda2);
+ *nit = 0;
+
+ N -= 2;
+ m = gsl_matrix_submatrix(&m.matrix, 2, 2, N, N);
+ }
+ else if (q != 0)
+ {
+ /*
+ * There is a zero element on the subdiagonal somewhere
+ * in the middle of the matrix - we can now operate
+ * separately on the two submatrices split by this
+ * element. q is the row index of the zero element.
+ */
+
+ /* operate on lower right (N - q)x(N - q) block first */
+ nev += schur_decomp(&m.matrix,
+ q,
+ N - 1,
+ eval,
+ evidx + nev,
+ nit,
+ w);
+
+ /* operate on upper left qxq block */
+ nev += schur_decomp(&m.matrix,
+ 0,
+ q - 1,
+ eval,
+ evidx + nev,
+ nit,
+ w);
+ N = 0;
+ }
+ }
+
+ if (N == 1)
+ {
+ GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda1);
+ *nit = 0;
+ }
+ else if (N == 2)
+ {
+ get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda1);
+ gsl_vector_complex_set(eval, evidx + nev++, lambda2);
+ *nit = 0;
+ }
+
+ return (nev);
+}
+
+/*
+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
+ s - sum of current shifts (a_1 + a_2)
+ t - product of shifts (a_1 * a_2)
+*/
+
+static inline int
+francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w,
+ double s, double t)
+{
+ const size_t N = H->size1;
+ double x, y, z;
+ double scale;
+ size_t i;
+ gsl_matrix_view m;
+ double tau_i;
+ size_t q, r;
+
+ 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);
+
+ scale = fabs(x) + fabs(y) + fabs(z);
+ if (scale != 0.0)
+ {
+ /* scale to prevent overflow or underflow */
+ x /= scale;
+ y /= scale;
+ z /= scale;
+ }
+
+ for (i = 0; i < N - 2; ++i)
+ {
+ gsl_vector_set(w->v3, 0, x);
+ gsl_vector_set(w->v3, 1, y);
+ gsl_vector_set(w->v3, 2, z);
+ tau_i = gsl_linalg_householder_transform(w->v3);
+
+ if (tau_i != 0.0)
+ {
+ /* 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, w->v3, &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, w->v3, &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);
+ }
+
+ scale = fabs(x) + fabs(y) + fabs(z);
+ if (scale != 0.0)
+ {
+ /* scale to prevent overflow or underflow */
+ x /= scale;
+ y /= scale;
+ z /= scale;
+ }
+ }
+
+ gsl_vector_set(w->v2, 0, x);
+ gsl_vector_set(w->v2, 1, y);
+ tau_i = gsl_linalg_householder_transform(w->v2);
+
+ m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
+ gsl_linalg_householder_hm(tau_i, w->v2, &m.matrix);
+
+ m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
+ gsl_linalg_householder_mh(tau_i, w->v2, &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/cosine/tmp2/gsl-1.8/linalg/Makefile.am ./linalg/Makefile.am
--- /home/cosine/tmp2/gsl-1.8/linalg/Makefile.am 2004-09-13 12:17:04.000000000 -0600
+++ ./linalg/Makefile.am 2006-05-20 13:41:10.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
diff -urN /home/cosine/tmp2/gsl-1.8/linalg/gsl_linalg.h ./linalg/gsl_linalg.h
--- /home/cosine/tmp2/gsl-1.8/linalg/gsl_linalg.h 2005-11-09 14:13:14.000000000 -0700
+++ ./linalg/gsl_linalg.h 2006-05-20 13:41:10.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/cosine/tmp2/gsl-1.8/linalg/hessenberg.c ./linalg/hessenberg.c
--- /home/cosine/tmp2/gsl-1.8/linalg/hessenberg.c 1969-12-31 17:00:00.000000000 -0700
+++ ./linalg/hessenberg.c 2006-05-20 13:41:10.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;
+ }
+}
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-19 15:13 ` C J Kenneth Tan -- OptimaNumerics
@ 2006-05-19 16:57 ` Patrick Alken
0 siblings, 0 replies; 16+ messages in thread
From: Patrick Alken @ 2006-05-19 16:57 UTC (permalink / raw)
To: gsl-discuss
Hi,
Code performance is very important of course. I haven't done
any benchmarks against LAPACK yet, since I want to make sure
the algorithm is stable first. It should be running in O(N^3)
time right now since thats what Golub & Van Loan say and I
implemented their algorithm pretty much as they outlined it.
I will optimize what I can once its working perfectly.
Patrick Alken
On Fri, May 19, 2006 at 12:51:04PM +0000, C J Kenneth Tan -- OptimaNumerics wrote:
> Patrick,
>
> On 2006-05-18 10:40 -0600 Patrick.Alken@colorado.edu wrote:
> > On Thu, May 18, 2006 at 11:21:33AM -0500, Linas Vepstas wrote:
> > > On Wed, May 17, 2006 at 02:45:10PM -0600, Patrick.Alken@colorado.edu wrote:
> > > >
> > > > I have recently needed to compute eigenvalues of unsymmetric
> > > > matrices,
> > >
> > > When I hit ths problem, someone recommended LAPACK, which
> > > I then used happily (routine DGEEV). This mailing list has
> > > seen occasional discusions about how to wrap/merge/whatever
> > > LAPACK routines with gsl (a technical problem is that LAPACK
> > > uses the fortran calling convention, making it ugly in C.)
> >
> > LAPACK I believe uses the same QR algorithm which I have been
> > implementing, so if I can make it as robust as LAPACK's
> > implementation it should be a far better alternative than
> > linking to fortran code.
>
> Are you taking into consideration performance of the code? How do you
> currently compare against NETLIB LAPACK?
>
>
> Kenneth Tan
> -----------------------------------------------------------------------
> C J Kenneth Tan, PhD
> OptimaNumerics Ltd Telephone: +44 798 941 7838
> E-mail: cjtan@OptimaNumerics.com Telephone: +44 207 099 4428
> Web: http://www.OptimaNumerics.com Facsimile: +44 207 100 4572
> -----------------------------------------------------------------------
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-18 16:43 Patrick.Alken
@ 2006-05-19 15:13 ` C J Kenneth Tan -- OptimaNumerics
2006-05-19 16:57 ` Patrick Alken
0 siblings, 1 reply; 16+ messages in thread
From: C J Kenneth Tan -- OptimaNumerics @ 2006-05-19 15:13 UTC (permalink / raw)
To: Patrick.Alken; +Cc: gsl-discuss
Patrick,
On 2006-05-18 10:40 -0600 Patrick.Alken@colorado.edu wrote:
> On Thu, May 18, 2006 at 11:21:33AM -0500, Linas Vepstas wrote:
> > On Wed, May 17, 2006 at 02:45:10PM -0600, Patrick.Alken@colorado.edu wrote:
> > >
> > > I have recently needed to compute eigenvalues of unsymmetric
> > > matrices,
> >
> > When I hit ths problem, someone recommended LAPACK, which
> > I then used happily (routine DGEEV). This mailing list has
> > seen occasional discusions about how to wrap/merge/whatever
> > LAPACK routines with gsl (a technical problem is that LAPACK
> > uses the fortran calling convention, making it ugly in C.)
>
> LAPACK I believe uses the same QR algorithm which I have been
> implementing, so if I can make it as robust as LAPACK's
> implementation it should be a far better alternative than
> linking to fortran code.
Are you taking into consideration performance of the code? How do you
currently compare against NETLIB LAPACK?
Kenneth Tan
-----------------------------------------------------------------------
C J Kenneth Tan, PhD
OptimaNumerics Ltd Telephone: +44 798 941 7838
E-mail: cjtan@OptimaNumerics.com Telephone: +44 207 099 4428
Web: http://www.OptimaNumerics.com Facsimile: +44 207 100 4572
-----------------------------------------------------------------------
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
@ 2006-05-18 23:46 Patrick.Alken
2006-05-20 23:58 ` Patrick Alken
0 siblings, 1 reply; 16+ messages in thread
From: Patrick.Alken @ 2006-05-18 23:46 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 386 bytes --]
Hi again,
I apologize for sending yet another patch, but this time the eigenvalue solver
is in much better shape. I have tested it on randomly generated 30x30 matrices
with excellent results - it converges every time within about 10-15 iterations
and matches the eigenvalues given by octave. I have attached the new patch
(remember to run automake before compiling).
Patrick Alken
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: gsl.unsymm.patch --]
[-- Type: text/x-patch; name="gsl.unsymm.patch", Size: 23541 bytes --]
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-16 16:38:14.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-18 17:08:34.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-18 17:07:18.000000000 -0600
@@ -0,0 +1,479 @@
+/* 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 <math.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,
+ double s, double t);
+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: On output, A contains the upper block triangular Schur
+ decomposition.
+*/
+
+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 p, q;
+ gsl_complex lambda1, lambda2; /* eigenvalues */
+ gsl_matrix_view m;
+ size_t nev; /* number of eigenvalues found so far */
+ size_t nit; /* number of iterations */
+ double s, t; /* shifts */
+
+ 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))
+ {
+ /*
+ * Calculate shift parameters for the QR step
+ */
+
+ if ((nit % 10) == 0)
+ {
+ /*
+ * adjust the shifts to accelerate convergence -
+ * see Numerical Recipes in C, eq 11.6.27
+ */
+ s = 1.5*(fabs(gsl_matrix_get(&m.matrix, N - 1, N - 2)) +
+ fabs(gsl_matrix_get(&m.matrix, N - 2, N - 3)));
+ t = fabs(gsl_matrix_get(&m.matrix, N - 1, N - 2)) +
+ fabs(gsl_matrix_get(&m.matrix, N - 2, N - 3));
+ t *= t;
+ }
+ else
+ {
+ /* s = a_1 + a_2 = m_{mm} + m_{nn} */
+ s = gsl_matrix_get(&m.matrix, N - 2, N - 2) +
+ gsl_matrix_get(&m.matrix, N - 1, N - 1);
+
+ /* t = a_1 * a_2 = m_{mm} * m_{nn} - m_{mn} * m_{nm} */
+ t = (gsl_matrix_get(&m.matrix, N - 2, N - 2) *
+ gsl_matrix_get(&m.matrix, N - 1, N - 1)) -
+ (gsl_matrix_get(&m.matrix, N - 2, N - 1) *
+ gsl_matrix_get(&m.matrix, N - 1, N - 2));
+ }
+
+ francis_qrstep(&m.matrix, w, s, t);
+ zero_subdiag_small_elements(&m.matrix);
+
+ /* Find the largest q such that A_{q,q-1} is 0 */
+ q = N - 1;
+ while ((q > 0) && (gsl_matrix_get(&m.matrix, 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(&m.matrix, 0, 0, N, N);
+ }
+ else if (q == (N - 2))
+ {
+ gsl_matrix_view v;
+
+ /*
+ * The bottom right 2x2 block of m is an eigenvalue
+ * system
+ */
+
+ v = gsl_matrix_submatrix(&m.matrix, q, q, 2, 2);
+ get_2b2_eigenvalues(&v.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+
+ N -= 2;
+ m = gsl_matrix_submatrix(&m.matrix, 0, 0, N, N);
+ }
+
+ /* Find the smallest p such that A_{p+1, p} is 0 */
+ p = 0;
+ while ((p < (N - 1)) && (gsl_matrix_get(&m.matrix, p + 1, p) != 0.0))
+ {
+ ++p;
+ }
+
+ if (p == 0)
+ {
+ /* The first matrix element is an eigenvalue */
+ GSL_SET_COMPLEX(&lambda1,
+ gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+
+ --N;
+ if (N > 0)
+ {
+ m = gsl_matrix_submatrix(&m.matrix, 1, 1, N, N);
+ }
+ }
+ else if (p == 1)
+ {
+ /* The upper left 2x2 block is an eigenvalue system */
+ gsl_matrix_view v;
+
+ v = gsl_matrix_submatrix(&m.matrix, 0, 0, 2, 2);
+ get_2b2_eigenvalues(&v.matrix, &lambda1, &lambda2);
+ gsl_vector_complex_set(eval, nev++, lambda1);
+ gsl_vector_complex_set(eval, nev++, lambda2);
+
+ N -= 2;
+ if (N > 0)
+ {
+ m = gsl_matrix_submatrix(&m.matrix, 2, 2, N, N);
+ }
+ }
+ }
+
+ if (nit > UNSYMM_MAX_ITERATIONS)
+ {
+ GSL_ERROR("maximum iterations exceeded", GSL_EMAXITER);
+ }
+
+ 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);
+ }
+
+ printf("nit = %d\n", nit);
+
+ 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
+ s - sum of current shifts (a_1 + a_2)
+ t - product of shifts (a_1 * a_2)
+*/
+
+static inline int
+francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w,
+ double s, double t)
+{
+ const size_t N = H->size1;
+ double x, y, z;
+ double scale;
+ size_t i;
+ gsl_vector_view v;
+ gsl_matrix_view m;
+ double tau_i;
+ size_t q, r;
+
+ 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);
+
+ scale = fabs(x) + fabs(y) + fabs(z);
+ if (scale != 0.0)
+ {
+ /* scale to prevent overflow or underflow */
+ x /= scale;
+ y /= scale;
+ z /= scale;
+ }
+
+ 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);
+
+ if (tau_i != 0.0)
+ {
+
+ /* 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);
+ }
+
+ scale = fabs(x) + fabs(y) + fabs(z);
+ if (scale != 0.0)
+ {
+ /* scale to prevent overflow or underflow */
+ x /= scale;
+ y /= scale;
+ z /= scale;
+ }
+ }
+
+ 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-16 10:36:55.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-18 17:09:04.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
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-17 20:45 Patrick.Alken
2006-05-18 16:24 ` Linas Vepstas
@ 2006-05-18 19:27 ` Brian Gough
1 sibling, 0 replies; 16+ messages in thread
From: Brian Gough @ 2006-05-18 19:27 UTC (permalink / raw)
To: Patrick.Alken; +Cc: gsl-discuss
Patrick.Alken@colorado.edu writes:
> 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.
Good work!
I would definitely like to add this (and the eigenvector algorithm and
generalised eigenproblem)
I will do some testing on it and get back to you with the results.
--
Brian Gough
(GSL Maintainer)
http://www.gnu.org/software/gsl/
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
@ 2006-05-18 16:43 Patrick.Alken
2006-05-19 15:13 ` C J Kenneth Tan -- OptimaNumerics
0 siblings, 1 reply; 16+ messages in thread
From: Patrick.Alken @ 2006-05-18 16:43 UTC (permalink / raw)
To: gsl-discuss
On Thu, May 18, 2006 at 11:21:33AM -0500, Linas Vepstas wrote:
> On Wed, May 17, 2006 at 02:45:10PM -0600, Patrick.Alken@colorado.edu wrote:
> >
> > I have recently needed to compute eigenvalues of unsymmetric
> > matrices,
>
> When I hit ths problem, someone recommended LAPACK, which
> I then used happily (routine DGEEV). This mailing list has
> seen occasional discusions about how to wrap/merge/whatever
> LAPACK routines with gsl (a technical problem is that LAPACK
> uses the fortran calling convention, making it ugly in C.)
LAPACK I believe uses the same QR algorithm which I have been
implementing, so if I can make it as robust as LAPACK's
implementation it should be a far better alternative than
linking to fortran code.
>
> > + f *= FLOAT_RADIX;
>
> Is this right? isn't there a "DOUBLE_RADIX" ?
I don't believe there is any standard define for the radix of
the floating point unit. The 'float' and 'double' radixes are
always the same of course - it is the base which the FPU uses
to represent float/double numbers. I used "FLOAT" to refer
to the floating point unit rather than float vs double.
As a side note, I noticed there is an underflow issue in the
patch I submitted yesterday which shows up for larger matrices
8x8+. I will try to fix this asap.
Patrick Alken
^ permalink raw reply [flat|nested] 16+ messages in thread
* Re: unsymmetric eigenvalue problem
2006-05-17 20:45 Patrick.Alken
@ 2006-05-18 16:24 ` Linas Vepstas
2006-05-22 21:37 ` Brian Gough
2006-05-18 19:27 ` Brian Gough
1 sibling, 1 reply; 16+ messages in thread
From: Linas Vepstas @ 2006-05-18 16:24 UTC (permalink / raw)
To: Patrick.Alken; +Cc: gsl-discuss
On Wed, May 17, 2006 at 02:45:10PM -0600, Patrick.Alken@colorado.edu wrote:
>
> I have recently needed to compute eigenvalues of unsymmetric
> matrices,
When I hit ths problem, someone recommended LAPACK, which
I then used happily (routine DGEEV). This mailing list has
seen occasional discusions about how to wrap/merge/whatever
LAPACK routines with gsl (a technical problem is that LAPACK
uses the fortran calling convention, making it ugly in C.)
> + f *= FLOAT_RADIX;
Is this right? isn't there a "DOUBLE_RADIX" ?
--linas
^ permalink raw reply [flat|nested] 16+ messages in thread
* unsymmetric eigenvalue problem
@ 2006-05-17 20:45 Patrick.Alken
2006-05-18 16:24 ` Linas Vepstas
2006-05-18 19:27 ` Brian Gough
0 siblings, 2 replies; 16+ messages in thread
From: Patrick.Alken @ 2006-05-17 20:45 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 1705 bytes --]
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
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: eig.c --]
[-- Type: text/x-csrc; name="eig.c", Size: 1795 bytes --]
/*
* 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;
}
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: gsl.unsymm.patch --]
[-- Type: text/x-patch; name="gsl.unsymm.patch", Size: 37650 bytes --]
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
^ permalink raw reply [flat|nested] 16+ messages in thread
end of thread, other threads:[~2006-05-26 10:43 UTC | newest]
Thread overview: 16+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2006-05-18 21:36 unsymmetric eigenvalue problem Patrick.Alken
-- strict thread matches above, loose matches on Subject: below --
2006-05-26 10:43 Patrick.Alken
2006-05-27 21:37 ` Brian Gough
2006-05-18 23:46 Patrick.Alken
2006-05-20 23:58 ` Patrick Alken
2006-05-21 5:17 ` James Bergstra
2006-05-21 10:13 ` Jochen Küpper
2006-05-22 22:26 ` Brian Gough
2006-05-23 18:35 ` Brian Gough
2006-05-18 16:43 Patrick.Alken
2006-05-19 15:13 ` C J Kenneth Tan -- OptimaNumerics
2006-05-19 16:57 ` Patrick Alken
2006-05-17 20:45 Patrick.Alken
2006-05-18 16:24 ` Linas Vepstas
2006-05-22 21:37 ` Brian Gough
2006-05-18 19:27 ` Brian Gough
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).