195 SUBROUTINE zlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196 $ ihiz, z, ldz, info )
204 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
208 COMPLEX*16 H( ldh, * ), W( * ), Z( ldz, * )
215 parameter ( zero = ( 0.0d0, 0.0d0 ),
216 $ one = ( 1.0d0, 0.0d0 ) )
217 DOUBLE PRECISION RZERO, RONE, HALF
218 parameter ( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
219 DOUBLE PRECISION DAT1
220 parameter ( dat1 = 3.0d0 / 4.0d0 )
223 COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
225 DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226 $ safmin, smlnum, sx, t2, tst, ulp
227 INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
235 DOUBLE PRECISION DLAMCH
236 EXTERNAL zladiv, dlamch
242 DOUBLE PRECISION CABS1
245 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
248 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
258 IF( ilo.EQ.ihi )
THEN
259 w( ilo ) = h( ilo, ilo )
264 DO 10 j = ilo, ihi - 3
269 $ h( ihi, ihi-2 ) = zero
278 DO 20 i = ilo + 1, ihi
279 IF( dimag( h( i, i-1 ) ).NE.rzero )
THEN
283 sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
284 sc = dconjg( sc ) / abs( sc )
285 h( i, i-1 ) = abs( h( i, i-1 ) )
286 CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
287 CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
290 $
CALL zscal( ihiz-iloz+1, dconjg( sc ), z( iloz, i ), 1 )
299 safmin = dlamch(
'SAFE MINIMUM' )
300 safmax = rone / safmin
301 CALL dlabad( safmin, safmax )
302 ulp = dlamch(
'PRECISION' )
303 smlnum = safmin*( dble( nh ) / ulp )
316 itmax = 30 * max( 10, nh )
334 DO 130 its = 0, itmax
338 DO 40 k = i, l + 1, -1
339 IF( cabs1( h( k, k-1 ) ).LE.smlnum )
341 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
342 IF( tst.EQ.zero )
THEN
344 $ tst = tst + abs( dble( h( k-1, k-2 ) ) )
346 $ tst = tst + abs( dble( h( k+1, k ) ) )
352 IF( abs( dble( h( k, k-1 ) ) ).LE.ulp*tst )
THEN
353 ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
354 ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
355 aa = max( cabs1( h( k, k ) ),
356 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
357 bb = min( cabs1( h( k, k ) ),
358 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
360 IF( ba*( ab / s ).LE.max( smlnum,
361 $ ulp*( bb*( aa / s ) ) ) )
GO TO 50
382 IF( .NOT.wantt )
THEN
391 s = dat1*abs( dble( h( l+1, l ) ) )
393 ELSE IF( its.EQ.20 )
THEN
397 s = dat1*abs( dble( h( i, i-1 ) ) )
404 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
406 IF( s.NE.rzero )
THEN
407 x = half*( h( i-1, i-1 )-t )
409 s = max( s, cabs1( x ) )
410 y = s*sqrt( ( x / s )**2+( u / s )**2 )
411 IF( sx.GT.rzero )
THEN
412 IF( dble( x / sx )*dble( y )+dimag( x / sx )*
413 $ dimag( y ).LT.rzero )y = -y
415 t = t - u*zladiv( u, ( x+y ) )
421 DO 60 m = i - 1, l + 1, -1
430 h21 = dble( h( m+1, m ) )
431 s = cabs1( h11s ) + abs( h21 )
436 h10 = dble( h( m, m-1 ) )
437 IF( abs( h10 )*abs( h21 ).LE.ulp*
438 $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
444 h21 = dble( h( l+1, l ) )
445 s = cabs1( h11s ) + abs( h21 )
469 $
CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
470 CALL zlarfg( 2, v( 1 ), v( 2 ), 1, t1 )
482 sum = dconjg( t1 )*h( k, j ) + t2*h( k+1, j )
483 h( k, j ) = h( k, j ) - sum
484 h( k+1, j ) = h( k+1, j ) - sum*v2
490 DO 90 j = i1, min( k+2, i )
491 sum = t1*h( j, k ) + t2*h( j, k+1 )
492 h( j, k ) = h( j, k ) - sum
493 h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
500 DO 100 j = iloz, ihiz
501 sum = t1*z( j, k ) + t2*z( j, k+1 )
502 z( j, k ) = z( j, k ) - sum
503 z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
507 IF( k.EQ.m .AND. m.GT.l )
THEN
515 temp = temp / abs( temp )
516 h( m+1, m ) = h( m+1, m )*dconjg( temp )
518 $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
522 $
CALL zscal( i2-j, temp, h( j, j+1 ), ldh )
523 CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
525 CALL zscal( nz, dconjg( temp ), z( iloz, j ),
536 IF( dimag( temp ).NE.rzero )
THEN
541 $
CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
542 CALL zscal( i-i1, temp, h( i1, i ), 1 )
544 CALL zscal( nz, temp, z( iloz, i ), 1 )
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
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 dlabad(SMALL, LARGE)
DLABAD
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL