121 SUBROUTINE chetri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
130 INTEGER INFO, LDA, N, NB
134 COMPLEX A( lda, * ), WORK( n+nb+1,* )
142 parameter ( one = 1.0e+0,
143 $ cone = ( 1.0e+0, 0.0e+0 ),
144 $ zero = ( 0.0e+0, 0.0e+0 ) )
148 INTEGER I, IINFO, IP, K, CUT, NNB
152 COMPLEX AK, AKKP1, AKP1, D, T
153 COMPLEX U01_I_J, U01_IP1_J
154 COMPLEX 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(
'CHETRI2X', -info )
194 CALL csyconv( 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 ctrtri( 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) = conjg(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 ctrmm(
'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 cgemm(
'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 ctrmm(
'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 cheswapr( uplo, n, a, lda, i ,ip )
391 IF (i .GT. ip)
CALL cheswapr( uplo, n, a, lda, ip ,i )
396 $
CALL cheswapr( uplo, n, a, lda, i-1 ,ip )
398 $
CALL cheswapr( uplo, n, a, lda, ip ,i-1 )
408 CALL ctrtri( 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) = conjg(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 ctrmm(
'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 cgemm(
'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 ctrmm(
'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 cheswapr( uplo, n, a, lda, i ,ip )
574 IF (i .GT. ip)
CALL cheswapr( uplo, n, a, lda, ip ,i )
577 IF ( i .LT. ip)
CALL cheswapr( uplo, n, a, lda, i ,ip )
578 IF ( i .GT. ip)
CALL cheswapr( uplo, n, a, lda, ip ,i )
subroutine chetri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
CHETRI2X
subroutine cheswapr(UPLO, N, A, LDA, I1, I2)
CHESWAPR applies an elementary permutation on the rows and columns of a Hermitian matrix...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine csyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
CSYCONV
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
subroutine ctrtri(UPLO, DIAG, N, A, LDA, INFO)
CTRTRI
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM