zlaqr4.f

*> \brief \b ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZLAQR4 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr4.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr4.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr4.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
*                          IHIZ, Z, LDZ, WORK, LWORK, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
*       LOGICAL            WANTT, WANTZ
*       ..
*       .. Array Arguments ..
*       COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    ZLAQR4 implements one level of recursion for ZLAQR0.
*>    It is a complete implementation of the small bulge multi-shift
*>    QR algorithm.  It may be called by ZLAQR0 and, for large enough
*>    deflation window size, it may be called by ZLAQR3.  This
*>    subroutine is identical to ZLAQR0 except that it calls ZLAQR2
*>    instead of ZLAQR3.
*>
*>    ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
*>    and, optionally, the matrices T and Z from the Schur decomposition
*>    H = Z T Z**H, where T is an upper triangular matrix (the
*>    Schur form), and Z is the unitary matrix of Schur vectors.
*>
*>    Optionally Z may be postmultiplied into an input unitary
*>    matrix Q so that this routine can give the Schur factorization
*>    of a matrix A which has been reduced to the Hessenberg form H
*>    by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] WANTT
*> \verbatim
*>          WANTT is LOGICAL
*>          = .TRUE. : the full Schur form T is required;
*>          = .FALSE.: only eigenvalues are required.
*> \endverbatim
*>
*> \param[in] WANTZ
*> \verbatim
*>          WANTZ is LOGICAL
*>          = .TRUE. : the matrix of Schur vectors Z is required;
*>          = .FALSE.: Schur vectors are not required.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>           The order of the matrix H.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*>           It is assumed that H is already upper triangular in rows
*>           and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*>           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*>           previous call to ZGEBAL, and then passed to ZGEHRD when the
*>           matrix output by ZGEBAL is reduced to Hessenberg form.
*>           Otherwise, ILO and IHI should be set to 1 and N,
*>           respectively.  If N > 0, then 1 <= ILO <= IHI <= N.
*>           If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim
*>
*> \param[in,out] H
*> \verbatim
*>          H is COMPLEX*16 array, dimension (LDH,N)
*>           On entry, the upper Hessenberg matrix H.
*>           On exit, if INFO = 0 and WANTT is .TRUE., then H
*>           contains the upper triangular matrix T from the Schur
*>           decomposition (the Schur form). If INFO = 0 and WANT is
*>           .FALSE., then the contents of H are unspecified on exit.
*>           (The output value of H when INFO > 0 is given under the
*>           description of INFO below.)
*>
*>           This subroutine may explicitly set H(i,j) = 0 for i > j and
*>           j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*>          LDH is INTEGER
*>           The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim
*>
*> \param[out] W
*> \verbatim
*>          W is COMPLEX*16 array, dimension (N)
*>           The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
*>           in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
*>           stored in the same order as on the diagonal of the Schur
*>           form returned in H, with W(i) = H(i,i).
*> \endverbatim
*>
*> \param[in] ILOZ
*> \verbatim
*>          ILOZ is INTEGER
*> \endverbatim
*>
*> \param[in] IHIZ
*> \verbatim
*>          IHIZ is INTEGER
*>           Specify the rows of Z to which transformations must be
*>           applied if WANTZ is .TRUE..
*>           1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is COMPLEX*16 array, dimension (LDZ,IHI)
*>           If WANTZ is .FALSE., then Z is not referenced.
*>           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*>           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*>           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
*>           (The output value of Z when INFO > 0 is given under
*>           the description of INFO below.)
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>           The leading dimension of the array Z.  if WANTZ is .TRUE.
*>           then LDZ >= MAX(1,IHIZ).  Otherwise, LDZ >= 1.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension LWORK
*>           On exit, if LWORK = -1, WORK(1) returns an estimate of
*>           the optimal value for LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>           The dimension of the array WORK.  LWORK >= max(1,N)
*>           is sufficient, but LWORK typically as large as 6*N may
*>           be required for optimal performance.  A workspace query
*>           to determine the optimal workspace size is recommended.
*>
*>           If LWORK = -1, then ZLAQR4 does a workspace query.
*>           In this case, ZLAQR4 checks the input parameters and
*>           estimates the optimal workspace size for the given
*>           values of N, ILO and IHI.  The estimate is returned
*>           in WORK(1).  No error message related to LWORK is
*>           issued by XERBLA.  Neither H nor Z are accessed.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>             =  0:  successful exit
*>             > 0:  if INFO = i, ZLAQR4 failed to compute all of
*>                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
*>                and WI contain those eigenvalues which have been
*>                successfully computed.  (Failures are rare.)
*>
*>                If INFO > 0 and WANT is .FALSE., then on exit,
*>                the remaining unconverged eigenvalues are the eigen-
*>                values of the upper Hessenberg matrix rows and
*>                columns ILO through INFO of the final, output
*>                value of H.
*>
*>                If INFO > 0 and WANTT is .TRUE., then on exit
*>
*>           (*)  (initial value of H)*U  = U*(final value of H)
*>
*>                where U is a unitary matrix.  The final
*>                value of  H is upper Hessenberg and triangular in
*>                rows and columns INFO+1 through IHI.
*>
*>                If INFO > 0 and WANTZ is .TRUE., then on exit
*>
*>                  (final value of Z(ILO:IHI,ILOZ:IHIZ)
*>                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
*>
*>                where U is the unitary matrix in (*) (regard-
*>                less of the value of WANTT.)
*>
*>                If INFO > 0 and WANTZ is .FALSE., then Z is not
*>                accessed.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup laqr4
*
*> \par Contributors:
*  ==================
*>
*>       Karen Braman and Ralph Byers, Department of Mathematics,
*>       University of Kansas, USA
*
*> \par References:
*  ================
*>
*>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
*>       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
*>       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
*>       929--947, 2002.
*> \n
*>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
*>       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
*>       of Matrix Analysis, volume 23, pages 948--973, 2002.
*>
*  =====================================================================
      SUBROUTINE zlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
     $                   IHIZ, Z, LDZ, WORK, LWORK, INFO )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
      LOGICAL            WANTT, WANTZ
*     ..
*     .. Array Arguments ..
      COMPLEX*16         H( ldh, * ), W( * ), WORK( * ), Z( ldz, * )
*     ..
*
*  ================================================================
*
*     .. Parameters ..
*
*     ==== Matrices of order NTINY or smaller must be processed by
*     .    ZLAHQR because of insufficient subdiagonal scratch space.
*     .    (This is a hard limit.) ====
      INTEGER            NTINY
      parameter( ntiny = 15 )
*
*     ==== Exceptional deflation windows:  try to cure rare
*     .    slow convergence by varying the size of the
*     .    deflation window after KEXNW iterations. ====
      INTEGER            KEXNW
      parameter( kexnw = 5 )
*
*     ==== Exceptional shifts: try to cure rare slow convergence
*     .    with ad-hoc exceptional shifts every KEXSH iterations.
*     .    ====
      INTEGER            KEXSH
      parameter( kexsh = 6 )
*
*     ==== The constant WILK1 is used to form the exceptional
*     .    shifts. ====
      DOUBLE PRECISION   WILK1
      parameter( wilk1 = 0.75d0 )
      COMPLEX*16         ZERO, ONE
      parameter( zero = ( 0.0d0, 0.0d0 ),
     $                   one = ( 1.0d0, 0.0d0 ) )
      DOUBLE PRECISION   TWO
      parameter( two = 2.0d0 )
*     ..
*     .. Local Scalars ..
      COMPLEX*16         AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
      DOUBLE PRECISION   S
      INTEGER            I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
     $                   kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
     $                   lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
     $                   nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
      LOGICAL            SORTED
      CHARACTER          JBCMPZ*2
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      EXTERNAL           ilaenv
*     ..
*     .. Local Arrays ..
      COMPLEX*16         ZDUM( 1, 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           zlacpy, zlahqr, zlaqr2, zlaqr5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dcmplx, dimag, int, max, min, mod,
     $                   sqrt
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
*     ..
*     .. Executable Statements ..
      info = 0
*
*     ==== Quick return for N = 0: nothing to do. ====
*
      IF( n.EQ.0 ) THEN
         work( 1 ) = one
         RETURN
      END IF
*
      IF( n.LE.ntiny ) THEN
*
*        ==== Tiny matrices must use ZLAHQR. ====
*
         lwkopt = 1
         IF( lwork.NE.-1 )
     $      CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
     $                   ihiz, z, ldz, info )
      ELSE
*
*        ==== Use small bulge multi-shift QR with aggressive early
*        .    deflation on larger-than-tiny matrices. ====
*
*        ==== Hope for the best. ====
*
         info = 0
*
*        ==== Set up job flags for ILAENV. ====
*
         IF( wantt ) THEN
            jbcmpz( 1: 1 ) = 'S'
         ELSE
            jbcmpz( 1: 1 ) = 'E'
         END IF
         IF( wantz ) THEN
            jbcmpz( 2: 2 ) = 'V'
         ELSE
            jbcmpz( 2: 2 ) = 'N'
         END IF
*
*        ==== NWR = recommended deflation window size.  At this
*        .    point,  N .GT. NTINY = 15, so there is enough
*        .    subdiagonal workspace for NWR.GE.2 as required.
*        .    (In fact, there is enough subdiagonal space for
*        .    NWR.GE.4.) ====
*
         nwr = ilaenv( 13, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nwr = max( 2, nwr )
         nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
*
*        ==== NSR = recommended number of simultaneous shifts.
*        .    At this point N .GT. NTINY = 15, so there is at
*        .    enough subdiagonal workspace for NSR to be even
*        .    and greater than or equal to two as required. ====
*
         nsr = ilaenv( 15, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
         nsr = max( 2, nsr-mod( nsr, 2 ) )
*
*        ==== Estimate optimal workspace ====
*
*        ==== Workspace query call to ZLAQR2 ====
*
         CALL zlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
     $                ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
     $                ldh, work, -1 )
*
*        ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ====
*
         lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
*
*        ==== Quick return in case of workspace query. ====
*
         IF( lwork.EQ.-1 ) THEN
            work( 1 ) = dcmplx( lwkopt, 0 )
            RETURN
         END IF
*
*        ==== ZLAHQR/ZLAQR0 crossover point ====
*
         nmin = ilaenv( 12, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nibble = max( 0, nibble )
*
*        ==== Accumulate reflections during ttswp?  Use block
*        .    2-by-2 structure during matrix-matrix multiply? ====
*
         kacc22 = ilaenv( 16, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         kacc22 = max( 0, kacc22 )
         kacc22 = min( 2, kacc22 )
*
*        ==== NWMAX = the largest possible deflation window for
*        .    which there is sufficient workspace. ====
*
         nwmax = min( ( n-1 ) / 3, lwork / 2 )
         nw = nwmax
*
*        ==== NSMAX = the Largest number of simultaneous shifts
*        .    for which there is sufficient workspace. ====
*
         nsmax = min( ( n-3 ) / 6, 2*lwork / 3 )
         nsmax = nsmax - mod( nsmax, 2 )
*
*        ==== NDFL: an iteration count restarted at deflation. ====
*
         ndfl = 1
*
*        ==== ITMAX = iteration limit ====
*
         itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
*
*        ==== Last row and column in the active block ====
*
         kbot = ihi
*
*        ==== Main Loop ====
*
         DO 70 it = 1, itmax
*
*           ==== Done when KBOT falls below ILO ====
*
            IF( kbot.LT.ilo )
     $         GO TO 80
*
*           ==== Locate active block ====
*
            DO 10 k = kbot, ilo + 1, -1
               IF( h( k, k-1 ).EQ.zero )
     $            GO TO 20
   10       CONTINUE
            k = ilo
   20       CONTINUE
            ktop = k
*
*           ==== Select deflation window size:
*           .    Typical Case:
*           .      If possible and advisable, nibble the entire
*           .      active block.  If not, use size MIN(NWR,NWMAX)
*           .      or MIN(NWR+1,NWMAX) depending upon which has
*           .      the smaller corresponding subdiagonal entry
*           .      (a heuristic).
*           .
*           .    Exceptional Case:
*           .      If there have been no deflations in KEXNW or
*           .      more iterations, then vary the deflation window
*           .      size.   At first, because, larger windows are,
*           .      in general, more powerful than smaller ones,
*           .      rapidly increase the window to the maximum possible.
*           .      Then, gradually reduce the window size. ====
*
            nh = kbot - ktop + 1
            nwupbd = min( nh, nwmax )
            IF( ndfl.LT.kexnw ) THEN
               nw = min( nwupbd, nwr )
            ELSE
               nw = min( nwupbd, 2*nw )
            END IF
            IF( nw.LT.nwmax ) THEN
               IF( nw.GE.nh-1 ) THEN
                  nw = nh
               ELSE
                  kwtop = kbot - nw + 1
                  IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
     $                cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
               END IF
            END IF
            IF( ndfl.LT.kexnw ) THEN
               ndec = -1
            ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
               ndec = ndec + 1
               IF( nw-ndec.LT.2 )
     $            ndec = 0
               nw = nw - ndec
            END IF
*
*           ==== Aggressive early deflation:
*           .    split workspace under the subdiagonal into
*           .      - an nw-by-nw work array V in the lower
*           .        left-hand-corner,
*           .      - an NW-by-at-least-NW-but-more-is-better
*           .        (NW-by-NHO) horizontal work array along
*           .        the bottom edge,
*           .      - an at-least-NW-but-more-is-better (NHV-by-NW)
*           .        vertical work array along the left-hand-edge.
*           .        ====
*
            kv = n - nw + 1
            kt = nw + 1
            nho = ( n-nw-1 ) - kt + 1
            kwv = nw + 2
            nve = ( n-nw ) - kwv + 1
*
*           ==== Aggressive early deflation ====
*
            CALL zlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
     $                   ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
     $                   h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
     $                   lwork )
*
*           ==== Adjust KBOT accounting for new deflations. ====
*
            kbot = kbot - ld
*
*           ==== KS points to the shifts. ====
*
            ks = kbot - ls + 1
*
*           ==== Skip an expensive QR sweep if there is a (partly
*           .    heuristic) reason to expect that many eigenvalues
*           .    will deflate without it.  Here, the QR sweep is
*           .    skipped if many eigenvalues have just been deflated
*           .    or if the remaining active block is small.
*
            IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
     $          ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
*
*              ==== NS = nominal number of simultaneous shifts.
*              .    This may be lowered (slightly) if ZLAQR2
*              .    did not provide that many shifts. ====
*
               ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
               ns = ns - mod( ns, 2 )
*
*              ==== If there have been no deflations
*              .    in a multiple of KEXSH iterations,
*              .    then try exceptional shifts.
*              .    Otherwise use shifts provided by
*              .    ZLAQR2 above or from the eigenvalues
*              .    of a trailing principal submatrix. ====
*
               IF( mod( ndfl, kexsh ).EQ.0 ) THEN
                  ks = kbot - ns + 1
                  DO 30 i = kbot, ks + 1, -2
                     w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
                     w( i-1 ) = w( i )
   30             CONTINUE
               ELSE
*
*                 ==== Got NS/2 or fewer shifts? Use ZLAHQR
*                 .    on a trailing principal submatrix to
*                 .    get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
*                 .    there is enough space below the subdiagonal
*                 .    to fit an NS-by-NS scratch array.) ====
*
                  IF( kbot-ks+1.LE.ns / 2 ) THEN
                     ks = kbot - ns + 1
                     kt = n - ns + 1
                     CALL zlacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     CALL zlahqr( .false., .false., ns, 1, ns,
     $                            h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
     $                            1, inf )
                     ks = ks + inf
*
*                    ==== In case of a rare QR failure use
*                    .    eigenvalues of the trailing 2-by-2
*                    .    principal submatrix.  Scale to avoid
*                    .    overflows, underflows and subnormals.
*                    .    (The scale factor S can not be zero,
*                    .    because H(KBOT,KBOT-1) is nonzero.) ====
*
                     IF( ks.GE.kbot ) THEN
                        s = cabs1( h( kbot-1, kbot-1 ) ) +
     $                      cabs1( h( kbot, kbot-1 ) ) +
     $                      cabs1( h( kbot-1, kbot ) ) +
     $                      cabs1( h( kbot, kbot ) )
                        aa = h( kbot-1, kbot-1 ) / s
                        cc = h( kbot, kbot-1 ) / s
                        bb = h( kbot-1, kbot ) / s
                        dd = h( kbot, kbot ) / s
                        tr2 = ( aa+dd ) / two
                        det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
                        rtdisc = sqrt( -det )
                        w( kbot-1 ) = ( tr2+rtdisc )*s
                        w( kbot ) = ( tr2-rtdisc )*s
*
                        ks = kbot - 1
                     END IF
                  END IF
*
                  IF( kbot-ks+1.GT.ns ) THEN
*
*                    ==== Sort the shifts (Helps a little) ====
*
                     sorted = .false.
                     DO 50 k = kbot, ks + 1, -1
                        IF( sorted )
     $                     GO TO 60
                        sorted = .true.
                        DO 40 i = ks, k - 1
                           IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
     $                          THEN
                              sorted = .false.
                              swap = w( i )
                              w( i ) = w( i+1 )
                              w( i+1 ) = swap
                           END IF
   40                   CONTINUE
   50                CONTINUE
   60                CONTINUE
                  END IF
               END IF
*
*              ==== If there are only two shifts, then use
*              .    only one.  ====
*
               IF( kbot-ks+1.EQ.2 ) THEN
                  IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
     $                cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
                     w( kbot-1 ) = w( kbot )
                  ELSE
                     w( kbot ) = w( kbot-1 )
                  END IF
               END IF
*
*              ==== Use up to NS of the the smallest magnitude
*              .    shifts.  If there aren't NS shifts available,
*              .    then use them all, possibly dropping one to
*              .    make the number of shifts even. ====
*
               ns = min( ns, kbot-ks+1 )
               ns = ns - mod( ns, 2 )
               ks = kbot - ns + 1
*
*              ==== Small-bulge multi-shift QR sweep:
*              .    split workspace under the subdiagonal into
*              .    - a KDU-by-KDU work array U in the lower
*              .      left-hand-corner,
*              .    - a KDU-by-at-least-KDU-but-more-is-better
*              .      (KDU-by-NHo) horizontal work array WH along
*              .      the bottom edge,
*              .    - and an at-least-KDU-but-more-is-better-by-KDU
*              .      (NVE-by-KDU) vertical work WV arrow along
*              .      the left-hand-edge. ====
*
               kdu = 2*ns
               ku = n - kdu + 1
               kwh = kdu + 1
               nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
               kwv = kdu + 4
               nve = n - kdu - kwv + 1
*
*              ==== Small-bulge multi-shift QR sweep ====
*
               CALL zlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
     $                      w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
     $                      3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
     $                      nho, h( ku, kwh ), ldh )
            END IF
*
*           ==== Note progress (or the lack of it). ====
*
            IF( ld.GT.0 ) THEN
               ndfl = 1
            ELSE
               ndfl = ndfl + 1
            END IF
*
*           ==== End of main loop ====
   70    CONTINUE
*
*        ==== Iteration limit exceeded.  Set INFO to show where
*        .    the problem occurred and exit. ====
*
         info = kbot
   80    CONTINUE
      END IF
*
*     ==== Return the optimal value of LWORK. ====
*
      work( 1 ) = dcmplx( lwkopt, 0 )
*
*     ==== End of ZLAQR4 ====
*
      END