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

[Oliver Jennrich]

On 8/27/07, Brian Gough <address@hidden> wrote:
> At Wed, 22 Aug 2007 17:30:18 +0200,
> Andries E. Brouwer wrote:
> > > In a program that has called gsl_eigen_symm() successfully 10^9 times,
> > > I find that on one specific matrix this routine hangs in an infinite loop.
> > > Trying this same matrix on a different architecture, all is well.
> > >
> > > The machine with infinite loop is an x86_64:
> > > Intel(R) Core(TM)2 CPU 6300 @ 1.86GHz
> > >
> > > A test program that loops is given below.
> >
> > 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.
>
> Thanks for the bug report.  Can you send me a complete list of the
> compiler version, operating system version and options used to compile
> GSL so I can try to reproduce the problem.
>
> I expect it is caused by different underflow behaviour with some
> recent 64bit or SSE instructions -- I haven't been able to reproduce
> it compiling with -msse2 -fpmath=sse on a 32-bit machine.
>

Well, I'm wondering if testing a float on an exact zero is in
principle the right thing to do. It seems as all sort of rounding
errors are invited in that are very hard to track down.

Wouldn't it be better to have something like GSL_IEEE_TEST_ZERO( x ) with
#define GSL_IEEE_TEST_ZERO(x)  ( (x < GSL_MIN_IEEE) && ( x> -GSL_MIN_IEEE) )


-- 
Space -- the final frontier