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 )