dlaqr4.f

Go to the documentation of this file.

*> \brief \b DLAQR4 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 DLAQR4 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr4.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr4.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr4.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
*                          ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO )
* 
*       .. Scalar Arguments ..
*       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
*       LOGICAL            WANTT, WANTZ
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
*      $                   Z( LDZ, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    DLAQR4 implements one level of recursion for DLAQR0.
*>    It is a complete implementation of the small bulge multi-shift
*>    QR algorithm.  It may be called by DLAQR0 and, for large enough
*>    deflation window size, it may be called by DLAQR3.  This
*>    subroutine is identical to DLAQR0 except that it calls DLAQR2
*>    instead of DLAQR3.
*>
*>    DLAQR4 computes the eigenvalues of a Hessenberg matrix H
*>    and, optionally, the matrices T and Z from the Schur decomposition
*>    H = Z T Z**T, where T is an upper quasi-triangular matrix (the
*>    Schur form), and Z is the orthogonal matrix of Schur vectors.
*>
*>    Optionally Z may be postmultiplied into an input orthogonal
*>    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 orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*> \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 .GE. 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.GT.1,
*>           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*>           previous call to DGEBAL, and then passed to DGEHRD when the
*>           matrix output by DGEBAL is reduced to Hessenberg form.
*>           Otherwise, ILO and IHI should be set to 1 and N,
*>           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
*>           If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim
*>
*> \param[in,out] H
*> \verbatim
*>          H is DOUBLE PRECISION 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 quasi-triangular matrix T from the Schur
*>           decomposition (the Schur form); 2-by-2 diagonal blocks
*>           (corresponding to complex conjugate pairs of eigenvalues)
*>           are returned in standard form, with H(i,i) = H(i+1,i+1)
*>           and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
*>           .FALSE., then the contents of H are unspecified on exit.
*>           (The output value of H when INFO.GT.0 is given under the
*>           description of INFO below.)
*>
*>           This subroutine may explicitly set H(i,j) = 0 for i.GT.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 .GE. max(1,N).
*> \endverbatim
*>
*> \param[out] WR
*> \verbatim
*>          WR is DOUBLE PRECISION array, dimension (IHI)
*> \endverbatim
*>
*> \param[out] WI
*> \verbatim
*>          WI is DOUBLE PRECISION array, dimension (IHI)
*>           The real and imaginary parts, respectively, of the computed
*>           eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
*>           and WI(ILO:IHI). If two eigenvalues are computed as a
*>           complex conjugate pair, they are stored in consecutive
*>           elements of WR and WI, say the i-th and (i+1)th, with
*>           WI(i) .GT. 0 and WI(i+1) .LT. 0. 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
*>           WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
*>           block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
*>           WI(i+1) = -WI(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 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION 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.GT.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.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION 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 .GE. 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 DLAQR4 does a workspace query.
*>           In this case, DLAQR4 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
*>           .GT. 0:  if INFO = i, DLAQR4 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 .GT. 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 .GT. 0 and WANTT is .TRUE., then on exit
*>
*>           (*)  (initial value of H)*U  = U*(final value of H)
*>
*>                where U is a orthogonal matrix.  The final
*>                value of  H is upper Hessenberg and triangular in
*>                rows and columns INFO+1 through IHI.
*>
*>                If INFO .GT. 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 orthogonal matrix in (*) (regard-
*>                less of the value of WANTT.)
*>
*>                If INFO .GT. 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. 
*
*> \date September 2012
*
*> \ingroup doubleOTHERauxiliary
*
*> \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 dlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
     $                   iloz, ihiz, z, ldz, work, lwork, info )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
      LOGICAL            WANTT, WANTZ
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   H( ldh, * ), WI( * ), WORK( * ), WR( * ),
     $                   z( ldz, * )
*     ..
*
*  ================================================================
*     .. Parameters ..
*
*     ==== Matrices of order NTINY or smaller must be processed by
*     .    DLAHQR because of insufficient subdiagonal scratch space.
*     .    (This is a hard limit.) ====
      INTEGER            NTINY
      parameter                ( ntiny = 11 )
*
*     ==== 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 constants WILK1 and WILK2 are used to form the
*     .    exceptional shifts. ====
      DOUBLE PRECISION   WILK1, WILK2
      parameter                ( wilk1 = 0.75d0, wilk2 = -0.4375d0 )
      DOUBLE PRECISION   ZERO, ONE
      parameter                ( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   AA, BB, CC, CS, DD, SN, SS, SWAP
      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 ..
      DOUBLE PRECISION   ZDUM( 1, 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, int, max, min, mod
*     ..
*     .. 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 DLAHQR. ====
*
         lwkopt = 1
         IF( lwork.NE.-1 )
     $      CALL dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
     $                   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 = 11, so there is enough
*        .    subdiagonal workspace for NWR.GE.2 as required.
*        .    (In fact, there is enough subdiagonal space for
*        .    NWR.GE.3.) ====
*
         nwr = ilaenv( 13, 'DLAQR4', 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 = 11, so there is at
*        .    enough subdiagonal workspace for NSR to be even
*        .    and greater than or equal to two as required. ====
*
         nsr = ilaenv( 15, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
         nsr = max( 2, nsr-mod( nsr, 2 ) )
*
*        ==== Estimate optimal workspace ====
*
*        ==== Workspace query call to DLAQR2 ====
*
         CALL dlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
     $                ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
     $                n, h, ldh, work, -1 )
*
*        ==== Optimal workspace = MAX(DLAQR5, DLAQR2) ====
*
         lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
*
*        ==== Quick return in case of workspace query. ====
*
         IF( lwork.EQ.-1 ) THEN
            work( 1 ) = dble( lwkopt )
            RETURN
         END IF
*
*        ==== DLAHQR/DLAQR0 crossover point ====
*
         nmin = ilaenv( 12, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'DLAQR4', 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, 'DLAQR4', 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+6 ) / 9, 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 80 it = 1, itmax
*
*           ==== Done when KBOT falls below ILO ====
*
            IF( kbot.LT.ilo )
     $         GO TO 90
*
*           ==== 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( abs( h( kwtop, kwtop-1 ) ).GT.
     $                abs( 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 dlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
     $                   ihiz, z, ldz, ls, ld, wr, wi, 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 DLAQR2
*              .    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
*              .    DLAQR2 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, max( ks+1, ktop+2 ), -2
                     ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
                     aa = wilk1*ss + h( i, i )
                     bb = ss
                     cc = wilk2*ss
                     dd = aa
                     CALL dlanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
     $                            wr( i ), wi( i ), cs, sn )
   30             CONTINUE
                  IF( ks.EQ.ktop ) THEN
                     wr( ks+1 ) = h( ks+1, ks+1 )
                     wi( ks+1 ) = zero
                     wr( ks ) = wr( ks+1 )
                     wi( ks ) = wi( ks+1 )
                  END IF
               ELSE
*
*                 ==== Got NS/2 or fewer shifts? Use DLAHQR
*                 .    on a trailing principal submatrix to
*                 .    get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
*                 .    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 dlacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     CALL dlahqr( .false., .false., ns, 1, ns,
     $                            h( kt, 1 ), ldh, wr( ks ), wi( 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.  ====
*
                     IF( ks.GE.kbot ) THEN
                        aa = h( kbot-1, kbot-1 )
                        cc = h( kbot, kbot-1 )
                        bb = h( kbot-1, kbot )
                        dd = h( kbot, kbot )
                        CALL dlanv2( aa, bb, cc, dd, wr( kbot-1 ),
     $                               wi( kbot-1 ), wr( kbot ),
     $                               wi( kbot ), cs, sn )
                        ks = kbot - 1
                     END IF
                  END IF
*
                  IF( kbot-ks+1.GT.ns ) THEN
*
*                    ==== Sort the shifts (Helps a little)
*                    .    Bubble sort keeps complex conjugate
*                    .    pairs together. ====
*
                     sorted = .false.
                     DO 50 k = kbot, ks + 1, -1
                        IF( sorted )
     $                     GO TO 60
                        sorted = .true.
                        DO 40 i = ks, k - 1
                           IF( abs( wr( i ) )+abs( wi( i ) ).LT.
     $                         abs( wr( i+1 ) )+abs( wi( i+1 ) ) ) THEN
                              sorted = .false.
*
                              swap = wr( i )
                              wr( i ) = wr( i+1 )
                              wr( i+1 ) = swap
*
                              swap = wi( i )
                              wi( i ) = wi( i+1 )
                              wi( i+1 ) = swap
                           END IF
   40                   CONTINUE
   50                CONTINUE
   60                CONTINUE
                  END IF
*
*                 ==== Shuffle shifts into pairs of real shifts
*                 .    and pairs of complex conjugate shifts
*                 .    assuming complex conjugate shifts are
*                 .    already adjacent to one another. (Yes,
*                 .    they are.)  ====
*
                  DO 70 i = kbot, ks + 2, -2
                     IF( wi( i ).NE.-wi( i-1 ) ) THEN
*
                        swap = wr( i )
                        wr( i ) = wr( i-1 )
                        wr( i-1 ) = wr( i-2 )
                        wr( i-2 ) = swap
*
                        swap = wi( i )
                        wi( i ) = wi( i-1 )
                        wi( i-1 ) = wi( i-2 )
                        wi( i-2 ) = swap
                     END IF
   70             CONTINUE
               END IF
*
*              ==== If there are only two shifts and both are
*              .    real, then use only one.  ====
*
               IF( kbot-ks+1.EQ.2 ) THEN
                  IF( wi( kbot ).EQ.zero ) THEN
                     IF( abs( wr( kbot )-h( kbot, kbot ) ).LT.
     $                   abs( wr( kbot-1 )-h( kbot, kbot ) ) ) THEN
                        wr( kbot-1 ) = wr( kbot )
                     ELSE
                        wr( kbot ) = wr( kbot-1 )
                     END IF
                  END IF
               END IF
*
*              ==== Use up to NS of the the smallest magnatiude
*              .    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 = 3*ns - 3
               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 dlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
     $                      wr( ks ), wi( 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 ====
   80    CONTINUE
*
*        ==== Iteration limit exceeded.  Set INFO to show where
*        .    the problem occurred and exit. ====
*
         info = kbot
   90    CONTINUE
      END IF
*
*     ==== Return the optimal value of LWORK. ====
*
      work( 1 ) = dble( lwkopt )
*
*     ==== End of DLAQR4 ====
*
      END