Re: Eigenvalues of a Hermitian matrix

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

Andrew W Steiner writes:
 > 	The following code fails to calculate the eigenvalues of the
 > second matrix. (It's only 70 lines, and even that would be
 > shortened if I wouldn't have tried to be so pedantic with my
 > style). Is my usage correct, or is there possibly an error in my
 > compilation of the libraries? I am using Red Hat Linux 6.2 with g++
 > (version egcs-2.91.66).
 > 	Could someone possibly double check this for me? The code is
 > supposed to print out and calculate the eigenvalues of two simple
 > matrices but never leaves the second call to gsl_eigen_herm.

Thanks for the bug report.  Your program is valid. The eigenvalue
routine was going into an infinite loop.

The following patch to linalg/qrstep.c fixes it.  I will check it into
CVS.

After applying the patch,

$ ./a.out 
....
(0 0) (0 0) (-1 0) (0 0) 
(0 0) (1 0) (0 0) (1 0) 
(-1 0) (0 0) (0 0) (0 0) 
(0 0) (1 0) (0 0) (0 0) 
1 -1 1.61803 -0.618034 


Index: qrstep.c
===================================================================
RCS file: /cvs/gsl/gsl/eigen/qrstep.c,v
retrieving revision 1.2
diff -u -r1.2 qrstep.c
--- qrstep.c    2001/06/12 12:21:46     1.2
+++ qrstep.c    2001/08/02 17:15:23
@@ -60,18 +60,13 @@
 
   double mu;
 
-  if (dt > 0)
+  if (dt >= 0)
     {
       mu = tb - (tab * tab) / (dt + hypot (dt, tab));
     }
-  else if (dt < 0)
-    {
-      mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
-    }
   else
     {
-      /* FIXME: what is the best choice of mu for dt == 0 ? */
-      mu = tb;
+      mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
     }
 
   return mu;