Date: Sat, 6 Dec 2008 08:07:58 -0700 From: Pauli Virtanen To: "lapack@cs.utk.edu" Cc: "scipy-dev@scipy.org" Subject: [Lapack] Invalid (?) output from DGGEV for some matrices Dear Lapack authors, DGGEV returns for some matrices an ALPHAI vector such that ALPHAI(i) = 0 but ALPHAI(i+1) < 0. The DGGEV documentation does not mention such outputs. Does this imply that the i, i+1 entries describe a complex eigenvalue pair, or does this indicate a bug in LAPACK or in the compiler? The problematic matrix pairs are of the form [ 1 0 0 0 ] [ 0 0 1 0 ] c1 = omega^2 - 9 [ 0 1 0 0 ] [ 0 0 0 1 ] c2 = 2 omega [ 0 0 c1 0 ] [ 1 0 0 -c2 ] [ 0 0 0 c1 ] [ 0 1 c2 0 ] A test program `2dof.f` illustrating this is attached. However, what it outputs appears to depend on compiler, compiler version, and optimization flags: * LAPACK 3.2, compiled with gfortran 4.3.2-1ubuntu11, OPTS = -O0: No problematic outputs. BLAS provided by Atlas 3.6.0. * LAPACK 3.2, compiled with gfortran 4.3.2-1ubuntu11, OPTS = -O1 (2/3): Missing (?) complex pair for omega = 1.5656565656565657 Info: 0 I / ALPHAR(I) / ALPHAI(I) / BETA(I): 1 0.0000000000000000 1.7080358288156292 0.37410519259457370 2 0.0000000000000000 -1.7080358288156292 0.37410519259457370 3 0.0000000000000000 0.0000000000000000 0.0000000000000000 4 0.0000000000000000 -2.62037428091521910E-016 1.82688066063807497E-016 Missing (?) complex pair for omega = 1.9191919191919191 Info: 0 I / ALPHAR(I) / ALPHAI(I) / BETA(I): 1 0.0000000000000000 0.0000000000000000 0.0000000000000000 2 0.0000000000000000 -1.46001785734970358E-017 2.96800344717906839E-018 3 0.0000000000000000 3.6846073726768025 3.4091227092990981 4 0.0000000000000000 -3.6846073726768025 3.4091227092990981 BLAS provided by Atlas 3.6.0. * LAPACK 3.2, compiled with gfortran 4.1.2 20061115, OPTS = -O0: Missing (?) complex pair for omega = 1.86868686868687 BLAS provided by Atlas 3.6.0; same with Netlib BLAS. * LAPACK 3.1.1, compiled with g77 3.4.6, OPTS = -O0/O1/O2/O3: Missing (?) complex pair for omega = 1.86868687 BLAS provided by Atlas 3.6.0; same with Netlib BLAS. What appears to matter for gfortran 4.3.2 is which optimization flags were used to compile the file `dhgeqz.f`; recompiling this file can make the problem appear (-O1) and disappear (-O0). Output from LAPACK double eigenvalue tests, `dgd-O0/O1.out` is also attached. It corresponds to gfortran 4.3.2-1ubuntu11, for -O0 and -O1. There are some failures, but they are the same for both flags. This issue is present at least in LAPACK libraries shipped by Ubuntu and Debian. It may be that also other Linux distributions are affected. *** We ran into this problem while wrapping DGGEV for use in the higher-level language Python: the ambiguity in the return value caused an error in repacking the complex eigenvectors in our code. [1] (Cc'd to the developer list.) If this is a LAPACK or compiler bug, suggestions on how the DGGEV return values should be interpreted are appreciated: should we treat this as an error condition, or should we attempt to treat the result as a complex eigenvalue pair? Best regards, Pauli Virtanen .. [1] http://scipy.org/scipy/scipy/ticket/709