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

[Andries E. Brouwer]

On Wed, Aug 29, 2007 at 08:31:02PM +0100, Brian Gough wrote:
> 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.

Very good, that solves the loop. Thanks!

Remains my question: what about the accuracy?
Is there any guarantee, or are the eigenvalues just reasonable guesses
with unknown error bounds?

I see 4.251127, 2.291957, 1.294160 where the exact values in 5 decimal
places are 4.24264, 2.30278, 1.30278.

Andries