160 SUBROUTINE ssytri_3x( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
169 INTEGER INFO, LDA, N, NB
173 REAL A( lda, * ), E( * ), WORK( n+nb+1, * )
180 parameter ( one = 1.0e+0, zero = 0.0e+0 )
184 INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
185 REAL 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(
'SSYTRI_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 strtri( 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 strmm(
'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 sgemm(
'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 strmm(
'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 ssyswapr( uplo, n, a, lda, i ,ip )
438 IF (i .GT. ip)
CALL ssyswapr( uplo, n, a, lda, ip ,i )
448 CALL strtri( 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 strmm(
'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 sgemm(
'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 strmm(
'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 ssyswapr( uplo, n, a, lda, i ,ip )
634 IF (i .GT. ip)
CALL ssyswapr( uplo, n, a, lda, ip ,i )
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
subroutine strtri(UPLO, DIAG, N, A, LDA, INFO)
STRTRI
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ssyswapr(UPLO, N, A, LDA, I1, I2)
SSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix...
subroutine ssytri_3x(UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO)
SSYTRI_3X