247 SUBROUTINE zlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
248 $ ihiz, z, ldz, work, lwork, info )
256 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
260 COMPLEX*16 H( ldh, * ), W( * ), WORK( * ), Z( ldz, * )
271 parameter ( ntiny = 11 )
277 parameter ( kexnw = 5 )
283 parameter ( kexsh = 6 )
287 DOUBLE PRECISION WILK1
288 parameter ( wilk1 = 0.75d0 )
290 parameter ( zero = ( 0.0d0, 0.0d0 ),
291 $ one = ( 1.0d0, 0.0d0 ) )
293 parameter ( two = 2.0d0 )
296 COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
298 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
299 $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
300 $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
301 $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
310 COMPLEX*16 ZDUM( 1, 1 )
316 INTRINSIC abs, dble, dcmplx, dimag, int, max, min, mod,
320 DOUBLE PRECISION CABS1
323 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
335 IF( n.LE.ntiny )
THEN
341 $
CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
342 $ ihiz, z, ldz, info )
371 nwr = ilaenv( 13,
'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
373 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
380 nsr = ilaenv( 15,
'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
381 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
382 nsr = max( 2, nsr-mod( nsr, 2 ) )
388 CALL zlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
389 $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
394 lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
398 IF( lwork.EQ.-1 )
THEN
399 work( 1 ) = dcmplx( lwkopt, 0 )
405 nmin = ilaenv( 12,
'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
406 nmin = max( ntiny, nmin )
410 nibble = ilaenv( 14,
'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
411 nibble = max( 0, nibble )
416 kacc22 = ilaenv( 16,
'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
417 kacc22 = max( 0, kacc22 )
418 kacc22 = min( 2, kacc22 )
423 nwmax = min( ( n-1 ) / 3, lwork / 2 )
429 nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
430 nsmax = nsmax - mod( nsmax, 2 )
438 itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
455 DO 10 k = kbot, ilo + 1, -1
456 IF( h( k, k-1 ).EQ.zero )
480 nwupbd = min( nh, nwmax )
481 IF( ndfl.LT.kexnw )
THEN
482 nw = min( nwupbd, nwr )
484 nw = min( nwupbd, 2*nw )
486 IF( nw.LT.nwmax )
THEN
487 IF( nw.GE.nh-1 )
THEN
490 kwtop = kbot - nw + 1
491 IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
492 $ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
495 IF( ndfl.LT.kexnw )
THEN
497 ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd )
THEN
517 nho = ( n-nw-1 ) - kt + 1
519 nve = ( n-nw ) - kwv + 1
523 CALL zlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
524 $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
525 $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
542 IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
543 $ ktop+1.GT.min( nmin, nwmax ) ) ) )
THEN
549 ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
550 ns = ns - mod( ns, 2 )
559 IF( mod( ndfl, kexsh ).EQ.0 )
THEN
561 DO 30 i = kbot, ks + 1, -2
562 w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
573 IF( kbot-ks+1.LE.ns / 2 )
THEN
576 CALL zlacpy(
'A', ns, ns, h( ks, ks ), ldh,
578 CALL zlahqr( .false., .false., ns, 1, ns,
579 $ h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
590 IF( ks.GE.kbot )
THEN
591 s = cabs1( h( kbot-1, kbot-1 ) ) +
592 $ cabs1( h( kbot, kbot-1 ) ) +
593 $ cabs1( h( kbot-1, kbot ) ) +
594 $ cabs1( h( kbot, kbot ) )
595 aa = h( kbot-1, kbot-1 ) / s
596 cc = h( kbot, kbot-1 ) / s
597 bb = h( kbot-1, kbot ) / s
598 dd = h( kbot, kbot ) / s
599 tr2 = ( aa+dd ) / two
600 det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
601 rtdisc = sqrt( -det )
602 w( kbot-1 ) = ( tr2+rtdisc )*s
603 w( kbot ) = ( tr2-rtdisc )*s
609 IF( kbot-ks+1.GT.ns )
THEN
614 DO 50 k = kbot, ks + 1, -1
619 IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
635 IF( kbot-ks+1.EQ.2 )
THEN
636 IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
637 $ cabs1( w( kbot-1 )-h( kbot, kbot ) ) )
THEN
638 w( kbot-1 ) = w( kbot )
640 w( kbot ) = w( kbot-1 )
649 ns = min( ns, kbot-ks+1 )
650 ns = ns - mod( ns, 2 )
667 nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
669 nve = n - kdu - kwv + 1
673 CALL zlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
674 $ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
675 $ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
676 $ nho, h( ku, kwh ), ldh )
699 work( 1 ) = dcmplx( lwkopt, 0 )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine zlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
ZLAQR5 performs a single small-bulge multi-shift QR sweep.
subroutine zlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.