1 SUBROUTINE clahqr2( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
11 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
14 COMPLEX H( LDH, * ), W( * ), Z( LDZ, * )
99 parameter( zero = ( 0.0e+0, 0.0e+0 ) )
101 parameter( rzero = 0.0e+0, rone = 1.0e+0 )
103 parameter( dat1 = 0.75e+0, dat2 = -0.4375e+0 )
106 INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
107 REAL CS, OVFL, S, SMLNUM, TST1, ULP, UNFL
108 COMPLEX CDUM, H00, H10, H11, H12, H21, H22, H33, H33S,
109 $ h43h34, h44, h44s, sn, sum, t1, t2, t3, v1, v2,
118 EXTERNAL slamch, clanhs
121 EXTERNAL slabad, ccopy,
clanv2, clarfg, crot
124 INTRINSIC abs, real, conjg, aimag,
max,
min
130 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
140 IF( ilo.EQ.ihi )
THEN
141 w( ilo ) = h( ilo, ilo )
151 unfl = slamch(
'Safe minimum' )
153 CALL slabad( unfl, ovfl )
154 ulp = slamch(
'Precision' )
155 smlnum = unfl*( nh / ulp )
190 DO 20 k = i, l + 1, -1
191 tst1 = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
193 $ tst1 = clanhs(
'1', i-l+1, h( l, l ), ldh, rwork )
194 IF( cabs1( h( k, k-1 ) ).LE.
max( ulp*tst1, smlnum ) )
215 IF( .NOT.wantt )
THEN
220 IF( its.EQ.10 .OR. its.EQ.20 )
THEN
225 s = cabs1( h( i, i-1 ) ) + cabs1( h( i-1, i-2 ) )
235 h43h34 = h( i, i-1 )*h( i-1, i )
240 DO 40 m = i - 2, l, -1
252 v1 = ( h33s*h44s-h43h34 ) / h21 + h12
253 v2 = h22 - h11 - h33s - h44s
255 s = cabs1( v1 ) + cabs1( v2 ) + abs( v3 )
266 tst1 = cabs1( v1 )*( cabs1( h00 )+cabs1( h11 )+
268 IF( cabs1( h10 )*( cabs1( v2 )+cabs1( v3 ) ).LE.ulp*tst1 )
288 $
CALL ccopy( nr, h( k, k-1 ), 1, v, 1 )
289 CALL clarfg( nr, v( 1 ), v( 2 ), 1, t1 )
294 $ h( k+2, k-1 ) = zero
295 ELSE IF( m.GT.l )
THEN
298 h( k, k-1 ) = h( k, k-1 ) - conjg( t1 )*h( k, k-1 )
310 sum = conjg( t1 )*h( k, j ) +
311 $ conjg( t2 )*h( k+1, j ) +
312 $ conjg( t3 )*h( k+2, j )
313 h( k, j ) = h( k, j ) - sum
314 h( k+1, j ) = h( k+1, j ) - sum*v2
315 h( k+2, j ) = h( k+2, j ) - sum*v3
321 DO 70 j = i1,
min( k+3, i )
322 sum = t1*h( j, k ) + t2*h( j, k+1 ) + t3*h( j, k+2 )
323 h( j, k ) = h( j, k ) - sum
324 h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
325 h( j, k+2 ) = h( j, k+2 ) - sum*conjg( v3 )
333 sum = t1*z( j, k ) + t2*z( j, k+1 ) +
335 z( j, k ) = z( j, k ) - sum
336 z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
337 z( j, k+2 ) = z( j, k+2 ) - sum*conjg( v3 )
340 ELSE IF( nr.EQ.2 )
THEN
346 sum = conjg( t1 )*h( k, j ) +
347 $ conjg( t2 )*h( k+1, j )
348 h( k, j ) = h( k, j ) - sum
349 h( k+1, j ) = h( k+1, j ) - sum*v2
355 DO 100 j = i1,
min( k+2, i )
356 sum = t1*h( j, k ) + t2*h( j, k+1 )
357 h( j, k ) = h( j, k ) - sum
358 h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
365 DO 110 j = iloz, ihiz
366 sum = t1*z( j, k ) + t2*z( j, k+1 )
367 z( j, k ) = z( j, k ) - sum
368 z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
402 ELSE IF( l.EQ.i-1 )
THEN
409 CALL clanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
410 $ h( i, i ), w( i-1 ), w( i ), cs, sn )
417 $
CALL crot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
419 CALL crot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
426 CALL crot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,