241 SUBROUTINE zlaqr0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
242 $ ihiz, z, ldz, work, lwork, info )
250 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
254 COMPLEX*16 H( ldh, * ), W( * ), WORK( * ), Z( ldz, * )
265 parameter ( ntiny = 11 )
271 parameter ( kexnw = 5 )
277 parameter ( kexsh = 6 )
281 DOUBLE PRECISION WILK1
282 parameter ( wilk1 = 0.75d0 )
284 parameter ( zero = ( 0.0d0, 0.0d0 ),
285 $ one = ( 1.0d0, 0.0d0 ) )
287 parameter ( two = 2.0d0 )
290 COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
292 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
293 $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
294 $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
295 $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
304 COMPLEX*16 ZDUM( 1, 1 )
310 INTRINSIC abs, dble, dcmplx, dimag, int, max, min, mod,
314 DOUBLE PRECISION CABS1
317 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
329 IF( n.LE.ntiny )
THEN
335 $
CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
336 $ ihiz, z, ldz, info )
365 nwr = ilaenv( 13,
'ZLAQR0', jbcmpz, n, ilo, ihi, lwork )
367 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
374 nsr = ilaenv( 15,
'ZLAQR0', jbcmpz, n, ilo, ihi, lwork )
375 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
376 nsr = max( 2, nsr-mod( nsr, 2 ) )
382 CALL zlaqr3( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
383 $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
388 lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
392 IF( lwork.EQ.-1 )
THEN
393 work( 1 ) = dcmplx( lwkopt, 0 )
399 nmin = ilaenv( 12,
'ZLAQR0', jbcmpz, n, ilo, ihi, lwork )
400 nmin = max( ntiny, nmin )
404 nibble = ilaenv( 14,
'ZLAQR0', jbcmpz, n, ilo, ihi, lwork )
405 nibble = max( 0, nibble )
410 kacc22 = ilaenv( 16,
'ZLAQR0', jbcmpz, n, ilo, ihi, lwork )
411 kacc22 = max( 0, kacc22 )
412 kacc22 = min( 2, kacc22 )
417 nwmax = min( ( n-1 ) / 3, lwork / 2 )
423 nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
424 nsmax = nsmax - mod( nsmax, 2 )
432 itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
449 DO 10 k = kbot, ilo + 1, -1
450 IF( h( k, k-1 ).EQ.zero )
474 nwupbd = min( nh, nwmax )
475 IF( ndfl.LT.kexnw )
THEN
476 nw = min( nwupbd, nwr )
478 nw = min( nwupbd, 2*nw )
480 IF( nw.LT.nwmax )
THEN
481 IF( nw.GE.nh-1 )
THEN
484 kwtop = kbot - nw + 1
485 IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
486 $ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
489 IF( ndfl.LT.kexnw )
THEN
491 ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd )
THEN
511 nho = ( n-nw-1 ) - kt + 1
513 nve = ( n-nw ) - kwv + 1
517 CALL zlaqr3( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
518 $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
519 $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
536 IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
537 $ ktop+1.GT.min( nmin, nwmax ) ) ) )
THEN
543 ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
544 ns = ns - mod( ns, 2 )
553 IF( mod( ndfl, kexsh ).EQ.0 )
THEN
555 DO 30 i = kbot, ks + 1, -2
556 w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
567 IF( kbot-ks+1.LE.ns / 2 )
THEN
570 CALL zlacpy(
'A', ns, ns, h( ks, ks ), ldh,
572 IF( ns.GT.nmin )
THEN
573 CALL zlaqr4( .false., .false., ns, 1, ns,
574 $ h( kt, 1 ), ldh, w( ks ), 1, 1,
575 $ zdum, 1, work, lwork, inf )
577 CALL zlahqr( .false., .false., ns, 1, ns,
578 $ h( kt, 1 ), ldh, w( ks ), 1, 1,
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 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.
subroutine zlaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine zlaqr3(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)
ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...