dlaqr0()

subroutine dlaqr0	(	logical	wantt,
		logical	wantz,
		integer	n,
		integer	ilo,
		integer	ihi,
		double precision, dimension( ldh, * )	h,
		integer	ldh,
		double precision, dimension( * )	wr,
		double precision, dimension( * )	wi,
		integer	iloz,
		integer	ihiz,
		double precision, dimension( ldz, * )	z,
		integer	ldz,
		double precision, dimension( * )	work,
		integer	lwork,
		integer	info )

DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.

Download DLAQR0 + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!>    DLAQR0 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.
!> 

Parameters

[in]	WANTT	!> WANTT is LOGICAL !> = .TRUE. : the full Schur form T is required; !> = .FALSE.: only eigenvalues are required. !>
[in]	WANTZ	!> WANTZ is LOGICAL !> = .TRUE. : the matrix of Schur vectors Z is required; !> = .FALSE.: Schur vectors are not required. !>
[in]	N	!> N is INTEGER !> The order of the matrix H. N >= 0. !>
[in]	ILO	!> ILO is INTEGER !>
[in]	IHI	!> 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 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 > 0, then 1 <= ILO <= IHI <= N. !> If N = 0, then ILO = 1 and IHI = 0. !>
[in,out]	H	!> 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) < 0. If INFO = 0 and WANTT 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. !>
[in]	LDH	!> LDH is INTEGER !> The leading dimension of the array H. LDH >= max(1,N). !>
[out]	WR	!> WR is DOUBLE PRECISION array, dimension (IHI) !>
[out]	WI	!> 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) > 0 and WI(i+1) < 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). !>
[in]	ILOZ	!> ILOZ is INTEGER !>
[in]	IHIZ	!> IHIZ is INTEGER !> Specify the rows of Z to which transformations must be !> applied if WANTZ is .TRUE.. !> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. !>
[in,out]	Z	!> 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 > 0 is given under !> the description of INFO below.) !>
[in]	LDZ	!> LDZ is INTEGER !> The leading dimension of the array Z. if WANTZ is .TRUE. !> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1. !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension LWORK !> On exit, if LWORK = -1, WORK(1) returns an estimate of !> the optimal value for LWORK. !>
[in]	LWORK	!> 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 DLAQR0 does a workspace query. !> In this case, DLAQR0 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. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> > 0: if INFO = i, DLAQR0 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 an orthogonal matrix. The final !> value of H is upper Hessenberg and quasi-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 orthogonal matrix in (*) (regard- !> less of the value of WANTT.) !> !> If INFO > 0 and WANTZ is .FALSE., then Z is not !> accessed. !>

Contributors:

Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

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.

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.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 252 of file dlaqr0.f.

*
*  -- 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 ..
      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 = 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 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, dlaqr3, dlaqr4,
     $                   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 = 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, 'DLAQR0', 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, 'DLAQR0', 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 DLAQR3 ====
*
         CALL dlaqr3( 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, DLAQR3) ====
*
         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, 'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'DLAQR0', 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, 'DLAQR0', 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 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 dlaqr3( 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 DLAQR3
*              .    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
*              .    DLAQR3 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 DLAQR4 or
*                 .    DLAHQR 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 dlacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     IF( ns.GT.nmin ) THEN
                        CALL dlaqr4( .false., .false., ns, 1, ns,
     $                               h( kt, 1 ), ldh, wr( ks ),
     $                               wi( ks ), 1, 1, zdum, 1, work,
     $                               lwork, inf )
                     ELSE
                        CALL dlahqr( .false., .false., ns, 1, ns,
     $                               h( kt, 1 ), ldh, wr( ks ),
     $                               wi( ks ), 1, 1, zdum, 1, inf )
                     END IF
                     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 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 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 DLAQR0 ====
*

Here is the call graph for this function:

Here is the caller graph for this function: