207 SUBROUTINE dlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
208 $ iloz, ihiz, z, ldz, info )
216 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
220 DOUBLE PRECISION H( ldh, * ), WI( * ), WR( * ), Z( ldz, * )
226 DOUBLE PRECISION ZERO, ONE, TWO
227 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
228 DOUBLE PRECISION DAT1, DAT2
229 parameter ( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
232 DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233 $ h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
234 $ safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
236 INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ
239 DOUBLE PRECISION V( 3 )
242 DOUBLE PRECISION DLAMCH
249 INTRINSIC abs, dble, max, min, sqrt
259 IF( ilo.EQ.ihi )
THEN
260 wr( ilo ) = h( ilo, ilo )
266 DO 10 j = ilo, ihi - 3
271 $ h( ihi, ihi-2 ) = zero
278 safmin = dlamch(
'SAFE MINIMUM' )
279 safmax = one / safmin
280 CALL dlabad( safmin, safmax )
281 ulp = dlamch(
'PRECISION' )
282 smlnum = safmin*( dble( nh ) / ulp )
295 itmax = 30 * max( 10, nh )
313 DO 140 its = 0, itmax
317 DO 30 k = i, l + 1, -1
318 IF( abs( h( k, k-1 ) ).LE.smlnum )
320 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
321 IF( tst.EQ.zero )
THEN
323 $ tst = tst + abs( h( k-1, k-2 ) )
325 $ tst = tst + abs( h( k+1, k ) )
331 IF( abs( h( k, k-1 ) ).LE.ulp*tst )
THEN
332 ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
333 ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
334 aa = max( abs( h( k, k ) ),
335 $ abs( h( k-1, k-1 )-h( k, k ) ) )
336 bb = min( abs( h( k, k ) ),
337 $ abs( h( k-1, k-1 )-h( k, k ) ) )
339 IF( ba*( ab / s ).LE.max( smlnum,
340 $ ulp*( bb*( aa / s ) ) ) )
GO TO 40
361 IF( .NOT.wantt )
THEN
370 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
371 h11 = dat1*s + h( l, l )
375 ELSE IF( its.EQ.20 )
THEN
379 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
380 h11 = dat1*s + h( i, i )
394 s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
405 tr = ( h11+h22 ) / two
406 det = ( h11-tr )*( h22-tr ) - h12*h21
407 rtdisc = sqrt( abs( det ) )
408 IF( det.GE.zero )
THEN
422 IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) )
THEN
436 DO 50 m = i - 2, l, -1
443 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
444 h21s = h( m+1, m ) / s
445 v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
446 $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
447 v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
448 v( 3 ) = h21s*h( m+2, m+1 )
449 s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
455 IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
456 $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
457 $ m ) )+abs( h( m+1, m+1 ) ) ) )
GO TO 60
476 $
CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
477 CALL dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
482 $ h( k+2, k-1 ) = zero
483 ELSE IF( m.GT.l )
THEN
488 h( k, k-1 ) = h( k, k-1 )*( one-t1 )
500 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
501 h( k, j ) = h( k, j ) - sum*t1
502 h( k+1, j ) = h( k+1, j ) - sum*t2
503 h( k+2, j ) = h( k+2, j ) - sum*t3
509 DO 80 j = i1, min( k+3, i )
510 sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
511 h( j, k ) = h( j, k ) - sum*t1
512 h( j, k+1 ) = h( j, k+1 ) - sum*t2
513 h( j, k+2 ) = h( j, k+2 ) - sum*t3
521 sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
522 z( j, k ) = z( j, k ) - sum*t1
523 z( j, k+1 ) = z( j, k+1 ) - sum*t2
524 z( j, k+2 ) = z( j, k+2 ) - sum*t3
527 ELSE IF( nr.EQ.2 )
THEN
533 sum = h( k, j ) + v2*h( k+1, j )
534 h( k, j ) = h( k, j ) - sum*t1
535 h( k+1, j ) = h( k+1, j ) - sum*t2
542 sum = h( j, k ) + v2*h( j, k+1 )
543 h( j, k ) = h( j, k ) - sum*t1
544 h( j, k+1 ) = h( j, k+1 ) - sum*t2
551 DO 120 j = iloz, ihiz
552 sum = z( j, k ) + v2*z( j, k+1 )
553 z( j, k ) = z( j, k ) - sum*t1
554 z( j, k+1 ) = z( j, k+1 ) - sum*t2
575 ELSE IF( l.EQ.i-1 )
THEN
582 CALL dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
583 $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
591 $
CALL drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
593 CALL drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
599 CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
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 dlabad(SMALL, LARGE)
DLABAD
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
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...