195 SUBROUTINE clahqr( 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 h( ldh, * ), w( * ), z( ldz, * )
215 parameter( itmax = 30 )
217 parameter( zero = ( 0.0e0, 0.0e0 ),
218 $ one = ( 1.0e0, 0.0e0 ) )
219 REAL rzero, rone, half
220 parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
222 parameter( dat1 = 3.0e0 / 4.0e0 )
225 COMPLEX cdum, h11, h11s, h22, sc, sum, t, t1, temp, u,
227 REAL aa, ab, ba, bb, h10, h21, rtemp, s, safmax,
228 $ safmin, smlnum, sx, t2, tst, ulp
229 INTEGER i, i1, i2, its,
j, jhi, jlo, k, l, m, nh, nz
246 INTRINSIC abs, aimag, conjg, max, min,
REAL, sqrt
249 cabs1( cdum ) = abs(
REAL( CDUM ) ) + abs( aimag( cdum ) )
259 IF( ilo.EQ.ihi )
THEN
260 w( ilo ) = h( ilo, ilo )
265 DO 10
j = ilo, ihi - 3
270 $ h( ihi, ihi-2 ) = zero
279 DO 20 i = ilo + 1, ihi
280 IF( aimag( h( i, i-1 ) ).NE.rzero )
THEN
284 sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
285 sc = conjg( sc ) / abs( sc )
286 h( i, i-1 ) = abs( h( i, i-1 ) )
287 CALL
cscal( jhi-i+1, sc, h( i, i ), ldh )
288 CALL
cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),
291 $ CALL
cscal( ihiz-iloz+1, conjg( sc ), z( iloz, i ), 1 )
300 safmin =
slamch(
'SAFE MINIMUM' )
301 safmax = rone / safmin
302 CALL
slabad( safmin, safmax )
303 ulp =
slamch(
'PRECISION' )
304 smlnum = safmin*(
REAL( NH ) / ulp )
331 DO 130 its = 0, itmax
335 DO 40 k = i, l + 1, -1
336 IF( cabs1( h( k, k-1 ) ).LE.smlnum )
338 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
339 IF( tst.EQ.zero )
THEN
341 $ tst = tst + abs(
REAL( H( K-1, K-2 ) ) )
343 $ tst = tst + abs(
REAL( H( K+1, K ) ) )
349 IF( abs(
REAL( H( K, K-1 ) ) ).LE.ulp*tst ) then
350 ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
351 ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
352 aa = max( cabs1( h( k, k ) ),
353 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
354 bb = min( cabs1( h( k, k ) ),
355 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
357 IF( ba*( ab / s ).LE.max( smlnum,
358 $ ulp*( bb*( aa / s ) ) ) )go to 50
379 IF( .NOT.wantt )
THEN
388 s = dat1*abs(
REAL( H( L+1, L ) ) )
390 ELSE IF( its.EQ.20 )
THEN
394 s = dat1*abs(
REAL( H( I, I-1 ) ) )
401 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
403 IF( s.NE.rzero )
THEN
404 x = half*( h( i-1, i-1 )-t )
406 s = max( s, cabs1( x ) )
407 y = s*sqrt( ( x / s )**2+( u / s )**2 )
408 IF( sx.GT.rzero )
THEN
409 IF(
REAL( x / sx )*
REAL( y )+aimag( x / sx )*
410 $ aimag( y ).LT.rzero )y = -y
412 t = t - u*
cladiv( u, ( x+y ) )
418 DO 60 m = i - 1, l + 1, -1
427 h21 =
REAL( H( M+1, M ) )
428 s = cabs1( h11s ) + abs( h21 )
433 h10 =
REAL( H( M, M-1 ) )
434 IF( abs( h10 )*abs( h21 ).LE.ulp*
435 $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
441 h21 =
REAL( H( L+1, L ) )
442 s = cabs1( h11s ) + abs( h21 )
466 $ CALL
ccopy( 2, h( k, k-1 ), 1, v, 1 )
467 CALL
clarfg( 2, v( 1 ), v( 2 ), 1, t1 )
479 sum = conjg( t1 )*h( k,
j ) + t2*h( k+1,
j )
480 h( k,
j ) = h( k,
j ) - sum
481 h( k+1,
j ) = h( k+1,
j ) - sum*v2
487 DO 90
j = i1, min( k+2, i )
488 sum = t1*h(
j, k ) + t2*h(
j, k+1 )
489 h(
j, k ) = h(
j, k ) - sum
490 h(
j, k+1 ) = h(
j, k+1 ) - sum*conjg( v2 )
497 DO 100
j = iloz, ihiz
498 sum = t1*z(
j, k ) + t2*z(
j, k+1 )
499 z(
j, k ) = z(
j, k ) - sum
500 z(
j, k+1 ) = z(
j, k+1 ) - sum*conjg( v2 )
504 IF( k.EQ.m .AND. m.GT.l )
THEN
512 temp = temp / abs( temp )
513 h( m+1, m ) = h( m+1, m )*conjg( temp )
515 $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
519 $ CALL
cscal( i2-
j, temp, h(
j,
j+1 ), ldh )
520 CALL
cscal(
j-i1, conjg( temp ), h( i1,
j ), 1 )
522 CALL
cscal( nz, conjg( temp ), z( iloz,
j ), 1 )
532 IF( aimag( temp ).NE.rzero )
THEN
537 $ CALL
cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
538 CALL
cscal( i-i1, temp, h( i1, i ), 1 )
540 CALL
cscal( nz, temp, z( iloz, i ), 1 )
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
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.