Date: Thu, 22 Jan 2009 02:04:47 -0700 From: Vasile Sima To: "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