Re: Eigenvalues of a Hermitian matrix

[Brian Gough <bjg at network-theory dot co dot uk>]

Brian Gough writes:
 > The following patch to linalg/qrstep.c fixes it.  I will check it into
 > CVS.

There was also a bug in the accumulation of the eigenvectors.  I've
added this to the test suite.

Index: symmv.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/symmv.c,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -r1.2 -r1.3
--- symmv.c	2001/07/01 21:48:38	1.2
+++ symmv.c	2001/08/02 19:10:53	1.3
@@ -192,10 +192,10 @@
                 
                 for (k = 0; k < N; k++)
                   {
-                    double qki = gsl_matrix_get (evec, k, i);
-                    double qkj = gsl_matrix_get (evec, k, i + 1);
-                    gsl_matrix_set (evec, k, i, qki * c - qkj * s);
-                    gsl_matrix_set (evec, k, i + 1, qki * s + qkj * c);
+                    double qki = gsl_matrix_get (evec, k, a + i);
+                    double qkj = gsl_matrix_get (evec, k, a + i + 1);
+                    gsl_matrix_set (evec, k, a + i, qki * c - qkj * s);
+                    gsl_matrix_set (evec, k, a + i + 1, qki * s + qkj * c);
                   }
               }
             
Index: hermv.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/hermv.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- hermv.c	2001/07/01 21:48:38	1.3
+++ hermv.c	2001/08/02 19:10:53	1.4
@@ -215,8 +215,8 @@
                 
                 for (k = 0; k < N; k++)
                   {
-                    gsl_complex qki = gsl_matrix_complex_get (evec, k, i);
-                    gsl_complex qkj = gsl_matrix_complex_get (evec, k, i + 1);
+                    gsl_complex qki = gsl_matrix_complex_get (evec, k, a + i);
+                    gsl_complex qkj = gsl_matrix_complex_get (evec, k, a + i + 1);
                     /* qki <= qki * c - qkj * s */
                     /* qkj <= qki * s + qkj * c */
                     gsl_complex x1 = gsl_complex_mul_real(qki, c);
@@ -228,8 +228,8 @@
                     gsl_complex qqki = gsl_complex_add(x1, y1);
                     gsl_complex qqkj = gsl_complex_add(x2, y2);
                     
-                    gsl_matrix_complex_set (evec, k, i, qqki);
-                    gsl_matrix_complex_set (evec, k, i + 1, qqkj);
+                    gsl_matrix_complex_set (evec, k, a + i, qqki);
+                    gsl_matrix_complex_set (evec, k, a + i + 1, qqkj);
                   }
               }