Invalid (?) return values from DGGEV

[pav]

Hi,

I've been puzzling over curious return values from DGGEV for some input data.

The DGGEV man page says:

Code: Select all

               (ALPHAR(j)  + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the gen‐
               eralized eigenvalues.  If ALPHAI(j)  is  zero,  then  the  j-th
               eigenvalue  is  real;  if  positive, then the j-th and (j+1)-st
               eigenvalues are a complex conjugate pair, with ALPHAI(j+1) neg‐
               ative.


However, for certain inputs I get an ALPHAI vector of the form

Code: Select all

ALPHAI = (..., 0, -1.12323e-12, ...)
BETA = (..., 0, 1.12323e-12, ...)


ie. the first eigenvalue of the complex pair has been rounded to zero, which, reading the man page, implies that it is a real eigenvalue. But does the following negative value reliably indicate a complex eigenvalue pair, and is the eigenvector stored in VR(I) now

• VR(I) is the real eigenvector corresponding to ALPHA(I), or,
• VR(I), VR(I+1) comprise the eigenvector corresponding to ALPHA(I) and ALPHA(I+1)

and is the above the intended behavior of LAPACK's DGGEV?

There's a small FORTRAN program reproducing the above condition. The parameters at which the above occurs appear to depend on either the Lapack version, the platform, or the compiler. (I've tested Lapack 3.1.1/gfortran as shipped with Ubuntu intrepid and Lapack 3.0/g77 as shipped with Debian etch). My computer gives the output

Code: Select all

 Checking for missing eigenvalue. If you see no
 further output, you don't have this problem.
 Missing (?) complex pair for omega =  0.25252525252525254     
           1   0.0000000000000000        2.7946960628771595       0.85923885162993430     
           2   0.0000000000000000       -2.7946960628771600       0.85923885162993441     
           3   0.0000000000000000        0.0000000000000000        0.0000000000000000     
           4   0.0000000000000000      -2.18645247810490925E-016  7.95804394604360255E-017
 Missing (?) complex pair for omega =   1.1111111111111112     
           1   0.0000000000000000        0.0000000000000000        0.0000000000000000     
           2   0.0000000000000000      -1.41964227796171439E-016  3.45318391936633162E-017
           3   0.0000000000000000        3.7392934692241413        1.9796259542951338     
           4   0.0000000000000000       -3.7392934692241413        1.9796259542951338     
 Missing (?) complex pair for omega =   1.5656565656565657     
           1   0.0000000000000000       5.67466488887219959E-016  1.24290226548307015E-016
           2   0.0000000000000000      -9.21955150411032210E-017  2.01932654625425201E-017
           3   0.0000000000000000        0.0000000000000000        0.0000000000000000     
           4   0.0000000000000000      -2.62037428091521910E-016  1.82688066063807497E-016
 Missing (?) complex pair for omega =   2.6767676767676769     
           1   0.0000000000000000        1.0350059973275265       0.18232312052566749     
           2   0.0000000000000000       -1.0350059973275265       0.18232312052566749     
           3   0.0000000000000000        0.0000000000000000        0.0000000000000000     
           4   0.0000000000000000      -1.21185328887430225E-016  3.74917111245487433E-016


This problem actually occurred in the Scipy Python library (which wraps LAPACK DGGEV) where the unexpected return value triggered an subsequent error. I'd like to know whether we should (1) treat the above as a LAPACK error condition, (2) try to restore the "missing" eigenvalue (we know it's there because the next ALPHAI is negative) and treat the eigenvector as complex, or (3) leave the "missing" eigenvalue as-is, but treat the eigenvector as complex. Advice would be appreciated.

(Sorry, need to paste the test program here; I don't seem able to attach files on this forum.)

Code: Select all

      PROGRAM TEST
      IMPLICIT NONE
      INTEGER I

      WRITE(*,*) 'Checking for missing eigenvalue. If you see no'
      WRITE(*,*) "further output, you don't have this problem."

      DO 10 I = 0, 99
         CALL CHECK(5.0D0 * I / 99)
 10   CONTINUE
C
C     On my machine, I get missing complex pairs for
C
C        Omega = 25/99, 110/99, 155/99, 265/99
C
C     Other people reported seeing the same at least for Omega = 165/99
C
      END 

      SUBROUTINE CHECK(OMEGA)
      IMPLICIT NONE
      INTEGER INFO, N, I, J
      DOUBLE PRECISION A(4,4), B(4,4), ALPHAR(4), ALPHAI(4),
     1                 BETA(4), VL(4,4), VR(4,4), WORK(8*4), OMEGA

      CALL MATRICES(A, B, OMEGA)
      INFO = 1
      CALL DGGEV('N', 'N', 4, A, 4, B, 4, ALPHAR, ALPHAI, BETA,
     1     VL, 4, VR, 4, WORK, 8*4, INFO)
 
      DO 10 I = 1, 3
         IF (ALPHAI(I) == 0 .AND. ALPHAI(I+1) < 0) THEN
            WRITE(*,*) 'Missing (?) complex pair for omega =', OMEGA
            DO 20 J = 1, 4
               WRITE(*,*) J, ALPHAR(J), ALPHAI(J), BETA(J)
 20         CONTINUE
         END IF
 10   CONTINUE

      RETURN
      END
     
      SUBROUTINE MATRICES(A, B, OMEGA)
      IMPLICIT NONE
      INTEGER I, J
      DOUBLE PRECISION A(4,4), B(4,4), OMEGA

     
      DO 30 I = 1, 4
        DO 40 J = 1, 4
           A(I,J) = 0
           B(I,J) = 0
 40     CONTINUE
 30   CONTINUE

      A(1,1) = 1
      A(2,2) = 1
      A(3,3) = -9 + OMEGA**2
      A(4,4) = -9 + OMEGA**2

      B(1,3) = 1
      B(2,4) = 1
      B(3,1) = 1
      B(4,2) = 1
      B(3,4) = -2*OMEGA
      B(4,3) = 2*OMEGA
     
      RETURN
      END