256 SUBROUTINE dlaqr0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
257 $ iloz, ihiz, z, ldz, work, lwork, info )
265 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
269 DOUBLE PRECISION H( ldh, * ), WI( * ), WORK( * ), WR( * ),
281 parameter ( ntiny = 11 )
287 parameter ( kexnw = 5 )
293 parameter ( kexsh = 6 )
297 DOUBLE PRECISION WILK1, WILK2
298 parameter ( wilk1 = 0.75d0, wilk2 = -0.4375d0 )
299 DOUBLE PRECISION ZERO, ONE
300 parameter ( zero = 0.0d0, one = 1.0d0 )
303 DOUBLE PRECISION AA, BB, CC, CS, DD, SN, SS, SWAP
304 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
305 $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
306 $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
307 $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
316 DOUBLE PRECISION ZDUM( 1, 1 )
322 INTRINSIC abs, dble, int, max, min, mod
334 IF( n.LE.ntiny )
THEN
340 $
CALL dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
341 $ iloz, ihiz, z, ldz, info )
370 nwr = ilaenv( 13,
'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
372 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
379 nsr = ilaenv( 15,
'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
380 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
381 nsr = max( 2, nsr-mod( nsr, 2 ) )
387 CALL dlaqr3( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
388 $ ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
389 $ n, h, ldh, work, -1 )
393 lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
397 IF( lwork.EQ.-1 )
THEN
398 work( 1 ) = dble( lwkopt )
404 nmin = ilaenv( 12,
'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
405 nmin = max( ntiny, nmin )
409 nibble = ilaenv( 14,
'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
410 nibble = max( 0, nibble )
415 kacc22 = ilaenv( 16,
'DLAQR0', jbcmpz, n, ilo, ihi, lwork )
416 kacc22 = max( 0, kacc22 )
417 kacc22 = min( 2, kacc22 )
422 nwmax = min( ( n-1 ) / 3, lwork / 2 )
428 nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
429 nsmax = nsmax - mod( nsmax, 2 )
437 itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
454 DO 10 k = kbot, ilo + 1, -1
455 IF( h( k, k-1 ).EQ.zero )
479 nwupbd = min( nh, nwmax )
480 IF( ndfl.LT.kexnw )
THEN
481 nw = min( nwupbd, nwr )
483 nw = min( nwupbd, 2*nw )
485 IF( nw.LT.nwmax )
THEN
486 IF( nw.GE.nh-1 )
THEN
489 kwtop = kbot - nw + 1
490 IF( abs( h( kwtop, kwtop-1 ) ).GT.
491 $ abs( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
494 IF( ndfl.LT.kexnw )
THEN
496 ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd )
THEN
516 nho = ( n-nw-1 ) - kt + 1
518 nve = ( n-nw ) - kwv + 1
522 CALL dlaqr3( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
523 $ ihiz, z, ldz, ls, ld, wr, wi, h( kv, 1 ), ldh,
524 $ nho, h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh,
541 IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
542 $ ktop+1.GT.min( nmin, nwmax ) ) ) )
THEN
548 ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
549 ns = ns - mod( ns, 2 )
558 IF( mod( ndfl, kexsh ).EQ.0 )
THEN
560 DO 30 i = kbot, max( ks+1, ktop+2 ), -2
561 ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
562 aa = wilk1*ss + h( i, i )
566 CALL dlanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
567 $ wr( i ), wi( i ), cs, sn )
569 IF( ks.EQ.ktop )
THEN
570 wr( ks+1 ) = h( ks+1, ks+1 )
572 wr( ks ) = wr( ks+1 )
573 wi( ks ) = wi( ks+1 )
583 IF( kbot-ks+1.LE.ns / 2 )
THEN
586 CALL dlacpy(
'A', ns, ns, h( ks, ks ), ldh,
588 IF( ns.GT.nmin )
THEN
589 CALL dlaqr4( .false., .false., ns, 1, ns,
590 $ h( kt, 1 ), ldh, wr( ks ),
591 $ wi( ks ), 1, 1, zdum, 1, work,
594 CALL dlahqr( .false., .false., ns, 1, ns,
595 $ h( kt, 1 ), ldh, wr( ks ),
596 $ wi( ks ), 1, 1, zdum, 1, inf )
604 IF( ks.GE.kbot )
THEN
605 aa = h( kbot-1, kbot-1 )
606 cc = h( kbot, kbot-1 )
607 bb = h( kbot-1, kbot )
609 CALL dlanv2( aa, bb, cc, dd, wr( kbot-1 ),
610 $ wi( kbot-1 ), wr( kbot ),
611 $ wi( kbot ), cs, sn )
616 IF( kbot-ks+1.GT.ns )
THEN
623 DO 50 k = kbot, ks + 1, -1
628 IF( abs( wr( i ) )+abs( wi( i ) ).LT.
629 $ abs( wr( i+1 ) )+abs( wi( i+1 ) ) )
THEN
651 DO 70 i = kbot, ks + 2, -2
652 IF( wi( i ).NE.-wi( i-1 ) )
THEN
656 wr( i-1 ) = wr( i-2 )
661 wi( i-1 ) = wi( i-2 )
670 IF( kbot-ks+1.EQ.2 )
THEN
671 IF( wi( kbot ).EQ.zero )
THEN
672 IF( abs( wr( kbot )-h( kbot, kbot ) ).LT.
673 $ abs( wr( kbot-1 )-h( kbot, kbot ) ) )
THEN
674 wr( kbot-1 ) = wr( kbot )
676 wr( kbot ) = wr( kbot-1 )
686 ns = min( ns, kbot-ks+1 )
687 ns = ns - mod( ns, 2 )
704 nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
706 nve = n - kdu - kwv + 1
710 CALL dlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
711 $ wr( ks ), wi( ks ), h, ldh, iloz, ihiz, z,
712 $ ldz, work, 3, h( ku, 1 ), ldh, nve,
713 $ h( kwv, 1 ), ldh, nho, h( ku, kwh ), ldh )
736 work( 1 ) = dble( lwkopt )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine dlaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine dlaqr3(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
subroutine dlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine dlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
DLAQR5 performs a single small-bulge multi-shift QR sweep.
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form...