160 SUBROUTINE dsytri_3x( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
169 INTEGER INFO, LDA, N, NB
173 DOUBLE PRECISION A( lda, * ), E( * ), WORK( n+nb+1, * )
179 DOUBLE PRECISION ONE, ZERO
180 parameter ( one = 1.0d+0, zero = 0.0d+0 )
184 INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
185 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, U01_I_J, U01_IP1_J,
196 INTRINSIC abs, max, mod
203 upper = lsame( uplo,
'U' )
204 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
206 ELSE IF( n.LT.0 )
THEN
208 ELSE IF( lda.LT.max( 1, n ) )
THEN
215 CALL xerbla(
'DSYTRI_3X', -info )
224 work( k, 1 ) = e( k )
234 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
242 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
268 CALL dtrtri( uplo,
'U', n, a, lda, info )
274 IF( ipiv( k ).GT.0 )
THEN
276 work( k, invd ) = one / a( k, k )
277 work( k, invd+1 ) = zero
282 akp1 = a( k+1, k+1 ) / t
283 akkp1 = work( k+1, 1 ) / t
284 d = t*( ak*akp1-one )
285 work( k, invd ) = akp1 / d
286 work( k+1, invd+1 ) = ak / d
287 work( k, invd+1 ) = -akkp1 / d
288 work( k+1, invd ) = work( k, invd+1 )
301 IF( cut.LE.nnb )
THEN
306 DO i = cut+1-nnb, cut
307 IF( ipiv( i ).LT.0 ) icount = icount + 1
310 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
319 work( i, j ) = a( i, cut+j )
326 work( u11+i, i ) = one
328 work( u11+i, j ) = zero
331 work( u11+i, j ) = a( cut+i, cut+j )
339 IF( ipiv( i ).GT.0 )
THEN
341 work( i, j ) = work( i, invd ) * work( i, j )
345 u01_i_j = work( i, j )
346 u01_ip1_j = work( i+1, j )
347 work( i, j ) = work( i, invd ) * u01_i_j
348 $ + work( i, invd+1 ) * u01_ip1_j
349 work( i+1, j ) = work( i+1, invd ) * u01_i_j
350 $ + work( i+1, invd+1 ) * u01_ip1_j
360 DO WHILE ( i.LE.nnb )
361 IF( ipiv( cut+i ).GT.0 )
THEN
363 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
367 u11_i_j = work(u11+i,j)
368 u11_ip1_j = work(u11+i+1,j)
369 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
370 $ + work(cut+i,invd+1) * work(u11+i+1,j)
371 work( u11+i+1, j ) = work(cut+i+1,invd) * u11_i_j
372 $ + work(cut+i+1,invd+1) * u11_ip1_j
381 CALL dtrmm(
'L',
'U',
'T',
'U', nnb, nnb,
382 $ one, a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
387 a( cut+i, cut+j ) = work( u11+i, j )
393 CALL dgemm(
'T',
'N', nnb, nnb, cut, one, a( 1, cut+1 ),
394 $ lda, work, n+nb+1, zero, work(u11+1,1), n+nb+1 )
401 a( cut+i, cut+j ) = a( cut+i, cut+j ) + work(u11+i,j)
407 CALL dtrmm(
'L', uplo,
'T',
'U', cut, nnb,
408 $ one, a, lda, work, n+nb+1 )
415 a( i, cut+j ) = work( i, j )
435 ip = abs( ipiv( i ) )
437 IF (i .LT. ip)
CALL dsyswapr( uplo, n, a, lda, i ,ip )
438 IF (i .GT. ip)
CALL dsyswapr( uplo, n, a, lda, ip ,i )
448 CALL dtrtri( uplo,
'U', n, a, lda, info )
453 DO WHILE ( k .GE. 1 )
454 IF( ipiv( k ).GT.0 )
THEN
456 work( k, invd ) = one / a( k, k )
457 work( k, invd+1 ) = zero
461 ak = a( k-1, k-1 ) / t
463 akkp1 = work( k-1, 1 ) / t
464 d = t*( ak*akp1-one )
465 work( k-1, invd ) = akp1 / d
466 work( k, invd ) = ak / d
467 work( k, invd+1 ) = -akkp1 / d
468 work( k-1, invd+1 ) = work( k, invd+1 )
481 IF( (cut + nnb).GT.n )
THEN
486 DO i = cut + 1, cut+nnb
487 IF ( ipiv( i ).LT.0 ) icount = icount + 1
490 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
497 work( i, j ) = a( cut+nnb+i, cut+j )
504 work( u11+i, i) = one
506 work( u11+i, j ) = zero
509 work( u11+i, j ) = a( cut+i, cut+j )
517 IF( ipiv( cut+nnb+i ).GT.0 )
THEN
519 work( i, j ) = work( cut+nnb+i, invd) * work( i, j)
524 u01_ip1_j = work(i-1,j)
525 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
526 $ work(cut+nnb+i,invd+1)*u01_ip1_j
527 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
528 $ work(cut+nnb+i-1,invd)*u01_ip1_j
539 IF( ipiv( cut+i ).GT.0 )
THEN
541 work( u11+i, j ) = work( cut+i, invd)*work(u11+i,j)
546 u11_i_j = work( u11+i, j )
547 u11_ip1_j = work( u11+i-1, j )
548 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
549 $ + work(cut+i,invd+1) * u11_ip1_j
550 work( u11+i-1, j ) = work(cut+i-1,invd+1) * u11_i_j
551 $ + work(cut+i-1,invd) * u11_ip1_j
560 CALL dtrmm(
'L', uplo,
'T',
'U', nnb, nnb, one,
561 $ a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
567 a( cut+i, cut+j ) = work( u11+i, j )
571 IF( (cut+nnb).LT.n )
THEN
575 CALL dgemm(
'T',
'N', nnb, nnb, n-nnb-cut, one,
576 $ a( cut+nnb+1, cut+1 ), lda, work, n+nb+1,
577 $ zero, work( u11+1, 1 ), n+nb+1 )
584 a( cut+i, cut+j ) = a( cut+i, cut+j )+work(u11+i,j)
590 CALL dtrmm(
'L', uplo,
'T',
'U', n-nnb-cut, nnb, one,
591 $ a( cut+nnb+1, cut+nnb+1 ), lda, work,
598 a( cut+nnb+i, cut+j ) = work( i, j )
608 a( cut+i, cut+j ) = work( u11+i, j )
631 ip = abs( ipiv( i ) )
633 IF (i .LT. ip)
CALL dsyswapr( uplo, n, a, lda, i ,ip )
634 IF (i .GT. ip)
CALL dsyswapr( uplo, n, a, lda, ip ,i )
subroutine dsytri_3x(UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO)
DSYTRI_3X
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dsyswapr(UPLO, N, A, LDA, I1, I2)
DSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dtrtri(UPLO, DIAG, N, A, LDA, INFO)
DTRTRI