This represents known issue with the LAPACK software package that have not been yet solved. (As of Sunday October 26th 2008.)
- LAPACK does not have any exceptional shift strategy for QZ
- DQDS algorithm (useful to compute the singular values of a matrix)
1 -- LAPACK does not have any exceptional shift strategy for QZ[1.1-- The problems] LAPACK current shift strategy is known to fail on various pencils. Here are three pencils, (A1,B1), (A2,B2), (A3,B3) where LAPACK fails A1 = [ 0 0 1 ; 0 1 0 ; 1 0 0 ]; B1 = [ 2 0 0 ; 0 0 1 ; 0 1 0 ]; Pencil (A1,B1) has been sent to us by Zbigniew Leyk, Senior Lecturer, Department of Computer Science, Texas A & M University on Jan 30 2007. See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=317 A2 = [ 0 0 0 1 ; 1 0 0 0 ; 0 1 0 0 ; 0 0 1 0 ]; B2 = [ 1 0 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1 ]; Pencil (A2,B2) has been sent to us by Daniel Kressner in 2005 who quote Vasile Sima. A3 = [ 1 1 0; 0 0 -1; 1 0 0 ]; B3 = [ -1 0 -1; 0 -1 0; -1 0 -1 ]; Pencil (A3,B3) has been sent to us by Patrick Alken from University of Colorado at Boulder on Apr 24 2007. See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=382 A4 = [ 1 1e-9 0; 1e-9 1 0; 0 0 1 ]; B4 = [ 1 0 0; 0 1 0; 0 0 1 ]; Pencil (A4,B4) has been sent to us by Micheal Baudin from the Scilab project on Tu Oct 28th 2008. This one pencil make ZGGEV fail. DGGEV handles it correctly. http://bugzilla.scilab.org/show_bug.cgi?id=3652 [1.2-- The fixes] [1.2.a-- Mathworks' patch] Bobby Cheng sent the lapackers the Mathworks patch to fix LAPACK QZ. This patch works fine on the pencils (A1,B1) and (A2,B2). But fails on pencil (A3,B3). As observed by Jim Demmel, matlab qz works fine with pencil (A3,B3). Date: Thu, 15 Sep 2005 20:35:00 +0200 (CEST) From: Daniel Kressner Subject: Re: [Lapackers] Mathworks modif on LAPACK Hi, below are some comments on Mathworks' modifications to DGEBAL and DHGEQZ. > 1)Change SCLFAC from 8 to 2 to be the same as LINPACK. > > File: > cgebal.f > zgebal.f > sgebal.f > dgebal.f I support the point of view taken by Mathworks. Changing SCLFAC from 8 to 2 may result in a few extra iterations of the balancing algorithm, but it may also give a balanced matrix of slightly smaller norm and consequently one or two more accurate digits in the computed eigenvalues. I am not aware that the execution time needed by balancing could be a major concern when solving nonsymmetric eigenvalues. > 3)Change to implicit shift and exceptional shift calculations. > Some problem were not converging. > > Files: > zhgeqz.f (lines 602-672) > chgeqz.f (lines 602-672) > shgeqz.f (lines 657-667, 819-821) > dhgeqz.f (lines 657-667, 819-821) As far as I can see, Mathworks applied four changes to *HGEQZ (not all of them are marked in the code): (1) Exceptional shift strategy As also pointed out by Vasile Sima, the exceptional shift strategy for the QZ algorithm currently implemented in DHGEQZ seems to be less effective than the one for the QR algorithm. For example, it fails for the matrix pair [ 0 0 0 1 ] [ 1 0 0 0 ] A = [ 0 1 0 0 ], B = identity. [ 0 0 1 0 ] LAPACK's exceptional shift is zero (this is the same as the standard shift -> no convergence). Mathworks' exceptional shift is one (much better). (2) Choice of single real shift DHGEQZ applies a single shift QZ algorithm if the computed shifts are real. This seems to be a relict of Ward's combination shift QZ algorithm, which has no meaningful purpose anymore (applying two single shift QZ iterations is more expensive than applying one double shift QZ iteration, even in terms of flops). If two real shifts are computed, Mathworks' change has the effect that always the shift closer to the last diagonal element of A\B is used. If I remember correctly, this is in the spirit of what Wilkinson proposed for symmetric matrices. I support this change and expect that it can lead to slightly faster convergence. (3) Mathworks included B22 = -B22 on line 820 (in DHGEQZ) This is a bug fix and I highly recommend to include it. It seems that DHGEQZ computes wrong complex conjugate eigenvalue pairs if DLASV2 computes negative diagonal elements (I am not sure whether this is an unlikely event in the context of the QZ algorithm). (4) Mathworks uses the workspace to store the number of iterations needed for each eigenvalue. This is probably used for debugging and should not be included in LAPACK. As announced at the Householder meeting, Bo Kagstrom and I are currently working on an extension of Ralph Byers' multishift aggressive QR code to QZ. This code will take the modifs of Mathworks into account. However, I suggest to include these modifs (particularly (1) and (3)) already in the next update of LAPACK. With best regards, Daniel [1.2.b-- David Day's patch] I found in the lapackers archive a patch from David Day. From: David Day (Sandia) Date: Tue Sep 7 09:49:40 2004 Subject: [Lapackers] Re: [Fwd: Convergence failure of LAPACK's zggevx] Dear Jim Demmel: The reported problem does involve the Exceptional shifts in QZ. Changing the QZ exceptional shift to something similar to what I have advixed for real QR and real QZ does fix this particular bug report. In my version of zhgeqz.f lines 603 read ELSE * * Exceptional shift. Chosen for no particularly good reason. * ESHIFT = ESHIFT + DCONJG( ( ASCALE*A( ILAST-1, ILAST ) ) / $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) ) After these lines, I added TEMP = $ ABS(ASCALE*A(ILAST-1,ILAST-2)/(BSCALE*B(ILAST-1,ILAST-1)))+ $ ABS(ASCALE*A(ILAST ,ILAST-1)/(BSCALE*B(ILAST-2,ILAST-2))) TEMP2 = $ ABS(ASCALE*A(ILAST-1,ILAST-2)/(BSCALE*B(ILAST-2,ILAST-2)))+ $ ABS(ASCALE*A(ILAST ,ILAST-1)/(BSCALE*B(ILAST-1,ILAST-1))) ESHIFT = DCMPLX(MIN(TEMP, TEMP2),0) This changes the value of the shift, now set on line 616 SHIFT = ESHIFT Thank you for including me in this discussion. --David M. Day [1.2.c-- Patrick Alken's patch] Date: Wed, 2 May 2007 16:14:55 -0600 From: Patrick Alken
Subject: Re: [Lapack] dggev fails on non-singular system Hello, I have done some testing and I have found that the following patch fixes all the 3x3 and 4x4 matrix systems that I found failures on: Replace (in dhgeqz.f): 629 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST) ).LT. 630 $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN 631 ESHIFT = ESHIFT + H( ILAST-1, ILAST ) / 632 $ T( ILAST-1, ILAST-1 ) 633 ELSE 634 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) ) 635 END IF with: 629 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1) ).LT. 630 $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN 631 ESHIFT = ESHIFT + 0.736 * H( ILAST, ILAST-1 ) / 632 $ T( ILAST-1, ILAST-1 ) 633 ELSE 634 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) ) 635 END IF The idea here is the following: in the previous code, its certainly possible for H(ILAST-1,ILAST) to equal 0, in which case the test will pass but ESHIFT will not be modified at all. I have just replaced it with H(ILAST,ILAST-1) which is guaranteed not to be zero, since its a subdiagonal element and the previous tests will search for zeros on the subdiagonal. The 0.736 coefficient is arbitrary - without it, I still found a small number of convergence failures on 4x4 integer systems. Patrick Alken[1.3-- Fix] Wait for a new QZ with early agressive deflation and multishift.
2 -- DQDS algorithm (useful to compute the singular values of a matrix)Cleve Moler reported that LAPACK DLASQ1 routine (called by DBDSQR when only the singular values are needed) fails on the bidiagonal matrix: set_Moler_200.c. Osni Marques provided us with DLASQ1-6 routines that passes the matrix. I have done some investigation to find out why the 200 by 200 matrix that Cleve Moler sent to us weeks ago leads to a failure in Matlab's "norm". My understanding is that "norm" resorts to LAPACK's dbdsqr, which in turn calls dqds, i.e. dlasq1 (which calls dlasq2 etc), if only singular vectors are desired. Indeed, LAPACK 3.1.1's dlasq1 returns INFO=2 (Matlab prints "Warning: SVD did not converge at index = 2") for Cleve's matrix, meaning "current block of Z not diagonalized after 30*N iterations (in inner while loop)". This is related to the NBIG = 30*( N0-I0+1 ) limit that is set in dlasq2, where I0 and N0 indicate the active (sub)matrix dlasq2 is dealing with. It starts with I0=1 and N0=N if there is no initial splitting; I0 and N0 may change later if splittings occur. These are my findings: 1) Contrary to what I had originally thought, LAPACK 3.1.1's dlasq2 does not contain our last modifications, in particular a mechanism to reverse (flip) the qd-array in some adverse situations. The lack of this mechanism may greatly delay convergence. In fact, LAPACK 3.1.1's dqds works if we allow more than 30*N iterations: without flipping, convergence is achieved in 8073 iterations (to achieve this I changed 30*N to 300*N). In contrast, the latest dqds that we have produced works just fine on Cleve's matrix: convergence is achieved in 969 iterations. Interestingly, LAPACK 3.1.1 converges in 970 iterations if we first flip Cleve's matrix (similarly, "norm" also works fine). 2) In the computation of shifts (by dlasq4), we have TTYPE.EQ.-6, i.e. "Case 6, no information to guide us." This case is based on a variable G that is used for the computation of the shift (S = G*DMIN). Modifications were performed in the dqds routines to make them thread safe (by removing SAVE statements) but in LAPACK 3.1.1 G seems to be set to zero too often. I am afraid that for a sequence of TTYPE.EQ.-6 this may lead to shifts that are not so efficient. My records indicate that I submitted our latest version of dqds for inclusion in LAPACK on 07/11/06 but the changes have not been incorporated (LAPACK 3.1 is dated 11/12/06, and LAPACK 3.1.1 is dated 02/26/07). In any case, I am sending it attached so it can be compared with LAPACK 3.1.1's dqds. I am also attaching Cleve's 200 by 200 matrix as well as its flipped form in case someone (in the cc list) would like to experiment with them. Cheers, Osni The reason for which Osni code is not in LAPACK is that Osni code fail on some of the matrices in the LAPACK testing, for example: set_PB1_30.c. Indeed the double precision routines of Osni are fine, there are only a few problems with the single precision routines. Current status: Hello lapackers, I did not understand upfront what Bobby Cheng was saying with his: "MATLAB switched back to the old algorithm and lost some speed." So I asked Bobby Cheng for their files. MathWorks simply have hacked xBDSQR so that they use the QR algorithm (they simply have commented the lines that call DLASQ1). So they use QR even if they want only the singular values. Matlab does not use the dqds algorithm when computing a norm for example. Julien.