Re: generalized eigensystems

[Patrick Alken <patrick.alken@colorado.edu>]

[-- Attachment #1: Type: text/plain, Size: 1685 bytes --]

> 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 #2: testfin.c --]
[-- Type: text/plain, Size: 217 bytes --]

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>

int
main()
{
  double a;
  int s;

  gsl_ieee_env_setup();

  a = GSL_POSINF;
  s = gsl_finite(a);
  printf("s = %d\n", s);

  return 0;
}