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, * )
227 parameter( itmax = 30 )
228 DOUBLE PRECISION ZERO, ONE, TWO
229 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
230 DOUBLE PRECISION DAT1, DAT2
231 parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
234 DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
235 $ h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
236 $ safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
238 INTEGER I, I1, I2, ITS, J, K, L, M, NH, NR, NZ
241 DOUBLE PRECISION V( 3 )
244 DOUBLE PRECISION DLAMCH
251 INTRINSIC abs, dble, max, min, sqrt
261 IF( ilo.EQ.ihi )
THEN
262 wr( ilo ) = h( ilo, ilo )
268 DO 10 j = ilo, ihi - 3
273 $ h( ihi, ihi-2 ) = zero
280 safmin = dlamch(
'SAFE MINIMUM' )
281 safmax = one / safmin
282 CALL
dlabad( safmin, safmax )
283 ulp = dlamch(
'PRECISION' )
284 smlnum = safmin*( dble( nh ) / ulp )
311 DO 140 its = 0, itmax
315 DO 30 k = i, l + 1, -1
316 IF( abs( h( k, k-1 ) ).LE.smlnum )
318 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
319 IF( tst.EQ.zero )
THEN
321 $ tst = tst + abs( h( k-1, k-2 ) )
323 $ tst = tst + abs( h( k+1, k ) )
329 IF( abs( h( k, k-1 ) ).LE.ulp*tst )
THEN
330 ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
331 ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
332 aa = max( abs( h( k, k ) ),
333 $ abs( h( k-1, k-1 )-h( k, k ) ) )
334 bb = min( abs( h( k, k ) ),
335 $ abs( h( k-1, k-1 )-h( k, k ) ) )
337 IF( ba*( ab / s ).LE.max( smlnum,
338 $ ulp*( bb*( aa / s ) ) ) )go to 40
359 IF( .NOT.wantt )
THEN
368 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
369 h11 = dat1*s + h( l, l )
373 ELSE IF( its.EQ.20 )
THEN
377 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
378 h11 = dat1*s + h( i, i )
392 s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
403 tr = ( h11+h22 ) / two
404 det = ( h11-tr )*( h22-tr ) - h12*h21
405 rtdisc = sqrt( abs( det ) )
406 IF( det.GE.zero )
THEN
420 IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) )
THEN
434 DO 50 m = i - 2, l, -1
441 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
442 h21s = h( m+1, m ) / s
443 v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
444 $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
445 v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
446 v( 3 ) = h21s*h( m+2, m+1 )
447 s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
453 IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
454 $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
455 $ m ) )+abs( h( m+1, m+1 ) ) ) )go to 60
474 $ CALL
dcopy( nr, h( k, k-1 ), 1, v, 1 )
475 CALL
dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
480 $ h( k+2, k-1 ) = zero
481 ELSE IF( m.GT.l )
THEN
486 h( k, k-1 ) = h( k, k-1 )*( one-t1 )
498 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
499 h( k, j ) = h( k, j ) - sum*t1
500 h( k+1, j ) = h( k+1, j ) - sum*t2
501 h( k+2, j ) = h( k+2, j ) - sum*t3
507 DO 80 j = i1, min( k+3, i )
508 sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
509 h( j, k ) = h( j, k ) - sum*t1
510 h( j, k+1 ) = h( j, k+1 ) - sum*t2
511 h( j, k+2 ) = h( j, k+2 ) - sum*t3
519 sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
520 z( j, k ) = z( j, k ) - sum*t1
521 z( j, k+1 ) = z( j, k+1 ) - sum*t2
522 z( j, k+2 ) = z( j, k+2 ) - sum*t3
525 ELSE IF( nr.EQ.2 )
THEN
531 sum = h( k, j ) + v2*h( k+1, j )
532 h( k, j ) = h( k, j ) - sum*t1
533 h( k+1, j ) = h( k+1, j ) - sum*t2
540 sum = h( j, k ) + v2*h( j, k+1 )
541 h( j, k ) = h( j, k ) - sum*t1
542 h( j, k+1 ) = h( j, k+1 ) - sum*t2
549 DO 120 j = iloz, ihiz
550 sum = z( j, k ) + v2*z( j, k+1 )
551 z( j, k ) = z( j, k ) - sum*t1
552 z( j, k+1 ) = z( j, k+1 ) - sum*t2
573 ELSE IF( l.EQ.i-1 )
THEN
580 CALL
dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
581 $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
589 $ CALL
drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
591 CALL
drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
597 CALL
drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 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...
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT