This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

unsymmetric eigenvalue problem


Hello all,

  I have recently needed to compute eigenvalues of unsymmetric
matrices, and so I wrote a code based in GSL to do it since
one does not currently exist. I followed the QR algorithm
outlined in Golub & Van Loan, chapter 7.

  I have attached a patch file to this email containing all
my code. Right now only the eigenvalue problem is solved,
the code will not compute eigenvectors. The patch is against
gsl-1.8 so just cd into a fresh gsl-1.8 and

patch -p0 < gsl.unsymm.patch

You will probably need to run automake again in the toplevel
dir to update the Makefile.am's. I have also attached a
test program eig.c to illustrate the use of the solver.

The patch creates 3 new source files:

linalg/hessenberg.c
eigen/balance.c
eigen/unsymm.c

and 1 header file
eigen/balance.h

hessenberg.c contains gsl_linalg_hessenberg() which reduces
a matrix to upper Hessenberg form.

balance.c contains balance_matrix() to do the standard matrix
balancing before finding eigenvalues.

unsymm.c contains the Francis QR algorithm to find the eigenvalues
of a general unsymmetric real matrix:
  gsl_eigen_unsymm_alloc()
  gsl_eigen_unsymm_free()
  gsl_eigen_unsymm() --- Solves A x = \lambda x

I am happy to contribute this to GSL if the developers want to
use it. Also, if there is interest, I would be willing to look
at coding up an eigenvector solver too for the unsymmetric
eigenvalue program - I only needed eigenvalues for my current
project so I haven't looked at the eigenvector problem.

Also if there is interest I'd be interested in looking at the
generalized eigenvalue problem since I feel GSL is incomplete
without an Ax = sBx solver.

Please test and give feedback.

Thanks,
Patrick Alken
/*
 * tester program for GSL unsymmetric eigenvalue solver
 *
 * To compile: gcc -o eig eig.c -lgsl -lgslcblas -lm
 */

#include <stdio.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_vector.h>

#define N 4

int
main(void)

{
#if N == 4
  double a_data[] = { 1.0, 2.0, 3.0, 4.0,
                      3.4, 2.2, 1.0, 9.7,
                      6.7, 5.0, 33.0, 8.9, 
                      78.7, 0.0, 4.5, 3.3 };
#endif
#if N == 2
  double a_data[] = { 3.0, -2.0,
                      4.0, -1.0 };
#endif
#if N == 3
  double a_data[] = { 3.0, -2.0, 0.0,
                      4.0, -1.0, 0.0,
                      1.0, 0.0, 1.0 };
#endif
  gsl_matrix_view m;
  gsl_eigen_unsymm_workspace *w;
  gsl_vector_complex *ev;
  int i,j;

  m = gsl_matrix_view_array(a_data, N, N);

  printf("original matrix\n");
  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
      printf("%e ", gsl_matrix_get(&m.matrix, i, j));
    printf("\n");
  }

  w = gsl_eigen_unsymm_alloc(N);
  ev = gsl_vector_complex_alloc(N);

  gsl_eigen_unsymm(&m.matrix, ev, w);

  printf("Schur decomposed matrix\n");
  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
      printf("%e ", gsl_matrix_get(&m.matrix, i, j));
    printf("\n");
  }

  for (i = 0; i < N; ++i)
    {
      gsl_complex a = gsl_vector_complex_get(ev, i);
      if (GSL_IMAG(a) == 0.0)
        {
          printf("ev[%d] = %f\n",
                 i,
                 GSL_REAL(a));
        }
      else
        {
          printf("ev[%d] = %f + %fi\n",
                 i,
                 GSL_REAL(a),
                 GSL_IMAG(a));
        }
    }

  gsl_eigen_unsymm_free(w);
  gsl_vector_complex_free(ev);

  return 0;
}
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
--- /home/palken/tmp2/gsl-1.8/eigen/balance.c	1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.c	2006-05-17 11:57:52.000000000 -0600
@@ -0,0 +1,117 @@
+/* eigen/balance.c
+ * 
+ * Copyright (C) 2006 Patrick Alken
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <stdlib.h>
+#include <math.h>
+
+#include "balance.h"
+
+/*
+balance_matrix()
+  Balance a given matrix by applying a diagonal matrix
+similarity transformation so that the rows and columns
+of the new matrix have norms which are the same order of
+magnitude. This is necessary for the unsymmetric eigenvalue
+problem since the calculation can become numerically unstable
+for unbalanced matrices.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
+
+and
+
+Numerical Recipes in C, Section 11.5
+*/
+
+void
+balance_matrix(gsl_matrix * A)
+{
+  double row_norm,
+         col_norm;
+  int not_converged;
+  const size_t N = A->size1;
+
+  not_converged = 1;
+
+  while (not_converged)
+    {
+      size_t i, j;
+      double g, f, s;
+
+      not_converged = 0;
+
+      for (i = 0; i < N; ++i)
+        {
+          row_norm = 0.0;
+          col_norm = 0.0;
+
+          for (j = 0; j < N; ++j)
+            {
+              if (j != i)
+                {
+                  col_norm += fabs(gsl_matrix_get(A, j, i));
+                  row_norm += fabs(gsl_matrix_get(A, i, j));
+                }
+            }
+
+          if ((col_norm == 0.0) || (row_norm == 0.0))
+            {
+              continue;
+            }
+
+          g = row_norm / FLOAT_RADIX;
+          f = 1.0;
+          s = col_norm + row_norm;
+
+          /*
+           * find the integer power of the machine radix which
+           * comes closest to balancing the matrix
+           */
+          while (col_norm < g)
+            {
+              f *= FLOAT_RADIX;
+              col_norm *= FLOAT_RADIX_SQ;
+            }
+
+          g = row_norm * FLOAT_RADIX;
+
+          while (col_norm > g)
+            {
+              f /= FLOAT_RADIX;
+              col_norm /= FLOAT_RADIX_SQ;
+            }
+
+          if ((row_norm + col_norm) < 0.95 * s * f)
+            {
+              not_converged = 1;
+
+              g = 1.0 / f;
+
+              /* apply similarity transformation */
+              for (j = 0; j < N; ++j)
+                {
+                  gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) * g);
+                }
+              for (j = 0; j < N; ++j)
+                {
+                  gsl_matrix_set(A, j, i, gsl_matrix_get(A, j, i) * f);
+                }
+            }
+        }
+    }
+}
diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
--- /home/palken/tmp2/gsl-1.8/eigen/balance.h	1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/balance.h	2006-05-17 10:49:00.000000000 -0600
@@ -0,0 +1,32 @@
+/* eigen/balance.h
+ * 
+ * Copyright (C) 2006 Patrick Alken
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __GSL_BALANCE_H__
+#define __GSL_BALANCE_H__
+
+#include <gsl/gsl_matrix.h>
+
+/* floating point radix */
+#define FLOAT_RADIX       2.0
+
+#define FLOAT_RADIX_SQ    (FLOAT_RADIX * FLOAT_RADIX)
+
+void balance_matrix(gsl_matrix * A);
+
+#endif /* __GSL_BALANCE_H__ */
diff -urN /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
--- /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h	2005-06-26 07:25:34.000000000 -0600
+++ ./eigen/gsl_eigen.h	2006-05-17 12:02:26.000000000 -0600
@@ -59,6 +59,16 @@
 
 typedef struct {
   size_t size;
+  gsl_vector *v; /* temporary vector */
+} gsl_eigen_unsymm_workspace;
+
+gsl_eigen_unsymm_workspace * gsl_eigen_unsymm_alloc (const size_t n);
+void gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w);
+int gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+                      gsl_eigen_unsymm_workspace * w);
+
+typedef struct {
+  size_t size;
   double * d;
   double * sd;
   double * tau;
diff -urN /home/palken/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
--- /home/palken/tmp2/gsl-1.8/eigen/Makefile.am	2004-09-11 07:45:45.000000000 -0600
+++ ./eigen/Makefile.am	2006-05-17 12:00:33.000000000 -0600
@@ -3,7 +3,7 @@
 check_PROGRAMS = test
 
 pkginclude_HEADERS = gsl_eigen.h
-libgsleigen_la_SOURCES =  jacobi.c symm.c symmv.c herm.c hermv.c sort.c
+libgsleigen_la_SOURCES =  jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c balance.c
 INCLUDES= -I$(top_builddir)
 
 noinst_HEADERS =  qrstep.c 
diff -urN /home/palken/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
--- /home/palken/tmp2/gsl-1.8/eigen/unsymm.c	1969-12-31 17:00:00.000000000 -0700
+++ ./eigen/unsymm.c	2006-05-17 11:57:57.000000000 -0600
@@ -0,0 +1,378 @@
+/* eigen/unsymm.c
+ * 
+ * Copyright (C) 2006 Patrick Alken
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <stdlib.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_vector_complex.h>
+#include <gsl/gsl_matrix.h>
+
+#include "balance.h"
+
+/* maximum iterations before failure is reported */
+#define UNSYMM_MAX_ITERATIONS     100
+
+/*
+ * This module computes the eigenvalues of a real unsymmetric
+ * matrix, using the QR decomposition.
+ *
+ * See Golub & Van Loan, "Matrix Computations" (3rd ed),
+ * algorithm 7.5.2
+ */
+
+inline static void zero_subdiag_small_elements(gsl_matrix * A);
+static inline int francis_qrstep(gsl_matrix * H,
+                                 gsl_eigen_unsymm_workspace * w);
+static void get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+                                gsl_complex * e2);
+
+/*
+gsl_eigen_unsymm_alloc()
+
+Allocate a workspace for solving the unsymmetric
+eigenvalue/eigenvector problem
+
+Inputs: n - size of matrix
+
+Return: pointer to workspace
+*/
+
+gsl_eigen_unsymm_workspace *
+gsl_eigen_unsymm_alloc(const size_t n)
+{
+  gsl_eigen_unsymm_workspace *w;
+
+  if (n == 0)
+    {
+      GSL_ERROR_NULL ("matrix dimension must be positive integer",
+                      GSL_EINVAL);
+    }
+
+  w = ((gsl_eigen_unsymm_workspace *)
+       malloc (sizeof (gsl_eigen_unsymm_workspace)));
+
+  if (w == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
+    }
+
+  w->size = n;
+  w->v = gsl_vector_alloc(3);
+
+  return (w);
+}
+
+void
+gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
+{
+  gsl_vector_free(w->v);
+
+  free(w);
+}
+
+/*
+gsl_eigen_unsymm()
+
+Solve the unsymmetric eigenvalue problem
+
+A x = \lambda x
+
+for the eigenvalues \lambda using algorithm 7.5.2 of
+Golub & Van Loan, "Matrix Computations" (3rd ed)
+
+Inputs: A    - matrix
+        eval - where to store eigenvalues
+        w    - workspace
+
+Notes: the matrix A is destroyed by this routine
+*/
+
+int
+gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
+                  gsl_eigen_unsymm_workspace * w)
+{
+  /* check matrix and vector sizes */
+
+  if (A->size1 != A->size2)
+    {
+      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
+    }
+  else if (eval->size != A->size1)
+    {
+      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
+    }
+  else
+    {
+      size_t N;
+      size_t q;
+      gsl_complex lambda1, lambda2; /* eigenvalues */
+      gsl_matrix_view m;
+      size_t nev;
+      size_t nit;
+
+      N = A->size1;
+
+      /* special cases */
+      if (N == 1)
+      {
+        GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(A, 0, 0), 0.0);
+        gsl_vector_complex_set(eval, 0, lambda1);
+        return GSL_SUCCESS;
+      }
+
+      if (N == 2)
+      {
+        /*
+         * The 2x2 case is special since the matrix is already
+         * in upper quasi-triangular form so no Schur decomposition
+         * is necessary
+         */
+        get_2b2_eigenvalues(A, &lambda1, &lambda2);
+        gsl_vector_complex_set(eval, 0, lambda1);
+        gsl_vector_complex_set(eval, 1, lambda2);
+        return GSL_SUCCESS;
+      }
+
+      nev = 0;
+
+      /* balance the matrix */
+      balance_matrix(A);
+
+      /* compute the Hessenberg reduction of A */
+      gsl_linalg_hessenberg(A);
+
+      /* set small elements on the subdiagonal to 0 */
+      zero_subdiag_small_elements(A);
+
+      m = gsl_matrix_submatrix(A, 0, 0, N, N);
+
+      nit = 0;
+      while ((N > 2) && (nit++ < UNSYMM_MAX_ITERATIONS))
+        {
+          francis_qrstep(&m.matrix, w);
+          zero_subdiag_small_elements(&m.matrix);
+
+          q = N - 1;
+          while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
+            {
+              --q;
+            }
+
+          if (q == (N - 1))
+            {
+              /*
+               * the (N, N - 1) element of the matrix is 0 -
+               * m_{NN} is a real eigenvalue
+               */
+              GSL_SET_COMPLEX(&lambda1,
+                              gsl_matrix_get(&m.matrix, q, q), 0.0);
+              gsl_vector_complex_set(eval, nev++, lambda1);
+
+              --N;
+              m = gsl_matrix_submatrix(A, 0, 0, N, N);
+            }
+          else if (q == (N - 2))
+            {
+              /*
+               * The bottom right 2x2 block of m is an eigenvalue
+               * system
+               */
+
+              m = gsl_matrix_submatrix(A, q, q, 2, 2);
+              get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+              gsl_vector_complex_set(eval, nev++, lambda1);
+              gsl_vector_complex_set(eval, nev++, lambda2);
+
+              N -= 2;
+              m = gsl_matrix_submatrix(A, 0, 0, N, N);
+            }
+        }
+
+      if (N == 1)
+        {
+          GSL_SET_COMPLEX(&lambda1,
+                          gsl_matrix_get(&m.matrix, 0, 0), 0.0);
+          gsl_vector_complex_set(eval, nev++, lambda1);
+        }
+      else if (N == 2)
+        {
+          get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
+          gsl_vector_complex_set(eval, nev++, lambda1);
+          gsl_vector_complex_set(eval, nev++, lambda2);
+        }
+
+      return GSL_SUCCESS;
+    }
+} /* gsl_eigen_unsymm() */
+
+/********************************************
+ *           INTERNAL ROUTINES              *
+ ********************************************/
+
+/*
+zero_subdiag_small_elements()
+  Sets to zero all elements on the subdiaganal of a matrix A
+which satisfy
+
+|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
+*/
+
+inline static void
+zero_subdiag_small_elements(gsl_matrix * A)
+{
+  const size_t N = A->size1;
+  size_t i;
+  double dpel = gsl_matrix_get(A, 0, 0);
+
+  for (i = 1; i < N; ++i)
+    {
+      double sel = gsl_matrix_get(A, i, i - 1);
+      double del = gsl_matrix_get(A, i, i);
+
+      if (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel)))
+        {
+          gsl_matrix_set(A, i, i - 1, 0.0);
+        }
+
+      dpel = del;
+    }
+}
+
+/*
+francis_qrstep()
+  Perform a Francis QR step.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed),
+algorithm 7.5.1
+
+Inputs: H - unreduced upper Hessenberg matrix
+        w - workspace
+*/
+
+static inline int
+francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w)
+{
+  const size_t N = H->size1;
+  double s, t, x, y, z;
+  size_t i;
+  gsl_vector_view v;
+  gsl_matrix_view m;
+  double tau_i;
+  size_t q, r;
+
+  /* s = a_1 + a_2 = H_{mm} + H_{nn} */
+  s = gsl_matrix_get(H, N - 2, N - 2) + gsl_matrix_get(H, N - 1, N - 1);
+
+  /* t = a_1 * a_2 = H_{mm} * H_{nn} - H_{mn} * H_{nm} */
+  t = gsl_matrix_get(H, N - 2, N - 2) * gsl_matrix_get(H, N - 1, N - 1) -
+      gsl_matrix_get(H, N - 2, N - 1) * gsl_matrix_get(H, N - 1, N - 2);
+  x = gsl_matrix_get(H, 0, 0) * gsl_matrix_get(H, 0, 0) +
+      gsl_matrix_get(H, 0, 1) * gsl_matrix_get(H, 1, 0) -
+      s*gsl_matrix_get(H, 0, 0) + t;
+  y = gsl_matrix_get(H, 1, 0) *
+      (gsl_matrix_get(H, 0, 0) + gsl_matrix_get(H, 1, 1) - s);
+  z = gsl_matrix_get(H, 1, 0) * gsl_matrix_get(H, 2, 1);
+
+  v = gsl_vector_subvector(w->v, 0, 3);
+
+  for (i = 0; i < N - 2; ++i)
+    {
+      gsl_vector_set(&v.vector, 0, x);
+      gsl_vector_set(&v.vector, 1, y);
+      gsl_vector_set(&v.vector, 2, z);
+      tau_i = gsl_linalg_householder_transform(&v.vector);
+
+      /* q = max(1, i - 1) */
+      q = (1 > ((int)i - 1)) ? 0 : (i - 1);
+
+      /* apply left householder matrix (I - tau_i v v') to H */
+      m = gsl_matrix_submatrix(H, i, q, 3, N - q);
+      gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+      /* r = min(i + 3, N - 1) */
+      r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
+
+      /* apply right householder matrix (I - tau_i v v') to H */
+      m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
+      gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+      x = gsl_matrix_get(H, i + 1, i);
+      y = gsl_matrix_get(H, i + 2, i);
+      if (i < (N - 3))
+        {
+          z = gsl_matrix_get(H, i + 3, i);
+        }
+    }
+
+  v = gsl_vector_subvector(w->v, 0, 2);
+  gsl_vector_set(&v.vector, 0, x);
+  gsl_vector_set(&v.vector, 1, y);
+  tau_i = gsl_linalg_householder_transform(&v.vector);
+
+  m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
+  gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+  m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
+  gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+
+  return GSL_SUCCESS;
+}
+
+/*
+get_2b2_eigenvalues()
+  Compute the eigenvalues of a 2x2 real matrix
+
+Inputs: A  - 2x2 matrix
+        e1 - where to store eigenvalue 1
+        e2 - where to store eigenvalue 2
+*/
+
+static void
+get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
+                    gsl_complex * e2)
+{
+  double discr;      /* discriminant of characteristic poly */
+  double a, b, c, d; /* matrix values */
+
+  a = gsl_matrix_get(A, 0, 0);
+  b = gsl_matrix_get(A, 0, 1);
+  c = gsl_matrix_get(A, 1, 0);
+  d = gsl_matrix_get(A, 1, 1);
+
+  discr = (a + d)*(a + d) - 4.0*(a*d - b*c);
+  if (discr < 0.0)
+    {
+      GSL_SET_REAL(e1, 0.5*(a + d));
+      GSL_SET_REAL(e2, 0.5*(a + d));
+
+      GSL_SET_IMAG(e1, 0.5*sqrt(-discr));
+      GSL_SET_IMAG(e2, -0.5*sqrt(-discr));
+    }
+  else
+    {
+      GSL_SET_REAL(e1, 0.5*(a + d + sqrt(discr)));
+      GSL_SET_REAL(e2, 0.5*(a + d - sqrt(discr)));
+
+      GSL_SET_IMAG(e1, 0.0);
+      GSL_SET_IMAG(e2, 0.0);
+    }
+}
diff -urN /home/palken/tmp2/gsl-1.8/gsl.unsymm.patch ./gsl.unsymm.patch
--- /home/palken/tmp2/gsl-1.8/gsl.unsymm.patch	1969-12-31 17:00:00.000000000 -0700
+++ ./gsl.unsymm.patch	2006-05-17 12:07:19.000000000 -0600
@@ -0,0 +1,571 @@
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.c ./eigen/balance.c
+--- /home/palken/tmp2/gsl-1.8/eigen/balance.c	1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/balance.c	2006-05-17 11:57:52.000000000 -0600
+@@ -0,0 +1,117 @@
++/* eigen/balance.c
++ * 
++ * Copyright (C) 2006 Patrick Alken
++ * 
++ * This program is free software; you can redistribute it and/or modify
++ * it under the terms of the GNU General Public License as published by
++ * the Free Software Foundation; either version 2 of the License, or (at
++ * your option) any later version.
++ * 
++ * This program is distributed in the hope that it will be useful, but
++ * WITHOUT ANY WARRANTY; without even the implied warranty of
++ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
++ * General Public License for more details.
++ * 
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#include <stdlib.h>
++#include <math.h>
++
++#include "balance.h"
++
++/*
++balance_matrix()
++  Balance a given matrix by applying a diagonal matrix
++similarity transformation so that the rows and columns
++of the new matrix have norms which are the same order of
++magnitude. This is necessary for the unsymmetric eigenvalue
++problem since the calculation can become numerically unstable
++for unbalanced matrices.
++
++See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
++
++and
++
++Numerical Recipes in C, Section 11.5
++*/
++
++void
++balance_matrix(gsl_matrix * A)
++{
++  double row_norm,
++         col_norm;
++  int not_converged;
++  const size_t N = A->size1;
++
++  not_converged = 1;
++
++  while (not_converged)
++    {
++      size_t i, j;
++      double g, f, s;
++
++      not_converged = 0;
++
++      for (i = 0; i < N; ++i)
++        {
++          row_norm = 0.0;
++          col_norm = 0.0;
++
++          for (j = 0; j < N; ++j)
++            {
++              if (j != i)
++                {
++                  col_norm += fabs(gsl_matrix_get(A, j, i));
++                  row_norm += fabs(gsl_matrix_get(A, i, j));
++                }
++            }
++
++          if ((col_norm == 0.0) || (row_norm == 0.0))
++            {
++              continue;
++            }
++
++          g = row_norm / FLOAT_RADIX;
++          f = 1.0;
++          s = col_norm + row_norm;
++
++          /*
++           * find the integer power of the machine radix which
++           * comes closest to balancing the matrix
++           */
++          while (col_norm < g)
++            {
++              f *= FLOAT_RADIX;
++              col_norm *= FLOAT_RADIX_SQ;
++            }
++
++          g = row_norm * FLOAT_RADIX;
++
++          while (col_norm > g)
++            {
++              f /= FLOAT_RADIX;
++              col_norm /= FLOAT_RADIX_SQ;
++            }
++
++          if ((row_norm + col_norm) < 0.95 * s * f)
++            {
++              not_converged = 1;
++
++              g = 1.0 / f;
++
++              /* apply similarity transformation */
++              for (j = 0; j < N; ++j)
++                {
++                  gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) * g);
++                }
++              for (j = 0; j < N; ++j)
++                {
++                  gsl_matrix_set(A, j, i, gsl_matrix_get(A, j, i) * f);
++                }
++            }
++        }
++    }
++}
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/balance.h ./eigen/balance.h
+--- /home/palken/tmp2/gsl-1.8/eigen/balance.h	1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/balance.h	2006-05-17 10:49:00.000000000 -0600
+@@ -0,0 +1,32 @@
++/* eigen/balance.h
++ * 
++ * Copyright (C) 2006 Patrick Alken
++ * 
++ * This program is free software; you can redistribute it and/or modify
++ * it under the terms of the GNU General Public License as published by
++ * the Free Software Foundation; either version 2 of the License, or (at
++ * your option) any later version.
++ * 
++ * This program is distributed in the hope that it will be useful, but
++ * WITHOUT ANY WARRANTY; without even the implied warranty of
++ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
++ * General Public License for more details.
++ * 
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#ifndef __GSL_BALANCE_H__
++#define __GSL_BALANCE_H__
++
++#include <gsl/gsl_matrix.h>
++
++/* floating point radix */
++#define FLOAT_RADIX       2.0
++
++#define FLOAT_RADIX_SQ    (FLOAT_RADIX * FLOAT_RADIX)
++
++void balance_matrix(gsl_matrix * A);
++
++#endif /* __GSL_BALANCE_H__ */
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h ./eigen/gsl_eigen.h
+--- /home/palken/tmp2/gsl-1.8/eigen/gsl_eigen.h	2005-06-26 07:25:34.000000000 -0600
++++ ./eigen/gsl_eigen.h	2006-05-17 12:02:26.000000000 -0600
+@@ -59,6 +59,16 @@
+ 
+ typedef struct {
+   size_t size;
++  gsl_vector *v; /* temporary vector */
++} gsl_eigen_unsymm_workspace;
++
++gsl_eigen_unsymm_workspace * gsl_eigen_unsymm_alloc (const size_t n);
++void gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w);
++int gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
++                      gsl_eigen_unsymm_workspace * w);
++
++typedef struct {
++  size_t size;
+   double * d;
+   double * sd;
+   double * tau;
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/Makefile.am ./eigen/Makefile.am
+--- /home/palken/tmp2/gsl-1.8/eigen/Makefile.am	2004-09-11 07:45:45.000000000 -0600
++++ ./eigen/Makefile.am	2006-05-17 12:00:33.000000000 -0600
+@@ -3,7 +3,7 @@
+ check_PROGRAMS = test
+ 
+ pkginclude_HEADERS = gsl_eigen.h
+-libgsleigen_la_SOURCES =  jacobi.c symm.c symmv.c herm.c hermv.c sort.c
++libgsleigen_la_SOURCES =  jacobi.c symm.c symmv.c unsymm.c herm.c hermv.c sort.c balance.c
+ INCLUDES= -I$(top_builddir)
+ 
+ noinst_HEADERS =  qrstep.c 
+diff -urN /home/palken/tmp2/gsl-1.8/eigen/unsymm.c ./eigen/unsymm.c
+--- /home/palken/tmp2/gsl-1.8/eigen/unsymm.c	1969-12-31 17:00:00.000000000 -0700
++++ ./eigen/unsymm.c	2006-05-17 11:57:57.000000000 -0600
+@@ -0,0 +1,378 @@
++/* eigen/unsymm.c
++ * 
++ * Copyright (C) 2006 Patrick Alken
++ * 
++ * This program is free software; you can redistribute it and/or modify
++ * it under the terms of the GNU General Public License as published by
++ * the Free Software Foundation; either version 2 of the License, or (at
++ * your option) any later version.
++ * 
++ * This program is distributed in the hope that it will be useful, but
++ * WITHOUT ANY WARRANTY; without even the implied warranty of
++ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
++ * General Public License for more details.
++ * 
++ * You should have received a copy of the GNU General Public License
++ * along with this program; if not, write to the Free Software
++ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
++ */
++
++#include <stdlib.h>
++#include <gsl/gsl_eigen.h>
++#include <gsl/gsl_linalg.h>
++#include <gsl/gsl_math.h>
++#include <gsl/gsl_blas.h>
++#include <gsl/gsl_vector.h>
++#include <gsl/gsl_vector_complex.h>
++#include <gsl/gsl_matrix.h>
++
++#include "balance.h"
++
++/* maximum iterations before failure is reported */
++#define UNSYMM_MAX_ITERATIONS     100
++
++/*
++ * This module computes the eigenvalues of a real unsymmetric
++ * matrix, using the QR decomposition.
++ *
++ * See Golub & Van Loan, "Matrix Computations" (3rd ed),
++ * algorithm 7.5.2
++ */
++
++inline static void zero_subdiag_small_elements(gsl_matrix * A);
++static inline int francis_qrstep(gsl_matrix * H,
++                                 gsl_eigen_unsymm_workspace * w);
++static void get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
++                                gsl_complex * e2);
++
++/*
++gsl_eigen_unsymm_alloc()
++
++Allocate a workspace for solving the unsymmetric
++eigenvalue/eigenvector problem
++
++Inputs: n - size of matrix
++
++Return: pointer to workspace
++*/
++
++gsl_eigen_unsymm_workspace *
++gsl_eigen_unsymm_alloc(const size_t n)
++{
++  gsl_eigen_unsymm_workspace *w;
++
++  if (n == 0)
++    {
++      GSL_ERROR_NULL ("matrix dimension must be positive integer",
++                      GSL_EINVAL);
++    }
++
++  w = ((gsl_eigen_unsymm_workspace *)
++       malloc (sizeof (gsl_eigen_unsymm_workspace)));
++
++  if (w == 0)
++    {
++      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
++    }
++
++  w->size = n;
++  w->v = gsl_vector_alloc(3);
++
++  return (w);
++}
++
++void
++gsl_eigen_unsymm_free (gsl_eigen_unsymm_workspace * w)
++{
++  gsl_vector_free(w->v);
++
++  free(w);
++}
++
++/*
++gsl_eigen_unsymm()
++
++Solve the unsymmetric eigenvalue problem
++
++A x = \lambda x
++
++for the eigenvalues \lambda using algorithm 7.5.2 of
++Golub & Van Loan, "Matrix Computations" (3rd ed)
++
++Inputs: A    - matrix
++        eval - where to store eigenvalues
++        w    - workspace
++
++Notes: the matrix A is destroyed by this routine
++*/
++
++int
++gsl_eigen_unsymm (gsl_matrix * A, gsl_vector_complex * eval,
++                  gsl_eigen_unsymm_workspace * w)
++{
++  /* check matrix and vector sizes */
++
++  if (A->size1 != A->size2)
++    {
++      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
++    }
++  else if (eval->size != A->size1)
++    {
++      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
++    }
++  else
++    {
++      size_t N;
++      size_t q;
++      gsl_complex lambda1, lambda2; /* eigenvalues */
++      gsl_matrix_view m;
++      size_t nev;
++      size_t nit;
++
++      N = A->size1;
++
++      /* special cases */
++      if (N == 1)
++      {
++        GSL_SET_COMPLEX(&lambda1, gsl_matrix_get(A, 0, 0), 0.0);
++        gsl_vector_complex_set(eval, 0, lambda1);
++        return GSL_SUCCESS;
++      }
++
++      if (N == 2)
++      {
++        /*
++         * The 2x2 case is special since the matrix is already
++         * in upper quasi-triangular form so no Schur decomposition
++         * is necessary
++         */
++        get_2b2_eigenvalues(A, &lambda1, &lambda2);
++        gsl_vector_complex_set(eval, 0, lambda1);
++        gsl_vector_complex_set(eval, 1, lambda2);
++        return GSL_SUCCESS;
++      }
++
++      nev = 0;
++
++      /* balance the matrix */
++      balance_matrix(A);
++
++      /* compute the Hessenberg reduction of A */
++      gsl_linalg_hessenberg(A);
++
++      /* set small elements on the subdiagonal to 0 */
++      zero_subdiag_small_elements(A);
++
++      m = gsl_matrix_submatrix(A, 0, 0, N, N);
++
++      nit = 0;
++      while ((N > 2) && (nit++ < UNSYMM_MAX_ITERATIONS))
++        {
++          francis_qrstep(&m.matrix, w);
++          zero_subdiag_small_elements(&m.matrix);
++
++          q = N - 1;
++          while ((q > 0) && (gsl_matrix_get(A, q, q - 1) != 0.0))
++            {
++              --q;
++            }
++
++          if (q == (N - 1))
++            {
++              /*
++               * the (N, N - 1) element of the matrix is 0 -
++               * m_{NN} is a real eigenvalue
++               */
++              GSL_SET_COMPLEX(&lambda1,
++                              gsl_matrix_get(&m.matrix, q, q), 0.0);
++              gsl_vector_complex_set(eval, nev++, lambda1);
++
++              --N;
++              m = gsl_matrix_submatrix(A, 0, 0, N, N);
++            }
++          else if (q == (N - 2))
++            {
++              /*
++               * The bottom right 2x2 block of m is an eigenvalue
++               * system
++               */
++
++              m = gsl_matrix_submatrix(A, q, q, 2, 2);
++              get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
++              gsl_vector_complex_set(eval, nev++, lambda1);
++              gsl_vector_complex_set(eval, nev++, lambda2);
++
++              N -= 2;
++              m = gsl_matrix_submatrix(A, 0, 0, N, N);
++            }
++        }
++
++      if (N == 1)
++        {
++          GSL_SET_COMPLEX(&lambda1,
++                          gsl_matrix_get(&m.matrix, 0, 0), 0.0);
++          gsl_vector_complex_set(eval, nev++, lambda1);
++        }
++      else if (N == 2)
++        {
++          get_2b2_eigenvalues(&m.matrix, &lambda1, &lambda2);
++          gsl_vector_complex_set(eval, nev++, lambda1);
++          gsl_vector_complex_set(eval, nev++, lambda2);
++        }
++
++      return GSL_SUCCESS;
++    }
++} /* gsl_eigen_unsymm() */
++
++/********************************************
++ *           INTERNAL ROUTINES              *
++ ********************************************/
++
++/*
++zero_subdiag_small_elements()
++  Sets to zero all elements on the subdiaganal of a matrix A
++which satisfy
++
++|A_{i,i-1}| <= eps * (|A_{i,i}| + |A_{i-1,i-1}|)
++*/
++
++inline static void
++zero_subdiag_small_elements(gsl_matrix * A)
++{
++  const size_t N = A->size1;
++  size_t i;
++  double dpel = gsl_matrix_get(A, 0, 0);
++
++  for (i = 1; i < N; ++i)
++    {
++      double sel = gsl_matrix_get(A, i, i - 1);
++      double del = gsl_matrix_get(A, i, i);
++
++      if (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel)))
++        {
++          gsl_matrix_set(A, i, i - 1, 0.0);
++        }
++
++      dpel = del;
++    }
++}
++
++/*
++francis_qrstep()
++  Perform a Francis QR step.
++
++See Golub & Van Loan, "Matrix Computations" (3rd ed),
++algorithm 7.5.1
++
++Inputs: H - unreduced upper Hessenberg matrix
++        w - workspace
++*/
++
++static inline int
++francis_qrstep(gsl_matrix * H, gsl_eigen_unsymm_workspace * w)
++{
++  const size_t N = H->size1;
++  double s, t, x, y, z;
++  size_t i;
++  gsl_vector_view v;
++  gsl_matrix_view m;
++  double tau_i;
++  size_t q, r;
++
++  /* s = a_1 + a_2 = H_{mm} + H_{nn} */
++  s = gsl_matrix_get(H, N - 2, N - 2) + gsl_matrix_get(H, N - 1, N - 1);
++
++  /* t = a_1 * a_2 = H_{mm} * H_{nn} - H_{mn} * H_{nm} */
++  t = gsl_matrix_get(H, N - 2, N - 2) * gsl_matrix_get(H, N - 1, N - 1) -
++      gsl_matrix_get(H, N - 2, N - 1) * gsl_matrix_get(H, N - 1, N - 2);
++  x = gsl_matrix_get(H, 0, 0) * gsl_matrix_get(H, 0, 0) +
++      gsl_matrix_get(H, 0, 1) * gsl_matrix_get(H, 1, 0) -
++      s*gsl_matrix_get(H, 0, 0) + t;
++  y = gsl_matrix_get(H, 1, 0) *
++      (gsl_matrix_get(H, 0, 0) + gsl_matrix_get(H, 1, 1) - s);
++  z = gsl_matrix_get(H, 1, 0) * gsl_matrix_get(H, 2, 1);
++
++  v = gsl_vector_subvector(w->v, 0, 3);
++
++  for (i = 0; i < N - 2; ++i)
++    {
++      gsl_vector_set(&v.vector, 0, x);
++      gsl_vector_set(&v.vector, 1, y);
++      gsl_vector_set(&v.vector, 2, z);
++      tau_i = gsl_linalg_householder_transform(&v.vector);
++
++      /* q = max(1, i - 1) */
++      q = (1 > ((int)i - 1)) ? 0 : (i - 1);
++
++      /* apply left householder matrix (I - tau_i v v') to H */
++      m = gsl_matrix_submatrix(H, i, q, 3, N - q);
++      gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
++
++      /* r = min(i + 3, N - 1) */
++      r = ((i + 3) < (N - 1)) ? (i + 3) : (N - 1);
++
++      /* apply right householder matrix (I - tau_i v v') to H */
++      m = gsl_matrix_submatrix(H, 0, i, r + 1, 3);
++      gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
++
++      x = gsl_matrix_get(H, i + 1, i);
++      y = gsl_matrix_get(H, i + 2, i);
++      if (i < (N - 3))
++        {
++          z = gsl_matrix_get(H, i + 3, i);
++        }
++    }
++
++  v = gsl_vector_subvector(w->v, 0, 2);
++  gsl_vector_set(&v.vector, 0, x);
++  gsl_vector_set(&v.vector, 1, y);
++  tau_i = gsl_linalg_householder_transform(&v.vector);
++
++  m = gsl_matrix_submatrix(H, N - 2, N - 3, 2, 3);
++  gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
++
++  m = gsl_matrix_submatrix(H, 0, N - 2, N, 2);
++  gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
++
++  return GSL_SUCCESS;
++}
++
++/*
++get_2b2_eigenvalues()
++  Compute the eigenvalues of a 2x2 real matrix
++
++Inputs: A  - 2x2 matrix
++        e1 - where to store eigenvalue 1
++        e2 - where to store eigenvalue 2
++*/
++
++static void
++get_2b2_eigenvalues(gsl_matrix * A, gsl_complex * e1,
++                    gsl_complex * e2)
++{
++  double discr;      /* discriminant of characteristic poly */
++  double a, b, c, d; /* matrix values */
++
++  a = gsl_matrix_get(A, 0, 0);
++  b = gsl_matrix_get(A, 0, 1);
++  c = gsl_matrix_get(A, 1, 0);
++  d = gsl_matrix_get(A, 1, 1);
++
++  discr = (a + d)*(a + d) - 4.0*(a*d - b*c);
++  if (discr < 0.0)
++    {
++      GSL_SET_REAL(e1, 0.5*(a + d));
++      GSL_SET_REAL(e2, 0.5*(a + d));
++
++      GSL_SET_IMAG(e1, 0.5*sqrt(-discr));
++      GSL_SET_IMAG(e2, -0.5*sqrt(-discr));
++    }
++  else
++    {
++      GSL_SET_REAL(e1, 0.5*(a + d + sqrt(discr)));
++      GSL_SET_REAL(e2, 0.5*(a + d - sqrt(discr)));
++
++      GSL_SET_IMAG(e1, 0.0);
++      GSL_SET_IMAG(e2, 0.0);
++    }
++}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h ./linalg/gsl_linalg.h
--- /home/palken/tmp2/gsl-1.8/linalg/gsl_linalg.h	2005-11-09 14:13:14.000000000 -0700
+++ ./linalg/gsl_linalg.h	2006-05-17 12:01:43.000000000 -0600
@@ -113,6 +113,10 @@
                                        const gsl_vector_complex * v, 
                                        gsl_vector_complex * w);
 
+/* Hessenberg reduction */
+
+int gsl_linalg_hessenberg (gsl_matrix * A);
+
 /* Singular Value Decomposition
 
  * exceptions: 
diff -urN /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c ./linalg/hessenberg.c
--- /home/palken/tmp2/gsl-1.8/linalg/hessenberg.c	1969-12-31 17:00:00.000000000 -0700
+++ ./linalg/hessenberg.c	2006-05-16 15:11:12.000000000 -0600
@@ -0,0 +1,76 @@
+/* linalg/hessenberg.c
+ * 
+ * Copyright (C) 2006 Patrick Alken
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+/*
+gsl_linalg_hessenberg()
+  Compute the Householder reduction to Hessenberg form of a
+square matrix A.
+
+See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm
+7.4.2
+*/
+
+int
+gsl_linalg_hessenberg(gsl_matrix * A)
+{
+  if (A->size1 != A->size2)
+    {
+      GSL_ERROR ("Hessenberg reduction requires square matrix",
+                 GSL_ENOTSQR);
+    }
+  else
+    {
+      const size_t N = A->size1;
+      size_t i, j;
+      gsl_vector_view v;
+      gsl_vector *v_copy;
+      gsl_matrix_view m;
+      double tau_i; /* beta in algorithm 7.4.2 */
+
+      v_copy = gsl_vector_alloc(N);
+
+      for (i = 0; i < N - 2; ++i)
+        {
+          /* make a copy of A(i + 1:n, i) */
+
+          v = gsl_matrix_column(A, i);
+          gsl_vector_memcpy(v_copy, &v.vector);
+          v = gsl_vector_subvector(v_copy, i + 1, N - (i + 1));
+
+          /* compute householder transformation of A(i+1:n,i) */
+          tau_i = gsl_linalg_householder_transform(&v.vector);
+
+          /* apply left householder matrix (I - tau_i v v') to A */
+          m = gsl_matrix_submatrix(A, i + 1, i, N - (i + 1), N - i);
+          gsl_linalg_householder_hm(tau_i, &v.vector, &m.matrix);
+
+          /* apply right householder matrix (I - tau_i v v') to A */
+          m = gsl_matrix_submatrix(A, 0, i + 1, N, N - (i + 1));
+          gsl_linalg_householder_mh(tau_i, &v.vector, &m.matrix);
+        }
+
+      gsl_vector_free(v_copy);
+
+      return GSL_SUCCESS;
+    }
+}
diff -urN /home/palken/tmp2/gsl-1.8/linalg/Makefile.am ./linalg/Makefile.am
--- /home/palken/tmp2/gsl-1.8/linalg/Makefile.am	2004-09-13 12:17:04.000000000 -0600
+++ ./linalg/Makefile.am	2006-05-17 12:01:11.000000000 -0600
@@ -4,7 +4,7 @@
 
 INCLUDES= -I$(top_builddir)
 
-libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
+libgsllinalg_la_SOURCES = multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c qr.c qrpt.c lq.c ptlq.c svd.c householder.c householdercomplex.c hessenberg.c cholesky.c symmtd.c hermtd.c bidiag.c balance.c
 
 noinst_HEADERS =  givens.c apply_givens.c svdstep.c tridiag.h 
 

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]