Date: Thu, 22 Jan 2009 02:04:47 -0700
From: Vasile Sima
To: "lapack@cs.utk.edu" <lapack@cs.utk.edu>
Cc: Pascal Gahinet, Peter Benner
Subject: [Lapack] Bug in DHGEQZ / SHGEQZ

Hello,

[...]

Another issue with DHGEQZ routine is the possible non-convergence for certain
matrix pencils.  I reported

this problem in January 2003, and David Day proposed a slightly modified
version, which always converged in my tests.  The calculation of the
exceptional shifts is different.

I am attaching an archive containing a modified file dhgeqz.f, including
Day's modification.  The code differences compared to the posted LAPACK
version are shown below:

$ diff dhgeqz.f lapack-3/lapack-3.2/src.
629,634c629,635
<             WR = ABS( H( ILAST-1, ILAST-2 ) / T( ILAST-1, ILAST-1 ) ) +
<      $           ABS( H( ILAST,   ILAST-1 ) / T( ILAST-2, ILAST-2 ) )
<             ESHIFT = ABS( H( ILAST-1, ILAST-2 ) /
<      $                    T( ILAST-2, ILAST-2 ) ) +
<      $               ABS( H( ILAST, ILAST-1 ) / T( ILAST-1, ILAST-1 ) )
<             ESHIFT = MIN( ESHIFT, WR )
---
>             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
>      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
>                ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
>      $                  T( ILAST-1, ILAST-1 )
>             ELSE
>                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
>             END IF
780c780
<      $                    T( ILAST, ILAST+1 ), LDT, CL, SL )
---
>      $                    T( ILAST, ILAST+1 ), LDH, CL, SL )

The file in the archive also contains the previous lines for the
computation of the exceptional shifts, in comments.  Of course, there
is no guarantee that the modified code will always converge.  Perhaps,
a better approach would be to use the original strategy at iteration
count 10, and the modified strategy at iteration 20, or similar.
I made no such experiments.

Besides, the archive contains two test files (one for DGGES, and the
other for DHGEQZ) and associated data and results.  DGGES failed even
for problems of order 4, with

       0.    1.    0.    0.
       0.    0.    0.   -1.
A =   -1.    0.    0.    0.
       0.    0.   -1.    0.

and B = I_4 (or a similar problem with changed signs for the last two
rows of A and B).  While this is actually a standard eigenproblem,
DGGES should be able to solve it.  But, DGGES returns with the

error message

 INFO on exit from DGGES =  4

and therefore, no eigenvalue was found.  The non-convergence appear in
DHGEQZ.

The DHGEQZ.dat file contains the data for the example for which I found
the bug, but simpler examples can be constructed;  the bug is clear,
but not so easy to find.

[...]

Sincerely yours,
Vasile Sima