unsymmetric eigenvalue update

unsymmetric eigenvalue update

[Patrick Alken]

[-- Attachment #1: Type: text/plain, Size: 3238 bytes --]

Hi all,

  After reviewing the LAPACK and EISPACK eigenvalue routines,
I've changed my convergence criteria to declare failure after
30*N iterations. My original number of 30 was based on the
Numerical Recipes algorithm, which I now see was copied
from the EISPACK algorithm but with a typo of 30 instead of
30*N.

  I have run my code through Brian Gough's test program
on 12x12 and 200x200 random matrices, as well as some integer
matrix cases and find no significant disagreement with LAPACK.

  For those who expressed concern earlier about speed issues,
I'd like to share some things I found while studying this
problem:

1) In 1962 the Francis double-shift QR algorithm was published
   which typically converges in about 8n^3 flops if you only
   want eigenvalues. This is the algorithm which was implemented
   in EISPACK (hqr.f) and is also implemented in LAPACK as
   DLAHQR. Also this is the algorithm discussed in Golub &
   Van Loan and is presented in Numerical Recipes (which is just
   the EISPACK hqr routine converted to C).

2) In 1989, Bai and Demmel published their multi-shift QR
   algorithm and gave some numerical tests showing that it
   is generally faster than hqr from EISPACK. However they
   say its too complicated to derive an approximate number of
   flops required for the general problem. Furthermore,
   the multi-shift algorithm needs to call the double-shift
   algorithm to calculate the shifts. So it is impossible
   to dispense with the double-shift algorithm. This algorithm
   is implemented in LAPACK as DGEES (DGEES calls DLAHQR to
   get the shifts)

3) In a further significant development, in 2001 a multishift
   QR algorithm was proposed in the papers

   K. Braman et al, "The Multishift QR Algorithm: Part I:
   Maintaining Well Focused Shifts"

   and

   K. Braman et al, "The Multishift QR Algorithm: Part II:
   Aggressive Early Deflation"

   which introduces a way to deflate eigenvalues much more
   quickly than the standard test used in the previous two
   methods. The numerical results in these papers are very
   impressive and show convergence with far less flops than
   the LAPACK algorithm. Subsequent papers have extended this
   method to the QZ algorithm of the generalized eigenvalue
   problem, and as far as I can tell, this aggressive deflation
   method is considered the best QR algorithm currently
   available. If any experts are reading this please correct
   me!

So, the code in my gsl patch is equivalent to the HQR routine
in EISPACK, which has perfectly fine accuracy but is not
the fastest converging. I believe ultimately, the goal should
be to implement method #3 which will perform better than
LAPACK's DGEES, but is fairly complex and will take me some
time to learn/implement. Furthermore, the implementation of
methods 2 or 3 require the double-shift method 1 anyway, so
it was not a waste of time to code up the double-shift method.

I propose adding the double-shift method to gsl (assuming it
passes any further tests needed) which would be a perfectly
fine eigenvalue solver (EISPACK used it for a number of years)
until I (or someone else) can find the time to implement method 3.

Attached is the latest (final? :)) patch.

Patrick Alken

[-- Attachment #2: gsl.unsymm.patch --]
[-- Type: text/plain, Size: 25732 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-25 11:24:48.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-25 11:24:48.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-06-12 15:03:26.000000000 -0600
@@ -58,6 +58,18 @@
 int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);
 
 typedef struct {
+  size_t size;    /* size of matrices */
+  unsigned int max_iterations; /* max iterations since last eigenvalue found */
+  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;
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-25 11:24: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-06-12 15:05:30.000000000 -0600
@@ -0,0 +1,552 @@
+/* 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"
+
+/*
+ * 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 size_t 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->max_iterations = 30 * 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);
+
+      nit = 0;
+
+      /*
+       * compute Schur decomposition of A and store eigenvalues
+       * into eval
+       */
+      schur_decomp(A, 0, N - 1, eval, 0, &nit, w);
+
+      if (nit > w->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)++ < w->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);
+      q = zero_subdiag_small_elements(&m.matrix);
+
+      if (q == 0)
+        {
+          /* no small subdiagonal element found */
+          continue;
+        }
+
+      if (q == (N - 1))
+        {
+          /*
+           * the last subdiagonal 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
+        {
+          /*
+           * 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}|)
+
+Inputs: A - matrix (must be at least 3x3)
+
+Return: row index of small subdiagonal element or 0 if not found
+*/
+
+inline static size_t
+zero_subdiag_small_elements(gsl_matrix * A)
+{
+  const size_t N = A->size1;
+  size_t i;
+  double dpel = gsl_matrix_get(A, N - 2, N - 2);
+
+  for (i = N - 1; i > 0; --i)
+    {
+      double sel = gsl_matrix_get(A, i, i - 1);
+      double del = gsl_matrix_get(A, i, i);
+
+      if ((sel == 0.0) ||
+          (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel))))
+        {
+          gsl_matrix_set(A, i, i - 1, 0.0);
+          return (i);
+        }
+
+      dpel = del;
+    }
+
+  return (0);
+}
+
+/*
+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/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-25 11:24:48.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-25 11:24:48.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-25 11:24:48.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 
 

Re: unsymmetric eigenvalue update

[Brian Gough]

Patrick Alken writes:
 >   I have run my code through Brian Gough's test program
 > on 12x12 and 200x200 random matrices, as well as some integer
 > matrix cases and find no significant disagreement with LAPACK.
 > .....
 > I propose adding the double-shift method to gsl (assuming it
 > passes any further tests needed) which would be a perfectly
 > fine eigenvalue solver (EISPACK used it for a number of years)
 > until I (or someone else) can find the time to implement method 3.

Excellent work!  Useful summary.
I will add the code at the next opportunity.

-- 
Brian Gough

Re: unsymmetric eigenvalue update

[Brian Gough]

At Tue, 13 Jun 2006 09:56:20 -0600,
> I propose adding the double-shift method to gsl (assuming it
> passes any further tests needed) which would be a perfectly
> fine eigenvalue solver (EISPACK used it for a number of years)
> until I (or someone else) can find the time to implement method 3.
> 
> Attached is the latest (final? :)) patch.

The latest patch is in sources.redhat.com CVS with some extra
cleanups.  It works well.  I modified the balance() routine to return
the diagonal scaling elements D (for computing the eigenvectors later)
and moved the memory allocation up out of hessenberg() into the
workspace.

-- 
Brian Gough

Re: unsymmetric eigenvalue update

[Patrick Alken]

Thanks for the update. I have a new balance and hessenberg
routine too that I fixed up a while back. Its no longer necessary
to do any dynamic allocation for the hessenberg function.
I also modified the francis routine a little more to
work better on ill-conditioned matrices and reduce more roundoff
error. Also, I modified linalg.texi and eigen.texi to give
complete documentation on all the functions I've added.
I will make a new patch against the latest CVS and post it as soon as
I can.

An update regarding the high-performance eigensolver: I've
been working on implementing the high-performance small bulge
aggressive early deflation method, and currently have a completely
working code. Unfortunately, the performance is currently about
twice as slow as the fortran version which is going to be in the
next lapack release. I already made gsl_matrix_{get,set} macros
and disabled range checking which resulted in about a 30% speed
improvement. I also made optimized versions of the householder
transform routines specifically for 3-by-1 vectors, but it still
does not match the speed of the fortran version. I am still
investigating why.

Another issue here is that I had to port about 4 lapack routines
over to GSL to get the high-performance eigensolver to work.
Plus the code for the eigensolver is heavily based on the fortran
code which will be in lapack. I know there have been discussions
on this list in the past regarding lapack licensing issues and
I confess I don't fully understand whether it is acceptable/legal for
GSL to contain lapack ported code.

Due to the sheer amount of code needed to get the high-performance
eigensolver to work (currently about 4500 lines) mainly due to
the large lapack routines, I'm starting to wonder if its worth it.
Most of this code is only useful in terms of the eigenvalue problem
and won't be useful for more general linear algebra problems. So
all this code is there for just solving one problem and it might
be wiser to go with the simpler francis solver whose code would
be easier to understand/maintain. The rest of GSL seems to have
the goal of containing small and simple to understand algorithms
and all of this lapack code might just "uglify" the library. Any
comments on this?

Patrick Alken

On Mon, Aug 14, 2006 at 03:24:08PM +0100, Brian Gough wrote:
> At Tue, 13 Jun 2006 09:56:20 -0600,
> > I propose adding the double-shift method to gsl (assuming it
> > passes any further tests needed) which would be a perfectly
> > fine eigenvalue solver (EISPACK used it for a number of years)
> > until I (or someone else) can find the time to implement method 3.
> > 
> > Attached is the latest (final? :)) patch.
> 
> The latest patch is in sources.redhat.com CVS with some extra
> cleanups.  It works well.  I modified the balance() routine to return
> the diagonal scaling elements D (for computing the eigenvectors later)
> and moved the memory allocation up out of hessenberg() into the
> workspace.
> 
> -- 
> Brian Gough
> 

Re: unsymmetric eigenvalue update

[Patrick Alken]

As promised, I made a new patch with some updates. I have
posted the patch and a test program at:

http://ucsu.colorado.edu/~alken

The patch is against the latest CVS.

Changes over the CVS version include:

* updates to linalg.texi and eigen.texi
* support for computing the full Schur form T and Schur vectors Z
  (this will be needed later for computing eigenvectors) Use the
  -z option in the test program to test this

Patrick Alken

On Mon, Aug 14, 2006 at 08:59:25AM -0600, Patrick Alken wrote:
> Thanks for the update. I have a new balance and hessenberg
> routine too that I fixed up a while back. Its no longer necessary
> to do any dynamic allocation for the hessenberg function.
> I also modified the francis routine a little more to
> work better on ill-conditioned matrices and reduce more roundoff
> error. Also, I modified linalg.texi and eigen.texi to give
> complete documentation on all the functions I've added.
> I will make a new patch against the latest CVS and post it as soon as
> I can.
> 
> An update regarding the high-performance eigensolver: I've
> been working on implementing the high-performance small bulge
> aggressive early deflation method, and currently have a completely
> working code. Unfortunately, the performance is currently about
> twice as slow as the fortran version which is going to be in the
> next lapack release. I already made gsl_matrix_{get,set} macros
> and disabled range checking which resulted in about a 30% speed
> improvement. I also made optimized versions of the householder
> transform routines specifically for 3-by-1 vectors, but it still
> does not match the speed of the fortran version. I am still
> investigating why.
> 
> Another issue here is that I had to port about 4 lapack routines
> over to GSL to get the high-performance eigensolver to work.
> Plus the code for the eigensolver is heavily based on the fortran
> code which will be in lapack. I know there have been discussions
> on this list in the past regarding lapack licensing issues and
> I confess I don't fully understand whether it is acceptable/legal for
> GSL to contain lapack ported code.
> 
> Due to the sheer amount of code needed to get the high-performance
> eigensolver to work (currently about 4500 lines) mainly due to
> the large lapack routines, I'm starting to wonder if its worth it.
> Most of this code is only useful in terms of the eigenvalue problem
> and won't be useful for more general linear algebra problems. So
> all this code is there for just solving one problem and it might
> be wiser to go with the simpler francis solver whose code would
> be easier to understand/maintain. The rest of GSL seems to have
> the goal of containing small and simple to understand algorithms
> and all of this lapack code might just "uglify" the library. Any
> comments on this?
> 
> Patrick Alken
> 
> On Mon, Aug 14, 2006 at 03:24:08PM +0100, Brian Gough wrote:
> > At Tue, 13 Jun 2006 09:56:20 -0600,
> > > I propose adding the double-shift method to gsl (assuming it
> > > passes any further tests needed) which would be a perfectly
> > > fine eigenvalue solver (EISPACK used it for a number of years)
> > > until I (or someone else) can find the time to implement method 3.
> > > 
> > > Attached is the latest (final? :)) patch.
> > 
> > The latest patch is in sources.redhat.com CVS with some extra
> > cleanups.  It works well.  I modified the balance() routine to return
> > the diagonal scaling elements D (for computing the eigenvectors later)
> > and moved the memory allocation up out of hessenberg() into the
> > workspace.
> > 
> > -- 
> > Brian Gough
> > 
> 

Re: unsymmetric eigenvalue update

[Brian Gough]

At Mon, 14 Aug 2006 08:59:25 -0600,
Patrick Alken wrote:
> Thanks for the update. I have a new balance and hessenberg
> routine too that I fixed up a while back. Its no longer necessary
> to do any dynamic allocation for the hessenberg function.
> I also modified the francis routine a little more to
> work better on ill-conditioned matrices and reduce more roundoff
> error. Also, I modified linalg.texi and eigen.texi to give
> complete documentation on all the functions I've added.
> I will make a new patch against the latest CVS and post it as soon as
> I can.

Thanks, I downloaded it and am integrating it.

> Due to the sheer amount of code needed to get the high-performance
> eigensolver to work (currently about 4500 lines) mainly due to
> the large lapack routines, I'm starting to wonder if its worth it.
> Most of this code is only useful in terms of the eigenvalue problem
> and won't be useful for more general linear algebra problems. So
> all this code is there for just solving one problem and it might
> be wiser to go with the simpler francis solver whose code would
> be easier to understand/maintain. The rest of GSL seems to have
> the goal of containing small and simple to understand algorithms
> and all of this lapack code might just "uglify" the library. Any
> comments on this?

I think the francis solver is the right tradeoff of complexity vs
efficiency for GSL, it will keep the maintenance costs down and will
be fine for most users.  The bulge solver would be best as a separate
package if you want to distribute it so people can try it out
(http://www.gnu.org/software/gsl/#extensions)

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/