zlaqr4()

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

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

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

Purpose:

!>
!>    ZLAQR4 implements one level of recursion for ZLAQR0.
!>    It is a complete implementation of the small bulge multi-shift
!>    QR algorithm.  It may be called by ZLAQR0 and, for large enough
!>    deflation window size, it may be called by ZLAQR3.  This
!>    subroutine is identical to ZLAQR0 except that it calls ZLAQR2
!>    instead of ZLAQR3.
!>
!>    ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
!>    and, optionally, the matrices T and Z from the Schur decomposition
!>    H = Z T Z**H, where T is an upper triangular matrix (the
!>    Schur form), and Z is the unitary matrix of Schur vectors.
!>
!>    Optionally Z may be postmultiplied into an input unitary
!>    matrix Q so that this routine can give the Schur factorization
!>    of a matrix A which has been reduced to the Hessenberg form H
!>    by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.
!> 

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 ZGEBAL, and then passed to ZGEHRD when the !> matrix output by ZGEBAL is reduced to Hessenberg form. !> Otherwise, ILO and IHI should be set to 1 and N, !> respectively. If N > 0, then 1 <= ILO <= IHI <= N. !> If N = 0, then ILO = 1 and IHI = 0. !>
[in,out]	H	!> H is COMPLEX*16 array, dimension (LDH,N) !> On entry, the upper Hessenberg matrix H. !> On exit, if INFO = 0 and WANTT is .TRUE., then H !> contains the upper triangular matrix T from the Schur !> decomposition (the Schur form). If INFO = 0 and WANT is !> .FALSE., then the contents of H are unspecified on exit. !> (The output value of H when INFO > 0 is given under the !> description of INFO below.) !> !> This subroutine may explicitly set H(i,j) = 0 for i > j and !> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. !>
[in]	LDH	!> LDH is INTEGER !> The leading dimension of the array H. LDH >= max(1,N). !>
[out]	W	!> W is COMPLEX*16 array, dimension (N) !> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored !> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are !> stored in the same order as on the diagonal of the Schur !> form returned in H, with W(i) = H(i,i). !>
[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*16 array, dimension (LDZ,IHI) !> If WANTZ is .FALSE., then Z is not referenced. !> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is !> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the !> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). !> (The output value of Z when INFO > 0 is given under !> the description of INFO below.) !>
[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*16 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 ZLAQR4 does a workspace query. !> In this case, ZLAQR4 checks the input parameters and !> estimates the optimal workspace size for the given !> values of N, ILO and IHI. The estimate is returned !> in WORK(1). No error message related to LWORK is !> issued by XERBLA. Neither H nor Z are accessed. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> > 0: if INFO = i, ZLAQR4 failed to compute all of !> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR !> and WI contain those eigenvalues which have been !> successfully computed. (Failures are rare.) !> !> If INFO > 0 and WANT is .FALSE., then on exit, !> the remaining unconverged eigenvalues are the eigen- !> values of the upper Hessenberg matrix rows and !> columns ILO through INFO of the final, output !> value of H. !> !> If INFO > 0 and WANTT is .TRUE., then on exit !> !> (*) (initial value of H)*U = U*(final value of H) !> !> where U is a unitary matrix. The final !> value of H is upper Hessenberg and triangular in !> rows and columns INFO+1 through IHI. !> !> If INFO > 0 and WANTZ is .TRUE., then on exit !> !> (final value of Z(ILO:IHI,ILOZ:IHIZ) !> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U !> !> where U is the unitary matrix in (*) (regard- !> less of the value of WANTT.) !> !> If INFO > 0 and WANTZ is .FALSE., then Z is not !> accessed. !>

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 243 of file zlaqr4.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*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
*     ..
*
*  ================================================================
*
*     .. Parameters ..
*
*     ==== Matrices of order NTINY or smaller must be processed by
*     .    ZLAHQR because of insufficient subdiagonal scratch space.
*     .    (This is a hard limit.) ====
      INTEGER            NTINY
      parameter( ntiny = 15 )
*
*     ==== Exceptional deflation windows:  try to cure rare
*     .    slow convergence by varying the size of the
*     .    deflation window after KEXNW iterations. ====
      INTEGER            KEXNW
      parameter( kexnw = 5 )
*
*     ==== Exceptional shifts: try to cure rare slow convergence
*     .    with ad-hoc exceptional shifts every KEXSH iterations.
*     .    ====
      INTEGER            KEXSH
      parameter( kexsh = 6 )
*
*     ==== The constant WILK1 is used to form the exceptional
*     .    shifts. ====
      DOUBLE PRECISION   WILK1
      parameter( wilk1 = 0.75d0 )
      COMPLEX*16         ZERO, ONE
      parameter( zero = ( 0.0d0, 0.0d0 ),
     $                   one = ( 1.0d0, 0.0d0 ) )
      DOUBLE PRECISION   TWO
      parameter( two = 2.0d0 )
*     ..
*     .. Local Scalars ..
      COMPLEX*16         AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
      DOUBLE PRECISION   S
      INTEGER            I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
     $                   KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS,
     $                   LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS,
     $                   NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD
      LOGICAL            SORTED
      CHARACTER          JBCMPZ*2
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      EXTERNAL           ilaenv
*     ..
*     .. Local Arrays ..
      COMPLEX*16         ZDUM( 1, 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           zlacpy, zlahqr, zlaqr2, zlaqr5
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dcmplx, dimag, int, max, min, mod,
     $                   sqrt
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
*     ..
*     .. Executable Statements ..
      info = 0
*
*     ==== Quick return for N = 0: nothing to do. ====
*
      IF( n.EQ.0 ) THEN
         work( 1 ) = one
         RETURN
      END IF
*
      IF( n.LE.ntiny ) THEN
*
*        ==== Tiny matrices must use ZLAHQR. ====
*
         lwkopt = 1
         IF( lwork.NE.-1 )
     $      CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
     $                   ihiz, z, ldz, info )
      ELSE
*
*        ==== Use small bulge multi-shift QR with aggressive early
*        .    deflation on larger-than-tiny matrices. ====
*
*        ==== Hope for the best. ====
*
         info = 0
*
*        ==== Set up job flags for ILAENV. ====
*
         IF( wantt ) THEN
            jbcmpz( 1: 1 ) = 'S'
         ELSE
            jbcmpz( 1: 1 ) = 'E'
         END IF
         IF( wantz ) THEN
            jbcmpz( 2: 2 ) = 'V'
         ELSE
            jbcmpz( 2: 2 ) = 'N'
         END IF
*
*        ==== NWR = recommended deflation window size.  At this
*        .    point,  N .GT. NTINY = 15, so there is enough
*        .    subdiagonal workspace for NWR.GE.2 as required.
*        .    (In fact, there is enough subdiagonal space for
*        .    NWR.GE.4.) ====
*
         nwr = ilaenv( 13, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nwr = max( 2, nwr )
         nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
*
*        ==== NSR = recommended number of simultaneous shifts.
*        .    At this point N .GT. NTINY = 15, so there is at
*        .    enough subdiagonal workspace for NSR to be even
*        .    and greater than or equal to two as required. ====
*
         nsr = ilaenv( 15, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
         nsr = max( 2, nsr-mod( nsr, 2 ) )
*
*        ==== Estimate optimal workspace ====
*
*        ==== Workspace query call to ZLAQR2 ====
*
         CALL zlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
     $                ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
     $                ldh, work, -1 )
*
*        ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ====
*
         lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
*
*        ==== Quick return in case of workspace query. ====
*
         IF( lwork.EQ.-1 ) THEN
            work( 1 ) = dcmplx( lwkopt, 0 )
            RETURN
         END IF
*
*        ==== ZLAHQR/ZLAQR0 crossover point ====
*
         nmin = ilaenv( 12, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        ==== Nibble crossover point ====
*
         nibble = ilaenv( 14, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         nibble = max( 0, nibble )
*
*        ==== Accumulate reflections during ttswp?  Use block
*        .    2-by-2 structure during matrix-matrix multiply? ====
*
         kacc22 = ilaenv( 16, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
         kacc22 = max( 0, kacc22 )
         kacc22 = min( 2, kacc22 )
*
*        ==== NWMAX = the largest possible deflation window for
*        .    which there is sufficient workspace. ====
*
         nwmax = min( ( n-1 ) / 3, lwork / 2 )
         nw = nwmax
*
*        ==== NSMAX = the Largest number of simultaneous shifts
*        .    for which there is sufficient workspace. ====
*
         nsmax = min( ( n-3 ) / 6, 2*lwork / 3 )
         nsmax = nsmax - mod( nsmax, 2 )
*
*        ==== NDFL: an iteration count restarted at deflation. ====
*
         ndfl = 1
*
*        ==== ITMAX = iteration limit ====
*
         itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
*
*        ==== Last row and column in the active block ====
*
         kbot = ihi
*
*        ==== Main Loop ====
*
         DO 70 it = 1, itmax
*
*           ==== Done when KBOT falls below ILO ====
*
            IF( kbot.LT.ilo )
     $         GO TO 80
*
*           ==== Locate active block ====
*
            DO 10 k = kbot, ilo + 1, -1
               IF( h( k, k-1 ).EQ.zero )
     $            GO TO 20
   10       CONTINUE
            k = ilo
   20       CONTINUE
            ktop = k
*
*           ==== Select deflation window size:
*           .    Typical Case:
*           .      If possible and advisable, nibble the entire
*           .      active block.  If not, use size MIN(NWR,NWMAX)
*           .      or MIN(NWR+1,NWMAX) depending upon which has
*           .      the smaller corresponding subdiagonal entry
*           .      (a heuristic).
*           .
*           .    Exceptional Case:
*           .      If there have been no deflations in KEXNW or
*           .      more iterations, then vary the deflation window
*           .      size.   At first, because, larger windows are,
*           .      in general, more powerful than smaller ones,
*           .      rapidly increase the window to the maximum possible.
*           .      Then, gradually reduce the window size. ====
*
            nh = kbot - ktop + 1
            nwupbd = min( nh, nwmax )
            IF( ndfl.LT.kexnw ) THEN
               nw = min( nwupbd, nwr )
            ELSE
               nw = min( nwupbd, 2*nw )
            END IF
            IF( nw.LT.nwmax ) THEN
               IF( nw.GE.nh-1 ) THEN
                  nw = nh
               ELSE
                  kwtop = kbot - nw + 1
                  IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
     $                cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
               END IF
            END IF
            IF( ndfl.LT.kexnw ) THEN
               ndec = -1
            ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
               ndec = ndec + 1
               IF( nw-ndec.LT.2 )
     $            ndec = 0
               nw = nw - ndec
            END IF
*
*           ==== Aggressive early deflation:
*           .    split workspace under the subdiagonal into
*           .      - an nw-by-nw work array V in the lower
*           .        left-hand-corner,
*           .      - an NW-by-at-least-NW-but-more-is-better
*           .        (NW-by-NHO) horizontal work array along
*           .        the bottom edge,
*           .      - an at-least-NW-but-more-is-better (NHV-by-NW)
*           .        vertical work array along the left-hand-edge.
*           .        ====
*
            kv = n - nw + 1
            kt = nw + 1
            nho = ( n-nw-1 ) - kt + 1
            kwv = nw + 2
            nve = ( n-nw ) - kwv + 1
*
*           ==== Aggressive early deflation ====
*
            CALL zlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh,
     $                   iloz,
     $                   ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
     $                   h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
     $                   lwork )
*
*           ==== Adjust KBOT accounting for new deflations. ====
*
            kbot = kbot - ld
*
*           ==== KS points to the shifts. ====
*
            ks = kbot - ls + 1
*
*           ==== Skip an expensive QR sweep if there is a (partly
*           .    heuristic) reason to expect that many eigenvalues
*           .    will deflate without it.  Here, the QR sweep is
*           .    skipped if many eigenvalues have just been deflated
*           .    or if the remaining active block is small.
*
            IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
     $          ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
*
*              ==== NS = nominal number of simultaneous shifts.
*              .    This may be lowered (slightly) if ZLAQR2
*              .    did not provide that many shifts. ====
*
               ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
               ns = ns - mod( ns, 2 )
*
*              ==== If there have been no deflations
*              .    in a multiple of KEXSH iterations,
*              .    then try exceptional shifts.
*              .    Otherwise use shifts provided by
*              .    ZLAQR2 above or from the eigenvalues
*              .    of a trailing principal submatrix. ====
*
               IF( mod( ndfl, kexsh ).EQ.0 ) THEN
                  ks = kbot - ns + 1
                  DO 30 i = kbot, ks + 1, -2
                     w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
                     w( i-1 ) = w( i )
   30             CONTINUE
               ELSE
*
*                 ==== Got NS/2 or fewer shifts? Use ZLAHQR
*                 .    on a trailing principal submatrix to
*                 .    get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
*                 .    there is enough space below the subdiagonal
*                 .    to fit an NS-by-NS scratch array.) ====
*
                  IF( kbot-ks+1.LE.ns / 2 ) THEN
                     ks = kbot - ns + 1
                     kt = n - ns + 1
                     CALL zlacpy( 'A', ns, ns, h( ks, ks ), ldh,
     $                            h( kt, 1 ), ldh )
                     CALL zlahqr( .false., .false., ns, 1, ns,
     $                            h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
     $                            1, inf )
                     ks = ks + inf
*
*                    ==== In case of a rare QR failure use
*                    .    eigenvalues of the trailing 2-by-2
*                    .    principal submatrix.  Scale to avoid
*                    .    overflows, underflows and subnormals.
*                    .    (The scale factor S can not be zero,
*                    .    because H(KBOT,KBOT-1) is nonzero.) ====
*
                     IF( ks.GE.kbot ) THEN
                        s = cabs1( h( kbot-1, kbot-1 ) ) +
     $                      cabs1( h( kbot, kbot-1 ) ) +
     $                      cabs1( h( kbot-1, kbot ) ) +
     $                      cabs1( h( kbot, kbot ) )
                        aa = h( kbot-1, kbot-1 ) / s
                        cc = h( kbot, kbot-1 ) / s
                        bb = h( kbot-1, kbot ) / s
                        dd = h( kbot, kbot ) / s
                        tr2 = ( aa+dd ) / two
                        det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
                        rtdisc = sqrt( -det )
                        w( kbot-1 ) = ( tr2+rtdisc )*s
                        w( kbot ) = ( tr2-rtdisc )*s
*
                        ks = kbot - 1
                     END IF
                  END IF
*
                  IF( kbot-ks+1.GT.ns ) THEN
*
*                    ==== Sort the shifts (Helps a little) ====
*
                     sorted = .false.
                     DO 50 k = kbot, ks + 1, -1
                        IF( sorted )
     $                     GO TO 60
                        sorted = .true.
                        DO 40 i = ks, k - 1
                           IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
     $                          THEN
                              sorted = .false.
                              swap = w( i )
                              w( i ) = w( i+1 )
                              w( i+1 ) = swap
                           END IF
   40                   CONTINUE
   50                CONTINUE
   60                CONTINUE
                  END IF
               END IF
*
*              ==== If there are only two shifts, then use
*              .    only one.  ====
*
               IF( kbot-ks+1.EQ.2 ) THEN
                  IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
     $                cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
                     w( kbot-1 ) = w( kbot )
                  ELSE
                     w( kbot ) = w( kbot-1 )
                  END IF
               END IF
*
*              ==== Use up to NS of the the smallest magnitude
*              .    shifts.  If there aren't NS shifts available,
*              .    then use them all, possibly dropping one to
*              .    make the number of shifts even. ====
*
               ns = min( ns, kbot-ks+1 )
               ns = ns - mod( ns, 2 )
               ks = kbot - ns + 1
*
*              ==== Small-bulge multi-shift QR sweep:
*              .    split workspace under the subdiagonal into
*              .    - a KDU-by-KDU work array U in the lower
*              .      left-hand-corner,
*              .    - a KDU-by-at-least-KDU-but-more-is-better
*              .      (KDU-by-NHo) horizontal work array WH along
*              .      the bottom edge,
*              .    - and an at-least-KDU-but-more-is-better-by-KDU
*              .      (NVE-by-KDU) vertical work WV arrow along
*              .      the left-hand-edge. ====
*
               kdu = 2*ns
               ku = n - kdu + 1
               kwh = kdu + 1
               nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
               kwv = kdu + 4
               nve = n - kdu - kwv + 1
*
*              ==== Small-bulge multi-shift QR sweep ====
*
               CALL zlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
     $                      w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
     $                      3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
     $                      nho, h( ku, kwh ), ldh )
            END IF
*
*           ==== Note progress (or the lack of it). ====
*
            IF( ld.GT.0 ) THEN
               ndfl = 1
            ELSE
               ndfl = ndfl + 1
            END IF
*
*           ==== End of main loop ====
   70    CONTINUE
*
*        ==== Iteration limit exceeded.  Set INFO to show where
*        .    the problem occurred and exit. ====
*
         info = kbot
   80    CONTINUE
      END IF
*
*     ==== Return the optimal value of LWORK. ====
*
      work( 1 ) = dcmplx( lwkopt, 0 )
*
*     ==== End of ZLAQR4 ====
*

Here is the call graph for this function:

Here is the caller graph for this function: