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;
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
