249 SUBROUTINE claqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
250 $ ihiz, z, ldz, work, lwork, info )
258 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
262 COMPLEX H( ldh, * ), W( * ), WORK( * ), Z( ldz, * )
274 parameter ( ntiny = 11 )
280 parameter ( kexnw = 5 )
286 parameter ( kexsh = 6 )
291 parameter ( wilk1 = 0.75e0 )
293 parameter ( zero = ( 0.0e0, 0.0e0 ),
294 $ one = ( 1.0e0, 0.0e0 ) )
296 parameter ( two = 2.0e0 )
299 COMPLEX AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
301 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
302 $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
303 $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
304 $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
319 INTRINSIC abs, aimag, cmplx, int, max, min, mod,
REAL,
326 cabs1( cdum ) = abs(
REAL( CDUM ) ) + abs( AIMAG( cdum ) )
338 IF( n.LE.ntiny )
THEN
344 $
CALL clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
345 $ ihiz, z, ldz, info )
374 nwr = ilaenv( 13,
'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
376 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
383 nsr = ilaenv( 15,
'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
384 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
385 nsr = max( 2, nsr-mod( nsr, 2 ) )
391 CALL claqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
392 $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
397 lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
401 IF( lwork.EQ.-1 )
THEN
402 work( 1 ) = cmplx( lwkopt, 0 )
408 nmin = ilaenv( 12,
'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
409 nmin = max( ntiny, nmin )
413 nibble = ilaenv( 14,
'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
414 nibble = max( 0, nibble )
419 kacc22 = ilaenv( 16,
'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
420 kacc22 = max( 0, kacc22 )
421 kacc22 = min( 2, kacc22 )
426 nwmax = min( ( n-1 ) / 3, lwork / 2 )
432 nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
433 nsmax = nsmax - mod( nsmax, 2 )
441 itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
458 DO 10 k = kbot, ilo + 1, -1
459 IF( h( k, k-1 ).EQ.zero )
483 nwupbd = min( nh, nwmax )
484 IF( ndfl.LT.kexnw )
THEN
485 nw = min( nwupbd, nwr )
487 nw = min( nwupbd, 2*nw )
489 IF( nw.LT.nwmax )
THEN
490 IF( nw.GE.nh-1 )
THEN
493 kwtop = kbot - nw + 1
494 IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
495 $ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
498 IF( ndfl.LT.kexnw )
THEN
500 ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd )
THEN
520 nho = ( n-nw-1 ) - kt + 1
522 nve = ( n-nw ) - kwv + 1
526 CALL claqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
527 $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
528 $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
545 IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
546 $ ktop+1.GT.min( nmin, nwmax ) ) ) )
THEN
552 ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
553 ns = ns - mod( ns, 2 )
562 IF( mod( ndfl, kexsh ).EQ.0 )
THEN
564 DO 30 i = kbot, ks + 1, -2
565 w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
576 IF( kbot-ks+1.LE.ns / 2 )
THEN
579 CALL clacpy(
'A', ns, ns, h( ks, ks ), ldh,
581 CALL clahqr( .false., .false., ns, 1, ns,
582 $ h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
593 IF( ks.GE.kbot )
THEN
594 s = cabs1( h( kbot-1, kbot-1 ) ) +
595 $ cabs1( h( kbot, kbot-1 ) ) +
596 $ cabs1( h( kbot-1, kbot ) ) +
597 $ cabs1( h( kbot, kbot ) )
598 aa = h( kbot-1, kbot-1 ) / s
599 cc = h( kbot, kbot-1 ) / s
600 bb = h( kbot-1, kbot ) / s
601 dd = h( kbot, kbot ) / s
602 tr2 = ( aa+dd ) / two
603 det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
604 rtdisc = sqrt( -det )
605 w( kbot-1 ) = ( tr2+rtdisc )*s
606 w( kbot ) = ( tr2-rtdisc )*s
612 IF( kbot-ks+1.GT.ns )
THEN
617 DO 50 k = kbot, ks + 1, -1
622 IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
638 IF( kbot-ks+1.EQ.2 )
THEN
639 IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
640 $ cabs1( w( kbot-1 )-h( kbot, kbot ) ) )
THEN
641 w( kbot-1 ) = w( kbot )
643 w( kbot ) = w( kbot-1 )
652 ns = min( ns, kbot-ks+1 )
653 ns = ns - mod( ns, 2 )
670 nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
672 nve = n - kdu - kwv + 1
676 CALL claqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
677 $ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
678 $ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
679 $ nho, h( ku, kwh ), ldh )
702 work( 1 ) = cmplx( lwkopt, 0 )
subroutine clahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine claqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
CLAQR5 performs a single small-bulge multi-shift QR sweep.
subroutine claqr2(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)
CLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
subroutine claqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.