160 SUBROUTINE zhetri_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, * )
180 parameter ( one = 1.0d+0 )
181 COMPLEX*16 CONE, CZERO
182 parameter ( cone = ( 1.0d+0, 0.0d+0 ),
183 $ czero = ( 0.0d+0, 0.0d+0 ) )
187 INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
188 DOUBLE PRECISION AK, AKP1, T
189 COMPLEX*16 AKKP1, D, U01_I_J, U01_IP1_J, U11_I_J,
200 INTRINSIC abs, dconjg, dble, max
207 upper = lsame( uplo,
'U' )
208 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
210 ELSE IF( n.LT.0 )
THEN
212 ELSE IF( lda.LT.max( 1, n ) )
THEN
219 CALL xerbla(
'ZHETRI_3X', -info )
228 work( k, 1 ) = e( k )
238 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
246 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.czero )
272 CALL ztrtri( uplo,
'U', n, a, lda, info )
278 IF( ipiv( k ).GT.0 )
THEN
280 work( k, invd ) = one / dble( a( k, k ) )
281 work( k, invd+1 ) = czero
284 t = abs( work( k+1, 1 ) )
285 ak = dble( a( k, k ) ) / t
286 akp1 = dble( a( k+1, k+1 ) ) / t
287 akkp1 = work( k+1, 1 ) / t
288 d = t*( ak*akp1-cone )
289 work( k, invd ) = akp1 / d
290 work( k+1, invd+1 ) = ak / d
291 work( k, invd+1 ) = -akkp1 / d
292 work( k+1, invd ) = dconjg( work( k, invd+1 ) )
305 IF( cut.LE.nnb )
THEN
310 DO i = cut+1-nnb, cut
311 IF( ipiv( i ).LT.0 ) icount = icount + 1
314 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
323 work( i, j ) = a( i, cut+j )
330 work( u11+i, i ) = cone
332 work( u11+i, j ) = czero
335 work( u11+i, j ) = a( cut+i, cut+j )
343 IF( ipiv( i ).GT.0 )
THEN
345 work( i, j ) = work( i, invd ) * work( i, j )
349 u01_i_j = work( i, j )
350 u01_ip1_j = work( i+1, j )
351 work( i, j ) = work( i, invd ) * u01_i_j
352 $ + work( i, invd+1 ) * u01_ip1_j
353 work( i+1, j ) = work( i+1, invd ) * u01_i_j
354 $ + work( i+1, invd+1 ) * u01_ip1_j
364 DO WHILE ( i.LE.nnb )
365 IF( ipiv( cut+i ).GT.0 )
THEN
367 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
371 u11_i_j = work(u11+i,j)
372 u11_ip1_j = work(u11+i+1,j)
373 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
374 $ + work(cut+i,invd+1) * work(u11+i+1,j)
375 work( u11+i+1, j ) = work(cut+i+1,invd) * u11_i_j
376 $ + work(cut+i+1,invd+1) * u11_ip1_j
385 CALL ztrmm(
'L',
'U',
'C',
'U', nnb, nnb,
386 $ cone, a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
391 a( cut+i, cut+j ) = work( u11+i, j )
397 CALL zgemm(
'C',
'N', nnb, nnb, cut, cone, a( 1, cut+1 ),
398 $ lda, work, n+nb+1, czero, work(u11+1,1),
406 a( cut+i, cut+j ) = a( cut+i, cut+j ) + work(u11+i,j)
412 CALL ztrmm(
'L', uplo,
'C',
'U', cut, nnb,
413 $ cone, a, lda, work, n+nb+1 )
420 a( i, cut+j ) = work( i, j )
440 ip = abs( ipiv( i ) )
442 IF (i .LT. ip)
CALL zheswapr( uplo, n, a, lda, i ,ip )
443 IF (i .GT. ip)
CALL zheswapr( uplo, n, a, lda, ip ,i )
453 CALL ztrtri( uplo,
'U', n, a, lda, info )
458 DO WHILE ( k .GE. 1 )
459 IF( ipiv( k ).GT.0 )
THEN
461 work( k, invd ) = one / dble( a( k, k ) )
462 work( k, invd+1 ) = czero
465 t = abs( work( k-1, 1 ) )
466 ak = dble( a( k-1, k-1 ) ) / t
467 akp1 = dble( a( k, k ) ) / t
468 akkp1 = work( k-1, 1 ) / t
469 d = t*( ak*akp1-cone )
470 work( k-1, invd ) = akp1 / d
471 work( k, invd ) = ak / d
472 work( k, invd+1 ) = -akkp1 / d
473 work( k-1, invd+1 ) = dconjg( work( k, invd+1 ) )
486 IF( (cut + nnb).GT.n )
THEN
491 DO i = cut + 1, cut+nnb
492 IF ( ipiv( i ).LT.0 ) icount = icount + 1
495 IF( mod( icount, 2 ).EQ.1 ) nnb = nnb + 1
502 work( i, j ) = a( cut+nnb+i, cut+j )
509 work( u11+i, i) = cone
511 work( u11+i, j ) = czero
514 work( u11+i, j ) = a( cut+i, cut+j )
522 IF( ipiv( cut+nnb+i ).GT.0 )
THEN
524 work( i, j ) = work( cut+nnb+i, invd) * work( i, j)
529 u01_ip1_j = work(i-1,j)
530 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
531 $ work(cut+nnb+i,invd+1)*u01_ip1_j
532 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
533 $ work(cut+nnb+i-1,invd)*u01_ip1_j
544 IF( ipiv( cut+i ).GT.0 )
THEN
546 work( u11+i, j ) = work( cut+i, invd)*work(u11+i,j)
551 u11_i_j = work( u11+i, j )
552 u11_ip1_j = work( u11+i-1, j )
553 work( u11+i, j ) = work(cut+i,invd) * work(u11+i,j)
554 $ + work(cut+i,invd+1) * u11_ip1_j
555 work( u11+i-1, j ) = work(cut+i-1,invd+1) * u11_i_j
556 $ + work(cut+i-1,invd) * u11_ip1_j
565 CALL ztrmm(
'L', uplo,
'C',
'U', nnb, nnb, cone,
566 $ a( cut+1, cut+1 ), lda, work( u11+1, 1 ),
572 a( cut+i, cut+j ) = work( u11+i, j )
576 IF( (cut+nnb).LT.n )
THEN
580 CALL zgemm(
'C',
'N', nnb, nnb, n-nnb-cut, cone,
581 $ a( cut+nnb+1, cut+1 ), lda, work, n+nb+1,
582 $ czero, work( u11+1, 1 ), n+nb+1 )
589 a( cut+i, cut+j ) = a( cut+i, cut+j )+work(u11+i,j)
595 CALL ztrmm(
'L', uplo,
'C',
'U', n-nnb-cut, nnb, cone,
596 $ a( cut+nnb+1, cut+nnb+1 ), lda, work,
603 a( cut+nnb+i, cut+j ) = work( i, j )
613 a( cut+i, cut+j ) = work( u11+i, j )
636 ip = abs( ipiv( i ) )
638 IF (i .LT. ip)
CALL zheswapr( uplo, n, a, lda, i ,ip )
639 IF (i .GT. ip)
CALL zheswapr( uplo, n, a, lda, ip ,i )
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 zheswapr(UPLO, N, A, LDA, I1, I2)
ZHESWAPR applies an elementary permutation on the rows and columns of a Hermitian matrix...
subroutine zhetri_3x(UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO)
ZHETRI_3X