generalized eigensystems

generalized eigensystems

[Patrick Alken]

Hello all,

  The current CVS contains support for generalized eigensystems.
All the functions are documented in eigen.texi. I have tested
the code for randomized and integer matrix systems up to N = 1000
and compared to lapack and things look good. All my development/testing
was done on amd32/linux so if anyone is willing to test it out on
other archs I'd appreciate it. There is a test program in
eigen/testgen.c which will generate matrices, compute eigenvalues,
and test them against lapack's output. The header of testgen.c
explains how to compile it. Once its compiled, you can just do:

./testgen -n 100 -z

to test N = 100 matrices and compute the schur vectors and test them
too. Sometimes it will output a small discrepancy with lapack
(say around a 10^{-5} difference for an eigenvalue) and most of
these I have tracked down to an initial near-singular or ill-conditioned
matrix input in which case lapack's output can't be trusted either.
But you can send along any output from the test program to me and
I'll have a look at it.

Thanks,
Patrick Alken

Re: generalized eigensystems

[Brian Gough]

At Mon, 12 Mar 2007 23:03:49 -0600,
Patrick Alken wrote:
>   The current CVS contains support for generalized eigensystems.
> All the functions are documented in eigen.texi. I have tested
> the code for randomized and integer matrix systems up to N = 1000
> and compared to lapack and things look good. All my development/testing
> was done on amd32/linux so if anyone is willing to test it out on
> other archs I'd appreciate it.

Good work! I have done some memory testing with valgrind on testgen
and not found any problems. I will try some numerical tests as well
now.

-- 
Brian Gough

Re: generalized eigensystems

[Patrick Alken]

> Good work! I have done some memory testing with valgrind on testgen
> and not found any problems. I will try some numerical tests as well
> now.

I recently realized that when compiling gen.c with -O2 on gcc 3.4.6
it gives incorrect eigenvalue output on some of the integer test
cases. This seems to be a bug in the optimizer since its fine without
the -O2. Upgrading to gcc 4.1.2 seems to fix the problem.

Patrick Alken

Re: generalized eigensystems

[Leopold Palomo Avellaneda]

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

A Dimarts 20 Març 2007 20:17, Patrick Alken va escriure:
> > Good work! I have done some memory testing with valgrind on testgen
> > and not found any problems. I will try some numerical tests as well
> > now.
>
> I recently realized that when compiling gen.c with -O2 on gcc 3.4.6
> it gives incorrect eigenvalue output on some of the integer test
> cases. This seems to be a bug in the optimizer since its fine without
> the -O2. Upgrading to gcc 4.1.2 seems to fix the problem.
>

with 
./testgen -n 100 -z

how long should be the test in your computer? amd32/linux

because I have run the test and it has spent about 2 hours without any finish. 
I have put it run again and redirect the output to a log file ....

Regards,

Leo

Pd this night my computer will work ....

-- 
--
Linux User 152692
PGP: 0xF944807E
Catalonia

[-- Attachment #2: Type: application/pgp-signature, Size: 189 bytes --]

Re: generalized eigensystems

[Patrick Alken]

> with 
> ./testgen -n 100 -z
> 
> how long should be the test in your computer? amd32/linux
> 
> because I have run the test and it has spent about 2 hours without any finish. 
> I have put it run again and redirect the output to a log file ....

It runs forever - it just keeps generating random matrices and testing
them, so if there is no output after 2 hours thats a good sign.

If you want it to stop, you can do:

./testgen -n 100 -z -c 1000

and it will test 1000 matrices and then quit.

Thanks!
Patrick Alken

Re: generalized eigensystems

[Leopold Palomo Avellaneda]

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

A Dimecres 21 Març 2007 01:01, Patrick Alken va escriure:
> > with
> > ./testgen -n 100 -z
> >
> > how long should be the test in your computer? amd32/linux
> >
> > because I have run the test and it has spent about 2 hours without any
> > finish. I have put it run again and redirect the output to a log file
> > ....
>
> It runs forever - it just keeps generating random matrices and testing
> them, so if there is no output after 2 hours thats a good sign.
>
> If you want it to stop, you can do:
>
> ./testgen -n 100 -z -c 1000
>
> and it will test 1000 matrices and then quit.

ok,

after run all the night I have a log file (compressed with bz2) of 10Mb. Do 
you want it?

Also, I understand that the -n parameter show the dimension. Isn't it? How big 
have you test it?

Regards,

Leo

-- 
--
Linux User 152692
PGP: 0xF944807E
Catalonia

[-- Attachment #2: Type: application/pgp-signature, Size: 189 bytes --]

Re: generalized eigensystems

[Brian Gough]

At Mon, 12 Mar 2007 23:03:49 -0600,
Patrick Alken wrote:
> But you can send along any output from the test program to me and
> I'll have a look at it.

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?

-- 
Brian Gough

$ ./a.out
gsl: gen.c:373: ERROR: gen_schur_standardize2 failed on complex block
Default GSL error handler invoked.
Aborted
$ GSL_IEEE_MODE=double-precision ./a.out
GSL_IEEE_MODE="double-precision,trap-common"
solution:
-8.817336437946092624e-16  0.000000000000000000e+00  6.919269373264336664e+00
 1.343730449186862330e+01  0.000000000000000000e+00  1.349256613371829339e+01
 6.035613520034543528e-08  0.000000000000000000e+00  8.100863692583372355e+00
-9.140256551631864568e-08  0.000000000000000000e+00  1.226784413824448983e+01

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

int main ()
{
  double AA[16]   = {   8,  -5,  -5,  -7,
                      -10, -10, -10, -10,
                      -10, -10, -10, -10,
                      -10, -10, -10, -10 };
  double BB[16]   = {  7, -10,   -9,   -1,
                       4,   3,   -7,    9,
                      -4,   1,   -9,    9,
                      -7,   5,   -5,   -1} ;

  gsl_matrix_view A = gsl_matrix_view_array(AA, 4, 4);
  gsl_matrix_view B = gsl_matrix_view_array(BB, 4, 4);

  gsl_vector_complex * alpha = gsl_vector_complex_alloc(4);
  gsl_vector * beta = gsl_vector_alloc(4);

  gsl_eigen_gen_workspace * w = gsl_eigen_gen_alloc(4);
  gsl_eigen_gen_params (0, 0, 0, w);

  gsl_ieee_env_setup();
  gsl_eigen_gen (&A.matrix, &B.matrix, alpha, beta, w);

  printf("solution:\n");

  {
    gsl_vector_view a_real = gsl_vector_complex_real (alpha);
    gsl_vector_view a_imag = gsl_vector_complex_imag (alpha);
    int k;

    for (k = 0; k < 4; k++) { 
      printf("% .18e % .18e % .18e\n", 
             gsl_vector_get(&a_real.vector,k), 
             gsl_vector_get(&a_imag.vector,k),
             gsl_vector_get(beta,k));
    }
  }
}

Re: generalized eigensystems

[Patrick Alken]

> ok,
> 
> after run all the night I have a log file (compressed with bz2) of 10Mb. Do 
> you want it?
> 
> Also, I understand that the -n parameter show the dimension. Isn't it? How big 
> have you test it?

I've tested it overnight with -n 1000 and found no problems. Could
you put your file on the web somewhere where I could download it?
putstuff.putfile.com is a public file hosting site where its
easy to upload files.

Thanks,
Patrick Alken

Re: generalized eigensystems

[Leopold Palomo-Avellaneda]

A Dimecres 21 Març 2007 15:47, Patrick Alken va escriure:
> > ok,
> >
> > after run all the night I have a log file (compressed with bz2) of 10Mb.
> > Do you want it?
> >
> > Also, I understand that the -n parameter show the dimension. Isn't it?
> > How big have you test it?
>
> I've tested it overnight with -n 1000 and found no problems. Could
> you put your file on the web somewhere where I could download it?
> putstuff.putfile.com is a public file hosting site where its
> easy to upload files.

Well, 

this morning I have launch :
./testgen -n 1000 -z 2&1> test_big2.log &  

now the process is:


 
  PID USER      PR  NI  VIRT  RES  SHR S %CPU %MEM    TIME+  COMMAND
 9688 leo       25   0  126m  64m 2832 R 99.9  6.5 301:22.14 testgen

and test_big2.log is empty ...

Regards,

Leo

Re: generalized eigensystems

[Patrick Alken]

[-- 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;
}

Re: generalized eigensystems

[Patrick Alken]

> 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).

Well, I've found out that the algorithm is failing under -O2 when
computing the eigenvalues of a 2-by-2 generalized block. To
do this I had copied the equations directly from Moler's paper
which is the same thing EISPACK did, but LAPACK has a fancier
routine which computes scaling factors and checks for overflow/underflow
etc. Apparantly the bug shows up when looking at a block with
"almost" complex eigenvalues - ie: a real eigenvalue block with
a determinant very close to 0. In the no-optimization case it
computes the determinant as slightly positive, and with optimization
on it gets computed as slightly negative which must be the result
of some sort of underflow condition.

Looks like I'll have to port over the lapack routine...bleh. I should
have a patch ready within the next week or so.

Patrick Alken

Re: generalized eigensystems

[Brian Gough]

At Wed, 21 Mar 2007 12:02:06 -0600,
Patrick Alken wrote:
> 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.

You can use GSL_IEEE_MODE=double-precision,mask-all to avoid that, or
just call finite() directly.

It does look like the order of the #ifdefs in sys/infnan.c could be
improved there to pick up the right functions though.

-- 
Brian Gough

Re: generalized eigensystems

[Patrick Alken]

> 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).

Sorry for the delay - the last few weeks have been busy. I believe
I've solved all the -O2 problems and the single/double precision
issues now. The latest CVS has the updates.

Patrick Alken

Re: generalized eigensystems

[Brian Gough]

At Fri, 30 Mar 2007 15:33:59 -0600,
Patrick Alken wrote:
> Sorry for the delay - the last few weeks have been busy. I believe
> I've solved all the -O2 problems and the single/double precision
> issues now. The latest CVS has the updates.

Thanks. I've added some changes to make the gsl_isinf/isnan/finite
functions always default to the system versions where available. 

-- 
Brian Gough

Re: generalized eigensystems

[Patrick Alken]

The latest CVS has support for generalized nonsymmetric eigenvectors, as
well as a solver for generalized symmetric-definite eigensystems (which
trivially reduce to symmetric eigensystems). These updates make
the gsl eigen package basically feature-equivalent to LAPACK, with
the exceptions of left eigenvectors and condition numbers.

There is a new testgen program (the old one is under testgen2.c). I'm
still working on the testgen program but its functional for now.

Patrick Alken