LAPACK does not have any exceptional shift strategy for QZ

[langou]

LAPACK current shift strategy is known to fail on various pencils. Here are six known 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 Michael 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

See as well:
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4357

And probably also:
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4855

The fixes

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.

With best regards,
Daniel

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

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

[langou]

This issue was fixed in LAPACK 3.5. (This was SVN commit 1351 by Rodney James.)