139 COMPLEX*16 A( LDA, * ), WORK( * )
146 COMPLEX*16 CONE, CZERO
147 parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ),
148 $ czero = ( 0.0d+0, 0.0d+0 ) )
152 INTEGER J, K, KP, KSTEP
153 DOUBLE PRECISION AK, AKP1, D, T
154 COMPLEX*16 AKKP1, TEMP
159 EXTERNAL lsame, zdotc
165 INTRINSIC abs, dconjg, max, dble
172 upper = lsame( uplo,
'U' )
173 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
175 ELSE IF( n.LT.0 )
THEN
177 ELSE IF( lda.LT.max( 1, n ) )
THEN
181 CALL xerbla(
'ZHETRI_ROOK', -info )
196 DO 10 info = n, 1, -1
197 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
205 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
226 IF( ipiv( k ).GT.0 )
THEN
232 a( k, k ) = one / dble( a( k, k ) )
237 CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
238 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
240 a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
250 t = abs( a( k, k+1 ) )
251 ak = dble( a( k, k ) ) / t
252 akp1 = dble( a( k+1, k+1 ) ) / t
253 akkp1 = a( k, k+1 ) / t
254 d = t*( ak*akp1-one )
256 a( k+1, k+1 ) = ak / d
257 a( k, k+1 ) = -akkp1 / d
262 CALL zcopy( k-1, a( 1, k ), 1, work, 1 )
263 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
265 a( k, k ) = a( k, k ) - dble( zdotc( k-1, work, 1, a( 1,
267 a( k, k+1 ) = a( k, k+1 ) -
268 $ zdotc( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
269 CALL zcopy( k-1, a( 1, k+1 ), 1, work, 1 )
270 CALL zhemv( uplo, k-1, -cone, a, lda, work, 1, czero,
272 a( k+1, k+1 ) = a( k+1, k+1 ) -
273 $ dble( zdotc( k-1, work, 1, a( 1, k+1 ),
279 IF( kstep.EQ.1 )
THEN
288 $
CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
290 DO 40 j = kp + 1, k - 1
291 temp = dconjg( a( j, k ) )
292 a( j, k ) = dconjg( a( kp, j ) )
296 a( kp, k ) = dconjg( a( kp, k ) )
299 a( k, k ) = a( kp, kp )
313 $
CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
315 DO 50 j = kp + 1, k - 1
316 temp = dconjg( a( j, k ) )
317 a( j, k ) = dconjg( a( kp, j ) )
321 a( kp, k ) = dconjg( a( kp, k ) )
324 a( k, k ) = a( kp, kp )
328 a( k, k+1 ) = a( kp, k+1 )
339 $
CALL zswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
341 DO 60 j = kp + 1, k - 1
342 temp = dconjg( a( j, k ) )
343 a( j, k ) = dconjg( a( kp, j ) )
347 a( kp, k ) = dconjg( a( kp, k ) )
350 a( k, k ) = a( kp, kp )
374 IF( ipiv( k ).GT.0 )
THEN
380 a( k, k ) = one / dble( a( k, k ) )
385 CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
386 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
387 $ 1, czero, a( k+1, k ), 1 )
388 a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
398 t = abs( a( k, k-1 ) )
399 ak = dble( a( k-1, k-1 ) ) / t
400 akp1 = dble( a( k, k ) ) / t
401 akkp1 = a( k, k-1 ) / t
402 d = t*( ak*akp1-one )
403 a( k-1, k-1 ) = akp1 / d
405 a( k, k-1 ) = -akkp1 / d
410 CALL zcopy( n-k, a( k+1, k ), 1, work, 1 )
411 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
412 $ 1, czero, a( k+1, k ), 1 )
413 a( k, k ) = a( k, k ) - dble( zdotc( n-k, work, 1,
415 a( k, k-1 ) = a( k, k-1 ) -
416 $ zdotc( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
418 CALL zcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
419 CALL zhemv( uplo, n-k, -cone, a( k+1, k+1 ), lda, work,
420 $ 1, czero, a( k+1, k-1 ), 1 )
421 a( k-1, k-1 ) = a( k-1, k-1 ) -
422 $ dble( zdotc( n-k, work, 1, a( k+1, k-1 ),
428 IF( kstep.EQ.1 )
THEN
437 $
CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
439 DO 90 j = k + 1, kp - 1
440 temp = dconjg( a( j, k ) )
441 a( j, k ) = dconjg( a( kp, j ) )
445 a( kp, k ) = dconjg( a( kp, k ) )
448 a( k, k ) = a( kp, kp )
462 $
CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
464 DO 100 j = k + 1, kp - 1
465 temp = dconjg( a( j, k ) )
466 a( j, k ) = dconjg( a( kp, j ) )
470 a( kp, k ) = dconjg( a( kp, k ) )
473 a( k, k ) = a( kp, kp )
477 a( k, k-1 ) = a( kp, k-1 )
488 $
CALL zswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
490 DO 110 j = k + 1, kp - 1
491 temp = dconjg( a( j, k ) )
492 a( j, k ) = dconjg( a( kp, j ) )
496 a( kp, k ) = dconjg( a( kp, k ) )
499 a( k, k ) = a( kp, kp )
subroutine xerbla(srname, info)
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine zhemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZHEMV
subroutine zhetri_rook(uplo, n, a, lda, ipiv, work, info)
ZHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch...
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP