160 SUBROUTINE zsytri_3x( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
169 INTEGER INFO, LDA, N, NB
173 COMPLEX*16 A( lda, * ), E( * ), WORK( n+nb+1, * )
179 COMPLEX*16 CONE, CZERO
180 parameter ( cone = ( 1.0d+0, 0.0d+0 ),
181 $ czero = ( 0.0d+0, 0.0d+0 ) )
185 INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
186 COMPLEX*16 AK, AKKP1, AKP1, D, T, U01_I_J, U01_IP1_J,
197 INTRINSIC abs, max, mod
204 upper = lsame( uplo,
'U' )
205 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
207 ELSE IF( n.LT.0 )
THEN
209 ELSE IF( lda.LT.max( 1, n ) )
THEN
216 CALL xerbla(
'ZSYTRI_3X', -info )
225 work( k, 1 ) = e( k )
235 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
243 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
269 CALL ztrtri( uplo,
'U', n, a, lda, info )
275 IF( ipiv( k ).GT.0 )
THEN
277 work( k, invd ) = cone / a( k, k )
278 work( k, invd+1 ) = czero
283 akp1 = a( k+1, k+1 ) / t
284 akkp1 = work( k+1, 1 ) / t
285 d = t*( ak*akp1-cone )
286 work( k, invd ) = akp1 / d
287 work( k+1, invd+1 ) = ak / d
288 work( k, invd+1 ) = -akkp1 / d
289 work( k+1, invd ) = work( k, invd+1 )
302 IF( cut.LE.nnb )
THEN
307 DO i = cut+1-nnb, cut
308 IF( ipiv( i ).LT.0 ) icount = icount + 1
311 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
320 work( i, j ) = a( i, cut+j )
327 work( u11+i, i ) = cone
329 work( u11+i, j ) = czero
332 work( u11+i, j ) = a( cut+i, cut+j )
340 IF( ipiv( i ).GT.0 )
THEN
342 work( i, j ) = work( i, invd ) * work( i, j )
346 u01_i_j = work( i, j )
347 u01_ip1_j = work( i+1, j )
348 work( i, j ) = work( i, invd ) * u01_i_j
349 $ + work( i, invd+1 ) * u01_ip1_j
350 work( i+1, j ) = work( i+1, invd ) * u01_i_j
351 $ + work( i+1, invd+1 ) * u01_ip1_j
361 DO WHILE ( i.LE.nnb )
362 IF( ipiv( cut+i ).GT.0 )
THEN
364 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
368 u11_i_j = work(u11+i,j)
369 u11_ip1_j = work(u11+i+1,j)
370 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
371 $ + work(cut+i,invd+1) * work(u11+i+1,j)
372 work( u11+i+1, j ) = work(cut+i+1,invd) * u11_i_j
373 $ + work(cut+i+1,invd+1) * u11_ip1_j
382 CALL ztrmm(
'L',
'U',
'T',
'U', nnb, nnb,
383 $ cone, a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
388 a( cut+i, cut+j ) = work( u11+i, j )
394 CALL zgemm(
'T',
'N', nnb, nnb, cut, cone, a( 1, cut+1 ),
395 $ lda, work, n+nb+1, czero, work(u11+1,1),
403 a( cut+i, cut+j ) = a( cut+i, cut+j ) + work(u11+i,j)
409 CALL ztrmm(
'L', uplo,
'T',
'U', cut, nnb,
410 $ cone, a, lda, work, n+nb+1 )
417 a( i, cut+j ) = work( i, j )
437 ip = abs( ipiv( i ) )
439 IF (i .LT. ip)
CALL zsyswapr( uplo, n, a, lda, i ,ip )
440 IF (i .GT. ip)
CALL zsyswapr( uplo, n, a, lda, ip ,i )
450 CALL ztrtri( uplo,
'U', n, a, lda, info )
455 DO WHILE ( k .GE. 1 )
456 IF( ipiv( k ).GT.0 )
THEN
458 work( k, invd ) = cone / a( k, k )
459 work( k, invd+1 ) = czero
463 ak = a( k-1, k-1 ) / t
465 akkp1 = work( k-1, 1 ) / t
466 d = t*( ak*akp1-cone )
467 work( k-1, invd ) = akp1 / d
468 work( k, invd ) = ak / d
469 work( k, invd+1 ) = -akkp1 / d
470 work( k-1, invd+1 ) = work( k, invd+1 )
483 IF( (cut + nnb).GT.n )
THEN
488 DO i = cut + 1, cut+nnb
489 IF ( ipiv( i ).LT.0 ) icount = icount + 1
492 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
499 work( i, j ) = a( cut+nnb+i, cut+j )
506 work( u11+i, i) = cone
508 work( u11+i, j ) = czero
511 work( u11+i, j ) = a( cut+i, cut+j )
519 IF( ipiv( cut+nnb+i ).GT.0 )
THEN
521 work( i, j ) = work( cut+nnb+i, invd) * work( i, j)
526 u01_ip1_j = work(i-1,j)
527 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
528 $ work(cut+nnb+i,invd+1)*u01_ip1_j
529 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
530 $ work(cut+nnb+i-1,invd)*u01_ip1_j
541 IF( ipiv( cut+i ).GT.0 )
THEN
543 work( u11+i, j ) = work( cut+i, invd)*work(u11+i,j)
548 u11_i_j = work( u11+i, j )
549 u11_ip1_j = work( u11+i-1, j )
550 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
551 $ + work(cut+i,invd+1) * u11_ip1_j
552 work( u11+i-1, j ) = work(cut+i-1,invd+1) * u11_i_j
553 $ + work(cut+i-1,invd) * u11_ip1_j
562 CALL ztrmm(
'L', uplo,
'T',
'U', nnb, nnb, cone,
563 $ a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
569 a( cut+i, cut+j ) = work( u11+i, j )
573 IF( (cut+nnb).LT.n )
THEN
577 CALL zgemm(
'T',
'N', nnb, nnb, n-nnb-cut, cone,
578 $ a( cut+nnb+1, cut+1 ), lda, work, n+nb+1,
579 $ czero, work( u11+1, 1 ), n+nb+1 )
586 a( cut+i, cut+j ) = a( cut+i, cut+j )+work(u11+i,j)
592 CALL ztrmm(
'L', uplo,
'T',
'U', n-nnb-cut, nnb, cone,
593 $ a( cut+nnb+1, cut+nnb+1 ), lda, work,
600 a( cut+nnb+i, cut+j ) = work( i, j )
610 a( cut+i, cut+j ) = work( u11+i, j )
633 ip = abs( ipiv( i ) )
635 IF (i .LT. ip)
CALL zsyswapr( uplo, n, a, lda, i ,ip )
636 IF (i .GT. ip)
CALL zsyswapr( uplo, n, a, lda, ip ,i )
subroutine zsytri_3x(UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO)
ZSYTRI_3X
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ztrtri(UPLO, DIAG, N, A, LDA, INFO)
ZTRTRI
subroutine zsyswapr(UPLO, N, A, LDA, I1, I2)
ZSYSWAPR