121 SUBROUTINE zsytri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
130 INTEGER INFO, LDA, N, NB
134 COMPLEX*16 A( lda, * ), WORK( n+nb+1,* )
141 parameter ( one = ( 1.0d+0, 0.0d+0 ),
142 $ zero = ( 0.0d+0, 0.0d+0 ) )
146 INTEGER I, IINFO, IP, K, CUT, NNB
150 COMPLEX*16 AK, AKKP1, AKP1, D, T
151 COMPLEX*16 U01_I_J, U01_IP1_J
152 COMPLEX*16 U11_I_J, U11_IP1_J
170 upper = lsame( uplo,
'U' )
171 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
173 ELSE IF( n.LT.0 )
THEN
175 ELSE IF( lda.LT.max( 1, n ) )
THEN
183 CALL xerbla(
'ZSYTRI2X', -info )
192 CALL zsyconv( uplo,
'C', n, a, lda, ipiv, work, iinfo )
201 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
209 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
229 CALL ztrtri( uplo,
'U', n, a, lda, info )
234 DO WHILE ( k .LE. n )
235 IF( ipiv( k ).GT.0 )
THEN
237 work(k,invd) = 1/ a( k, k )
244 akp1 = a( k+1, k+1 ) / t
245 akkp1 = work(k+1,1) / t
246 d = t*( ak*akp1-one )
247 work(k,invd) = akp1 / d
248 work(k+1,invd+1) = ak / d
249 work(k,invd+1) = -akkp1 / d
250 work(k+1,invd) = -akkp1 / d
260 DO WHILE (cut .GT. 0)
262 IF (cut .LE. nnb)
THEN
268 IF (ipiv(i) .LT. 0) count=count+1
271 IF (mod(count,2) .EQ. 1) nnb=nnb+1
292 work(u11+i,j)=a(cut+i,cut+j)
299 DO WHILE (i .LE. cut)
300 IF (ipiv(i) > 0)
THEN
302 work(i,j)=work(i,invd)*work(i,j)
308 u01_ip1_j = work(i+1,j)
309 work(i,j)=work(i,invd)*u01_i_j+
310 $ work(i,invd+1)*u01_ip1_j
311 work(i+1,j)=work(i+1,invd)*u01_i_j+
312 $ work(i+1,invd+1)*u01_ip1_j
321 DO WHILE (i .LE. nnb)
322 IF (ipiv(cut+i) > 0)
THEN
324 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
329 u11_i_j = work(u11+i,j)
330 u11_ip1_j = work(u11+i+1,j)
331 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
332 $ work(cut+i,invd+1)*work(u11+i+1,j)
333 work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
334 $ work(cut+i+1,invd+1)*u11_ip1_j
342 CALL ztrmm(
'L',
'U',
'T',
'U',nnb, nnb,
343 $ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
347 a(cut+i,cut+j)=work(u11+i,j)
353 CALL zgemm(
'T',
'N',nnb,nnb,cut,one,a(1,cut+1),lda,
354 $ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
360 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
366 CALL ztrmm(
'L',uplo,
'T',
'U',cut, nnb,
367 $ one,a,lda,work,n+nb+1)
385 DO WHILE ( i .LE. n )
386 IF( ipiv(i) .GT. 0 )
THEN
388 IF (i .LT. ip)
CALL zsyswapr( uplo, n, a, lda, i ,ip )
389 IF (i .GT. ip)
CALL zsyswapr( uplo, n, a, lda, ip ,i )
394 $
CALL zsyswapr( uplo, n, a, lda, i-1 ,ip )
396 $
CALL zsyswapr( uplo, n, a, lda, ip ,i-1 )
406 CALL ztrtri( uplo,
'U', n, a, lda, info )
411 DO WHILE ( k .GE. 1 )
412 IF( ipiv( k ).GT.0 )
THEN
414 work(k,invd) = 1/ a( k, k )
420 ak = a( k-1, k-1 ) / t
422 akkp1 = work(k-1,1) / t
423 d = t*( ak*akp1-one )
424 work(k-1,invd) = akp1 / d
425 work(k,invd) = ak / d
426 work(k,invd+1) = -akkp1 / d
427 work(k-1,invd+1) = -akkp1 / d
437 DO WHILE (cut .LT. n)
439 IF (cut + nnb .GE. n)
THEN
445 IF (ipiv(i) .LT. 0) count=count+1
448 IF (mod(count,2) .EQ. 1) nnb=nnb+1
453 work(i,j)=a(cut+nnb+i,cut+j)
463 work(u11+i,j)=a(cut+i,cut+j)
471 IF (ipiv(cut+nnb+i) > 0)
THEN
473 work(i,j)=work(cut+nnb+i,invd)*work(i,j)
479 u01_ip1_j = work(i-1,j)
480 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
481 $ work(cut+nnb+i,invd+1)*u01_ip1_j
482 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
483 $ work(cut+nnb+i-1,invd)*u01_ip1_j
493 IF (ipiv(cut+i) > 0)
THEN
495 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
500 u11_i_j = work(u11+i,j)
501 u11_ip1_j = work(u11+i-1,j)
502 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
503 $ work(cut+i,invd+1)*u11_ip1_j
504 work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
505 $ work(cut+i-1,invd)*u11_ip1_j
513 CALL ztrmm(
'L',uplo,
'T',
'U',nnb, nnb,
514 $ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
518 a(cut+i,cut+j)=work(u11+i,j)
523 IF ( (cut+nnb) .LT. n )
THEN
527 CALL zgemm(
'T',
'N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
528 $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
535 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
541 CALL ztrmm(
'L',uplo,
'T',
'U', n-nnb-cut, nnb,
542 $ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
547 a(cut+nnb+i,cut+j)=work(i,j)
556 a(cut+i,cut+j)=work(u11+i,j)
569 DO WHILE ( i .GE. 1 )
570 IF( ipiv(i) .GT. 0 )
THEN
572 IF (i .LT. ip)
CALL zsyswapr( uplo, n, a, lda, i ,ip )
573 IF (i .GT. ip)
CALL zsyswapr( uplo, n, a, lda, ip ,i )
576 IF ( i .LT. ip)
CALL zsyswapr( uplo, n, a, lda, i ,ip )
577 IF ( i .GT. ip)
CALL zsyswapr( uplo, n, a, lda, ip ,i )
subroutine zsyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
ZSYCONV
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zsytri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
ZSYTRI2X
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