Re: [Bug-gsl] Re: infinite loop in gsl_eigen_symm

[Brian Gough]

At Wed, 22 Aug 2007 17:30:18 +0200,
Andries E. Brouwer wrote:
> Investigating a bit closer, I see in symm.c:gsl_eigen_symm()
> a loop while (b > 0) { ... } that hangs (flipflops between two states).
> If I change the test
>       if (sd[b - 1] == 0.0 || ...
> into
>       if ((sd[b - 1] > -5e-188 && sd[b-1] < 5e-188) || ...
> then all is fine.
> 
> So, it seems gsl has assumptions about the arithmetic on very small
> numbers that are false on this particular machine.
> 

I found the problem.  The qrstep uses the square of the offdiagonal
element (O(1e-200^2). This underflows to zero without extended
precision so no progress is made.  The following patch should fix the
problem. It will be included in the next release.  Thanks for the bug
report.

-- 
Brian Gough


Index: qrstep.c
===================================================================
RCS file: /home/gsl-cvs/gsl/eigen/qrstep.c,v
retrieving revision 1.4
diff -u -r1.4 qrstep.c
--- qrstep.c    26 Jul 2003 13:44:33 -0000      1.4
+++ qrstep.c    29 Aug 2007 19:20:55 -0000
@@ -1,7 +1,7 @@
 /* remove off-diagonal elements which are neglegible compared with the
    neighboring diagonal elements */
 
-inline static void
+static void
 chop_small_elements (const size_t N, const double d[], double sd[])
 {
   double d_i = d[0];
@@ -60,13 +60,17 @@
 
   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));
+      mu = tb - fabs(tab);
     }
   else
     {
-      mu = tb + (tab * tab) / ((-dt) + hypot (dt, tab));
+      mu = tb + tab * (tab / ((-dt) + hypot (dt, tab)));
     }
 
   return mu;