slaqr0()

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

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

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

Purpose:

!>
!>    SLAQR0 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 SGEBAL, and then passed to SGEHRD when the !> matrix output by SGEBAL 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 REAL 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 REAL array, dimension (IHI) !>
[out]	WI	!> WI is REAL 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 REAL 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 REAL 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 SLAQR0 does a workspace query. !> In this case, SLAQR0 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, SLAQR0 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. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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.

Definition at line 252 of file slaqr0.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 ..
      REAL               H( LDH, * ), WI( * ), WORK( * ), WR( * ),
     $                   Z( LDZ, * )
*     ..
*
*  ================================================================
*     .. Parameters ..
*
*     ==== Matrices of order NTINY or smaller must be processed by
*     .    SLAHQR 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. ====
      REAL               WILK1, WILK2
      parameter( wilk1 = 0.75e0, wilk2 = -0.4375e0 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     ..
*     .. Local Scalars ..
      REAL               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 ..
      REAL               ZDUM( 1, 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           slacpy, slahqr, slanv2, slaqr3, slaqr4,
     $                   slaqr5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, int, max, min, mod, real
*     ..
*     .. 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 SLAHQR. ====
*
         lwkopt = 1
         IF( lwork.NE.-1 )
     $      CALL slahqr( 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, 'SLAQR0', 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, 'SLAQR0', 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 SLAQR3 ====
*
         CALL slaqr3( 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(SLAQR5, SLAQR3) ====
*
         lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
*
*        ==== Quick return in case of workspace query. ====
*
         IF( lwork.EQ.-1 ) THEN
            work( 1 ) = real( lwkopt )
            RETURN
         END IF
*
*        ==== SLAHQR/SLAQR0 crossover point ====
*
         nmin = ilaenv( 12, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'SLAQR0', 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, 'SLAQR0', 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 slaqr3( 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 SLAQR3
*              .    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
*              .    SLAQR3 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 slanv2( 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 SLAQR4 or
*                 .    SLAHQR 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 slacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     IF( ns.GT.nmin ) THEN
                        CALL slaqr4( .false., .false., ns, 1, ns,
     $                               h( kt, 1 ), ldh, wr( ks ),
     $                               wi( ks ), 1, 1, zdum, 1, work,
     $                               lwork, inf )
                     ELSE
                        CALL slahqr( .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 slanv2( 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 slaqr5( 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 ) = real( lwkopt )
*
*     ==== End of SLAQR0 ====
*

Here is the call graph for this function:

Here is the caller graph for this function: