* generalized eigensystems @ 2007-03-13 5:03 Patrick Alken 2007-03-16 15:49 ` Brian Gough 2007-03-21 10:10 ` Brian Gough 0 siblings, 2 replies; 15+ messages in thread From: Patrick Alken @ 2007-03-13 5:03 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-13 5:03 generalized eigensystems Patrick Alken @ 2007-03-16 15:49 ` Brian Gough 2007-03-20 19:17 ` Patrick Alken 2007-03-21 10:10 ` Brian Gough 1 sibling, 1 reply; 15+ messages in thread From: Brian Gough @ 2007-03-16 15:49 UTC (permalink / raw) To: Patrick Alken; +Cc: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-16 15:49 ` Brian Gough @ 2007-03-20 19:17 ` Patrick Alken 2007-03-20 22:35 ` Leopold Palomo Avellaneda 0 siblings, 1 reply; 15+ messages in thread From: Patrick Alken @ 2007-03-20 19:17 UTC (permalink / raw) To: gsl-discuss > 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-20 19:17 ` Patrick Alken @ 2007-03-20 22:35 ` Leopold Palomo Avellaneda 2007-03-21 0:01 ` Patrick Alken 0 siblings, 1 reply; 15+ messages in thread From: Leopold Palomo Avellaneda @ 2007-03-20 22:35 UTC (permalink / raw) To: gsl-discuss [-- 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 --] ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-20 22:35 ` Leopold Palomo Avellaneda @ 2007-03-21 0:01 ` Patrick Alken 2007-03-21 7:32 ` Leopold Palomo Avellaneda 0 siblings, 1 reply; 15+ messages in thread From: Patrick Alken @ 2007-03-21 0:01 UTC (permalink / raw) To: gsl-discuss > 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 0:01 ` Patrick Alken @ 2007-03-21 7:32 ` Leopold Palomo Avellaneda 2007-03-21 14:47 ` Patrick Alken 0 siblings, 1 reply; 15+ messages in thread From: Leopold Palomo Avellaneda @ 2007-03-21 7:32 UTC (permalink / raw) To: gsl-discuss [-- 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 --] ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 7:32 ` Leopold Palomo Avellaneda @ 2007-03-21 14:47 ` Patrick Alken 2007-03-21 15:02 ` Leopold Palomo-Avellaneda 0 siblings, 1 reply; 15+ messages in thread From: Patrick Alken @ 2007-03-21 14:47 UTC (permalink / raw) To: gsl-discuss > 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 14:47 ` Patrick Alken @ 2007-03-21 15:02 ` Leopold Palomo-Avellaneda 0 siblings, 0 replies; 15+ messages in thread From: Leopold Palomo-Avellaneda @ 2007-03-21 15:02 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-13 5:03 generalized eigensystems Patrick Alken 2007-03-16 15:49 ` Brian Gough @ 2007-03-21 10:10 ` Brian Gough 2007-03-21 18:02 ` Patrick Alken ` (2 more replies) 1 sibling, 3 replies; 15+ messages in thread From: Brian Gough @ 2007-03-21 10:10 UTC (permalink / raw) To: Patrick Alken; +Cc: gsl-discuss 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)); } } } ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 10:10 ` Brian Gough @ 2007-03-21 18:02 ` Patrick Alken 2007-03-21 19:38 ` Brian Gough 2007-03-21 19:18 ` Patrick Alken 2007-03-30 21:34 ` Patrick Alken 2 siblings, 1 reply; 15+ messages in thread From: Patrick Alken @ 2007-03-21 18:02 UTC (permalink / raw) To: gsl-discuss [-- 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; } ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 18:02 ` Patrick Alken @ 2007-03-21 19:38 ` Brian Gough 0 siblings, 0 replies; 15+ messages in thread From: Brian Gough @ 2007-03-21 19:38 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 10:10 ` Brian Gough 2007-03-21 18:02 ` Patrick Alken @ 2007-03-21 19:18 ` Patrick Alken 2007-03-30 21:34 ` Patrick Alken 2 siblings, 0 replies; 15+ messages in thread From: Patrick Alken @ 2007-03-21 19:18 UTC (permalink / raw) To: gsl-discuss > 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-21 10:10 ` Brian Gough 2007-03-21 18:02 ` Patrick Alken 2007-03-21 19:18 ` Patrick Alken @ 2007-03-30 21:34 ` Patrick Alken 2007-04-03 14:32 ` Brian Gough 2007-04-23 22:15 ` Patrick Alken 2 siblings, 2 replies; 15+ messages in thread From: Patrick Alken @ 2007-03-30 21:34 UTC (permalink / raw) To: gsl-discuss > 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-30 21:34 ` Patrick Alken @ 2007-04-03 14:32 ` Brian Gough 2007-04-23 22:15 ` Patrick Alken 1 sibling, 0 replies; 15+ messages in thread From: Brian Gough @ 2007-04-03 14:32 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: generalized eigensystems 2007-03-30 21:34 ` Patrick Alken 2007-04-03 14:32 ` Brian Gough @ 2007-04-23 22:15 ` Patrick Alken 1 sibling, 0 replies; 15+ messages in thread From: Patrick Alken @ 2007-04-23 22:15 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 15+ messages in thread
end of thread, other threads:[~2007-04-23 22:15 UTC | newest] Thread overview: 15+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2007-03-13 5:03 generalized eigensystems Patrick Alken 2007-03-16 15:49 ` Brian Gough 2007-03-20 19:17 ` Patrick Alken 2007-03-20 22:35 ` Leopold Palomo Avellaneda 2007-03-21 0:01 ` Patrick Alken 2007-03-21 7:32 ` Leopold Palomo Avellaneda 2007-03-21 14:47 ` Patrick Alken 2007-03-21 15:02 ` Leopold Palomo-Avellaneda 2007-03-21 10:10 ` Brian Gough 2007-03-21 18:02 ` Patrick Alken 2007-03-21 19:38 ` Brian Gough 2007-03-21 19:18 ` Patrick Alken 2007-03-30 21:34 ` Patrick Alken 2007-04-03 14:32 ` Brian Gough 2007-04-23 22:15 ` Patrick Alken
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).