LAPACK Archives

[Lapack] Bug in DHGEQZ / SHGEQZ

Hello,

I am an enthusiastic user of LAPACK.  This time, I would like to draw
your attention to a bug in the file dhgeqz.f (and also in shgeqz.f)
included in the archive lapack.tgz, updated on November 18, 2008 on the
site http://www.netlib.org/lapack/.  Namely, the argument LDH in line
780 must be replaced by LDT.  The bug is minor, but can lead to bad
results if LDH differs from LDT.  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.  Note that the computed eigenvalues are
correct, but the retuned matrix pair is not in real Schur-triangular
form, so the claim in the comments of the original code is not true.
For some problems, it is important to have the real Schur-triangular
form available.

I used an IBM PC-compatible computer, the Compaq Visual Fortran 6.5
(September 2000), and the compiled or prebuilt LAPACK library
(lapack_win32.lib) from netlib.

Sincerely yours,
Vasile Sima

============================================================
Dr. Vasile Sima,  SMIEEE
National Institute for Research & Development in Informatics
Bd. Maresal Al. Averescu, Nr. 8-10
011455, Bucharest 1,  Romania
Phone:   (+40)-21-316.07.65  (in Romania: 021-316.07.65)
Fax:     (+40)-21-316.10.30  (in Romania: 021-316.10.30)
E-mail:  vsima@Domain.Removed
http://www.ici.ro/ici/homepage/vsima.html
============================================================
-------------- next part --------------
A non-text attachment was scrubbed...
Name: DHGEQZ_er.ZIP
Type: application/x-zip-compressed
Size: 15781 bytes
Desc: not available
Url : 
http://lists.cs.utk.edu/private/lapack/attachments/20090122/ea6d776f/DHGEQZ_er.bin
 

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or