Re: generalized eigensystems

[Patrick Alken <patrick dot alken at colorado dot edu>]

> Here is an interesting test case, it aborts with a "shouldn't get
> here" error in extended precision but works in double precision (via
> GSL_IEEE_MODE).
> 
> Might be worth also looking at the previous case you had that went
> away when changing the optimisation - did you try it in
> double-precision?

I could not reproduce the exact error you got, but that same test
case is giving incorrect eigenvalues when gsl is compiled with -O2
with gcc 4.1.2 with single precision. When GSL_IEEE_MODE is set to
double-precision, it gives correct output. If I compile libgsleigen
without -O2, it gives correct output in both cases...very odd. I'm still
looking into this. It may be an issue with not properly detecting
when a matrix element has become negligible.

When trying to debug this, I came across the following issue:

when GSL_IEEE_MODE is set to double-precision, gsl_finite(GSL_POSINF)
fails. I have attached a sample test program to demonstrate this:

> gcc --version
gcc (GCC) 4.1.2

> gcc -g -Wall -o testfin testfin.c -lgsl -lgslcblas
> setenv GSL_IEEE_MODE "double-precision"
> ./testfin
GSL_IEEE_MODE="double-precision,trap-common"
Floating exception (core dumped)
> gdb testfin core
...
Program terminated with signal 8, Arithmetic exception.
#0  0xb7ef3578 in gsl_finite (x=inf) at infnan.c:103
103       const double y = x - x;

...

> unsetenv GSL_IEEE_MODE
> ./testfin
s = 0

Is this a known issue? I'm trying to use gsl_finite() in testgen
to better sort the eigenvalues so I can compare them to lapack
if there are infinities present. Basically I'm getting segfaults
all over the place whenever I try to access nans/infs in
double-precision mode.

Patrick Alken

Attachment: testfin.c
Description: Text document