160 SUBROUTINE chetri_3x( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
169 INTEGER INFO, LDA, N, NB
173 COMPLEX A( lda, * ), E( * ), WORK( n+nb+1, * )
180 parameter ( one = 1.0e+0 )
182 parameter ( cone = ( 1.0e+0, 0.0e+0 ),
183 $ czero = ( 0.0e+0, 0.0e+0 ) )
187 INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
189 COMPLEX AKKP1, D, U01_I_J, U01_IP1_J, U11_I_J,
200 INTRINSIC abs, conjg, max, real
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(
'CHETRI_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 ctrtri( uplo,
'U', n, a, lda, info )
278 IF( ipiv( k ).GT.0 )
THEN
280 work( k, invd ) = one /
REAL( A( K, K ) )
281 work( k, invd+1 ) = czero
284 t = abs( work( k+1, 1 ) )
285 ak =
REAL( A( K, K ) ) / T
286 akp1 =
REAL( 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 ) = conjg( 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 ctrmm(
'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 cgemm(
'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 ctrmm(
'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 cheswapr( uplo, n, a, lda, i ,ip )
443 IF (i .GT. ip)
CALL cheswapr( uplo, n, a, lda, ip ,i )
453 CALL ctrtri( uplo,
'U', n, a, lda, info )
458 DO WHILE ( k .GE. 1 )
459 IF( ipiv( k ).GT.0 )
THEN
461 work( k, invd ) = one /
REAL( A( K, K ) )
462 work( k, invd+1 ) = czero
465 t = abs( work( k-1, 1 ) )
466 ak =
REAL( A( K-1, K-1 ) ) / T
467 akp1 =
REAL( 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 ) = conjg( 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 ctrmm(
'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 cgemm(
'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 ctrmm(
'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 cheswapr( uplo, n, a, lda, i ,ip )
639 IF (i .GT. ip)
CALL cheswapr( uplo, n, a, lda, ip ,i )
subroutine cheswapr(UPLO, N, A, LDA, I1, I2)
CHESWAPR applies an elementary permutation on the rows and columns of a Hermitian matrix...
subroutine chetri_3x(UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO)
CHETRI_3X
subroutine xerbla(SRNAME, INFO)
XERBLA
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