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

[Andries E. Brouwer]

On Tue, Aug 28, 2007 at 12:16:41AM +0200, Oliver Jennrich wrote:
>> 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.
> 
> 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) )

Note that my above "all is fine" only refers to the looping behaviour
that is avoided by this test. The results are bad, in the sense that
the computed eigenvalues already have errors in the second decimal
(e.g. 2.292 instead of 2.303).

I have not seen much documentation for these functions - don't know whether
there are any claims about a guaranteed accuracy. Is this basically a
function that you give a matrix and it returns some numbers and nobody
knows how far the real eigenvalues are from the returned numbers?

Andries

(It is easy to compute the exact eigenvalues in the example I gave.
They are 0 with multiplicity 21, +-3sqrt{3}, +-(1+-sqrt{13})/2.)