claqr4()

subroutine claqr4	(	logical	wantt,
		logical	wantz,
		integer	n,
		integer	ilo,
		integer	ihi,
		complex, dimension( ldh, * )	h,
		integer	ldh,
		complex, dimension( * )	w,
		integer	iloz,
		integer	ihiz,
		complex, dimension( ldz, * )	z,
		integer	ldz,
		complex, dimension( * )	work,
		integer	lwork,
		integer	info
	)

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

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

Purpose:

    CLAQR4 implements one level of recursion for CLAQR0.
    It is a complete implementation of the small bulge multi-shift
    QR algorithm.  It may be called by CLAQR0 and, for large enough
    deflation window size, it may be called by CLAQR3.  This
    subroutine is identical to CLAQR0 except that it calls CLAQR2
    instead of CLAQR3.

    CLAQR4 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.

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 CGEBAL, and then passed to CGEHRD when the matrix output by CGEBAL 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 COMPLEX 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.
[in]	LDH	LDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).
[out]	W	W is COMPLEX 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).
[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 COMPLEX 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 COMPLEX 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 CLAQR4 does a workspace query. In this case, CLAQR4 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, CLAQR4 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.

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 246 of file claqr4.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 ..
      COMPLEX            H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
*     ..
*
*
*  ================================================================
*
*     .. Parameters ..
*
*     ==== Matrices of order NTINY or smaller must be processed by
*     .    CLAHQR 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. ====
      REAL               WILK1
      parameter( wilk1 = 0.75e0 )
      COMPLEX            ZERO, ONE
      parameter( zero = ( 0.0e0, 0.0e0 ),
     $                   one = ( 1.0e0, 0.0e0 ) )
      REAL               TWO
      parameter( two = 2.0e0 )
*     ..
*     .. Local Scalars ..
      COMPLEX            AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
      REAL               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            ZDUM( 1, 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           clacpy, clahqr, claqr2, claqr5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, cmplx, int, max, min, mod, real,
     $                   sqrt
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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 CLAHQR. ====
*
         lwkopt = 1
         IF( lwork.NE.-1 )
     $      CALL clahqr( 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, 'CLAQR4', 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, 'CLAQR4', 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 CLAQR2 ====
*
         CALL claqr2( 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(CLAQR5, CLAQR2) ====
*
         lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
*
*        ==== Quick return in case of workspace query. ====
*
         IF( lwork.EQ.-1 ) THEN
            work( 1 ) = cmplx( lwkopt, 0 )
            RETURN
         END IF
*
*        ==== CLAHQR/CLAQR0 crossover point ====
*
         nmin = ilaenv( 12, 'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'CLAQR4', 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, 'CLAQR4', 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 claqr2( 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 CLAQR2
*              .    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
*              .    CLAQR2 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 CLAHQR
*                 .    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 clacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     CALL clahqr( .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 claqr5( 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 ) = cmplx( lwkopt, 0 )
*
*     ==== End of CLAQR4 ====
*

Here is the call graph for this function:

Here is the caller graph for this function: