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

[Patrick Alken]

I can confirm the infinite loop problem on:

Linux 2.6.18
gcc 4.1.1
intel pentium dual core 3.2 GHz

standard configure/make for GSL

After briefly looking at the code, I'm not sure what the problem is :).
For this particular matrix, the diagonal is all 0 (and stays very close
to 0 during the QR process. This makes the test in chop_small_elements
always fail.

My best guess is that the problem lies in trailing_eigenvalue. This
function may not handle the case where the lower 2x2 block has entries
which are all close to 0, but I'm not sure. I see that lapack (dsterf)
has a more complicated way of computing that eigenvalue shift.


On Mon, Aug 27, 2007 at 10:21:04PM +0100, Brian Gough 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.
> 
> -- 
> Brian Gough
> 
> GNU Scientific Library -
> http://www.gnu.org/software/gsl/
> 
> 
> _______________________________________________
> Bug-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/bug-gsl
>