121 SUBROUTINE zhetri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
130 INTEGER INFO, LDA, N, NB
134 COMPLEX*16 A( lda, * ), WORK( n+nb+1,* )
141 COMPLEX*16 CONE, ZERO
142 parameter ( one = 1.0d+0,
143 $ cone = ( 1.0d+0, 0.0d+0 ),
144 $ zero = ( 0.0d+0, 0.0d+0 ) )
148 INTEGER I, IINFO, IP, K, CUT, NNB
152 COMPLEX*16 AK, AKKP1, AKP1, D, T
153 COMPLEX*16 U01_I_J, U01_IP1_J
154 COMPLEX*16 U11_I_J, U11_IP1_J
172 upper = lsame( uplo,
'U' )
173 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
175 ELSE IF( n.LT.0 )
THEN
177 ELSE IF( lda.LT.max( 1, n ) )
THEN
185 CALL xerbla(
'ZHETRI2X', -info )
194 CALL zsyconv( uplo,
'C', n, a, lda, ipiv, work, iinfo )
203 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
211 IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
231 CALL ztrtri( uplo,
'U', n, a, lda, info )
236 DO WHILE ( k .LE. n )
237 IF( ipiv( k ).GT.0 )
THEN
239 work(k,invd) = one /
REAL ( A( K, K ) )
244 t = abs( work(k+1,1) )
245 ak =
REAL ( A( K, K ) ) / T
246 akp1 =
REAL ( A( K+1, K+1 ) ) / T
247 akkp1 = work(k+1,1) / t
248 d = t*( ak*akp1-one )
249 work(k,invd) = akp1 / d
250 work(k+1,invd+1) = ak / d
251 work(k,invd+1) = -akkp1 / d
252 work(k+1,invd) = dconjg(work(k,invd+1) )
262 DO WHILE (cut .GT. 0)
264 IF (cut .LE. nnb)
THEN
270 IF (ipiv(i) .LT. 0) count=count+1
273 IF (mod(count,2) .EQ. 1) nnb=nnb+1
294 work(u11+i,j)=a(cut+i,cut+j)
301 DO WHILE (i .LE. cut)
302 IF (ipiv(i) > 0)
THEN
304 work(i,j)=work(i,invd)*work(i,j)
310 u01_ip1_j = work(i+1,j)
311 work(i,j)=work(i,invd)*u01_i_j+
312 $ work(i,invd+1)*u01_ip1_j
313 work(i+1,j)=work(i+1,invd)*u01_i_j+
314 $ work(i+1,invd+1)*u01_ip1_j
323 DO WHILE (i .LE. nnb)
324 IF (ipiv(cut+i) > 0)
THEN
326 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
331 u11_i_j = work(u11+i,j)
332 u11_ip1_j = work(u11+i+1,j)
333 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
334 $ work(cut+i,invd+1)*work(u11+i+1,j)
335 work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
336 $ work(cut+i+1,invd+1)*u11_ip1_j
344 CALL ztrmm(
'L',
'U',
'C',
'U',nnb, nnb,
345 $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
349 a(cut+i,cut+j)=work(u11+i,j)
355 CALL zgemm(
'C',
'N',nnb,nnb,cut,cone,a(1,cut+1),lda,
356 $ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
362 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
368 CALL ztrmm(
'L',uplo,
'C',
'U',cut, nnb,
369 $ cone,a,lda,work,n+nb+1)
387 DO WHILE ( i .LE. n )
388 IF( ipiv(i) .GT. 0 )
THEN
390 IF (i .LT. ip)
CALL zheswapr( uplo, n, a, lda, i ,ip )
391 IF (i .GT. ip)
CALL zheswapr( uplo, n, a, lda, ip ,i )
396 $
CALL zheswapr( uplo, n, a, lda, i-1 ,ip )
398 $
CALL zheswapr( uplo, n, a, lda, ip ,i-1 )
408 CALL ztrtri( uplo,
'U', n, a, lda, info )
413 DO WHILE ( k .GE. 1 )
414 IF( ipiv( k ).GT.0 )
THEN
416 work(k,invd) = one /
REAL ( A( K, K ) )
421 t = abs( work(k-1,1) )
422 ak =
REAL ( A( K-1, K-1 ) ) / T
423 akp1 =
REAL ( A( K, K ) ) / T
424 akkp1 = work(k-1,1) / t
425 d = t*( ak*akp1-one )
426 work(k-1,invd) = akp1 / d
427 work(k,invd) = ak / d
428 work(k,invd+1) = -akkp1 / d
429 work(k-1,invd+1) = dconjg(work(k,invd+1) )
439 DO WHILE (cut .LT. n)
441 IF (cut + nnb .GE. n)
THEN
447 IF (ipiv(i) .LT. 0) count=count+1
450 IF (mod(count,2) .EQ. 1) nnb=nnb+1
455 work(i,j)=a(cut+nnb+i,cut+j)
465 work(u11+i,j)=a(cut+i,cut+j)
473 IF (ipiv(cut+nnb+i) > 0)
THEN
475 work(i,j)=work(cut+nnb+i,invd)*work(i,j)
481 u01_ip1_j = work(i-1,j)
482 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
483 $ work(cut+nnb+i,invd+1)*u01_ip1_j
484 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
485 $ work(cut+nnb+i-1,invd)*u01_ip1_j
495 IF (ipiv(cut+i) > 0)
THEN
497 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
502 u11_i_j = work(u11+i,j)
503 u11_ip1_j = work(u11+i-1,j)
504 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
505 $ work(cut+i,invd+1)*u11_ip1_j
506 work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
507 $ work(cut+i-1,invd)*u11_ip1_j
515 CALL ztrmm(
'L',uplo,
'C',
'U',nnb, nnb,
516 $ cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
520 a(cut+i,cut+j)=work(u11+i,j)
524 IF ( (cut+nnb) .LT. n )
THEN
528 CALL zgemm(
'C',
'N',nnb,nnb,n-nnb-cut,cone,a(cut+nnb+1,cut+1)
529 $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
536 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
542 CALL ztrmm(
'L',uplo,
'C',
'U', n-nnb-cut, nnb,
543 $ cone,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
548 a(cut+nnb+i,cut+j)=work(i,j)
557 a(cut+i,cut+j)=work(u11+i,j)
570 DO WHILE ( i .GE. 1 )
571 IF( ipiv(i) .GT. 0 )
THEN
573 IF (i .LT. ip)
CALL zheswapr( uplo, n, a, lda, i ,ip )
574 IF (i .GT. ip)
CALL zheswapr( uplo, n, a, lda, ip ,i )
577 IF ( i .LT. ip)
CALL zheswapr( uplo, n, a, lda, i ,ip )
578 IF ( i .GT. ip)
CALL zheswapr( 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 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 zhetri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
ZHETRI2X