subroutine zhetri2x	(	character	UPLO,
		integer	N,
		complex*16, dimension( lda, * )	A,
		integer	LDA,
		integer, dimension( * )	IPIV,
		complex*16, dimension( n+nb+1,* )	WORK,
		integer	NB,
		integer	INFO
	)

ZHETRI2X

Download ZHETRI2X + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
 A using the factorization A = U*D*U**H or A = L*D*L**H computed by
 ZHETRF.

Parameters

[in]	UPLO	UPLO is CHARACTER*1 Specifies whether the details of the factorization are stored as an upper or lower triangular matrix. = 'U': Upper triangular, form is A = U*D*U**H; = 'L': Lower triangular, form is A = L*D*L**H.
[in]	N	N is INTEGER The order of the matrix A. N >= 0.
[in,out]	A	A is COMPLEX*16 array, dimension (LDA,N) On entry, the NNB diagonal matrix D and the multipliers used to obtain the factor U or L as computed by ZHETRF. On exit, if INFO = 0, the (symmetric) inverse of the original matrix. If UPLO = 'U', the upper triangular part of the inverse is formed and the part of A below the diagonal is not referenced; if UPLO = 'L' the lower triangular part of the inverse is formed and the part of A above the diagonal is not referenced.
[in]	LDA	LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).
[in]	IPIV	IPIV is INTEGER array, dimension (N) Details of the interchanges and the NNB structure of D as determined by ZHETRF.
[out]	WORK	WORK is COMPLEX*16 array, dimension (N+NB+1,NB+3)
[in]	NB	NB is INTEGER Block size
[out]	INFO	INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its inverse could not be computed.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Date

November 2015

Definition at line 122 of file zhetri2x.f.

*
*  -- LAPACK computational routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      CHARACTER          uplo
      INTEGER            info, lda, n, nb
*     ..
*     .. Array Arguments ..
      INTEGER            ipiv( * )
      COMPLEX*16            a( lda, * ), work( n+nb+1,* )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               one
      COMPLEX*16            cone, zero
      parameter                ( one = 1.0d+0,
     $                   cone = ( 1.0d+0, 0.0d+0 ),
     $                   zero = ( 0.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            upper
      INTEGER            i, iinfo, ip, k, cut, nnb
      INTEGER            count
      INTEGER            j, u11, invd

      COMPLEX*16   ak, akkp1, akp1, d, t
      COMPLEX*16   u01_i_j, u01_ip1_j
      COMPLEX*16   u11_i_j, u11_ip1_j
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           zsyconv, xerbla, ztrtri
      EXTERNAL           zgemm, ztrmm, zheswapr
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      upper = lsame( uplo, 'U' )
      IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -4
      END IF
*
*     Quick return if possible
*
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZHETRI2X', -info )
         RETURN
      END IF
      IF( n.EQ.0 )
     $   RETURN
*
*     Convert A
*     Workspace got Non-diag elements of D
*
      CALL zsyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
*
*     Check that the diagonal matrix D is nonsingular.
*
      IF( upper ) THEN
*
*        Upper triangular storage: examine D from bottom to top
*
         DO info = n, 1, -1
            IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
     $         RETURN
         END DO
      ELSE
*
*        Lower triangular storage: examine D from top to bottom.
*
         DO info = 1, n
            IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
     $         RETURN
         END DO
      END IF
      info = 0
*
*  Splitting Workspace
*     U01 is a block (N,NB+1) 
*     The first element of U01 is in WORK(1,1)
*     U11 is a block (NB+1,NB+1)
*     The first element of U11 is in WORK(N+1,1)
      u11 = n
*     INVD is a block (N,2)
*     The first element of INVD is in WORK(1,INVD)
      invd = nb+2

      IF( upper ) THEN
*
*        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
*
        CALL ztrtri( uplo, 'U', n, a, lda, info )
*
*       inv(D) and inv(D)*inv(U)
* 
        k=1
        DO WHILE ( k .LE. n )
         IF( ipiv( k ).GT.0 ) THEN
*           1 x 1 diagonal NNB
             work(k,invd) = one / REAL ( A( K, K ) )
             work(k,invd+1) = 0
            k=k+1
         ELSE
*           2 x 2 diagonal NNB
             t = abs( work(k+1,1) )
             ak = REAL ( A( K, K ) ) / t
             akp1 = REAL ( A( K+1, K+1 ) ) / t
             akkp1 = work(k+1,1)  / t
             d = t*( ak*akp1-one )
             work(k,invd) = akp1 / d
             work(k+1,invd+1) = ak / d
             work(k,invd+1) = -akkp1 / d  
             work(k+1,invd) = dconjg(work(k,invd+1) )
            k=k+2
         END IF
        END DO
*
*       inv(U**H) = (inv(U))**H
*
*       inv(U**H)*inv(D)*inv(U)
*
        cut=n
        DO WHILE (cut .GT. 0)
           nnb=nb
           IF (cut .LE. nnb) THEN
              nnb=cut
           ELSE
              count = 0
*             count negative elements, 
              DO i=cut+1-nnb,cut
                  IF (ipiv(i) .LT. 0) count=count+1
              END DO
*             need a even number for a clear cut
              IF (mod(count,2) .EQ. 1) nnb=nnb+1
           END IF

           cut=cut-nnb
*
*          U01 Block 
*
           DO i=1,cut
             DO j=1,nnb
              work(i,j)=a(i,cut+j)
             END DO
           END DO
*
*          U11 Block
*
           DO i=1,nnb
             work(u11+i,i)=cone
             DO j=1,i-1
                work(u11+i,j)=zero
             END DO
             DO j=i+1,nnb
                work(u11+i,j)=a(cut+i,cut+j)
             END DO
           END DO
*
*          invD*U01
*
           i=1
           DO WHILE (i .LE. cut)
             IF (ipiv(i) > 0) THEN
                DO j=1,nnb
                    work(i,j)=work(i,invd)*work(i,j)
                END DO
                i=i+1
             ELSE
                DO j=1,nnb
                   u01_i_j = work(i,j)
                   u01_ip1_j = work(i+1,j)
                   work(i,j)=work(i,invd)*u01_i_j+
     $                      work(i,invd+1)*u01_ip1_j
                   work(i+1,j)=work(i+1,invd)*u01_i_j+
     $                      work(i+1,invd+1)*u01_ip1_j
                END DO
                i=i+2
             END IF
           END DO
*
*        invD1*U11
*
           i=1
           DO WHILE (i .LE. nnb)
             IF (ipiv(cut+i) > 0) THEN
                DO j=i,nnb
                    work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
                END DO
                i=i+1
             ELSE
                DO j=i,nnb
                   u11_i_j = work(u11+i,j)
                   u11_ip1_j = work(u11+i+1,j)
                work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
     $                      work(cut+i,invd+1)*work(u11+i+1,j)
                work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
     $                      work(cut+i+1,invd+1)*u11_ip1_j
                END DO
                i=i+2
             END IF
           END DO
*    
*       U11**H*invD1*U11->U11
*
        CALL ztrmm('L','U','C','U',nnb, nnb,
     $             cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
*
         DO i=1,nnb
            DO j=i,nnb
              a(cut+i,cut+j)=work(u11+i,j)
            END DO
         END DO
*
*          U01**H*invD*U01->A(CUT+I,CUT+J)
*
         CALL zgemm('C','N',nnb,nnb,cut,cone,a(1,cut+1),lda,
     $              work,n+nb+1, zero, work(u11+1,1), n+nb+1)
*
*        U11 =  U11**H*invD1*U11 + U01**H*invD*U01
*
         DO i=1,nnb
            DO j=i,nnb
              a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
            END DO
         END DO
*
*        U01 =  U00**H*invD0*U01
*
         CALL ztrmm('L',uplo,'C','U',cut, nnb,
     $             cone,a,lda,work,n+nb+1)

*
*        Update U01
*
         DO i=1,cut
           DO j=1,nnb
            a(i,cut+j)=work(i,j)
           END DO
         END DO
*
*      Next Block
*
       END DO
*
*        Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
*  
            i=1
            DO WHILE ( i .LE. n )
               IF( ipiv(i) .GT. 0 ) THEN
                  ip=ipiv(i)
                 IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
                 IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 i=i+1
                 IF ( (i-1) .LT. ip) 
     $                  CALL zheswapr( uplo, n, a, lda, i-1 ,ip )
                 IF ( (i-1) .GT. ip) 
     $                  CALL zheswapr( uplo, n, a, lda, ip ,i-1 )
              ENDIF
               i=i+1
            END DO
      ELSE
*
*        LOWER...
*
*        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
*
         CALL ztrtri( uplo, 'U', n, a, lda, info )
*
*       inv(D) and inv(D)*inv(U)
* 
        k=n
        DO WHILE ( k .GE. 1 )
         IF( ipiv( k ).GT.0 ) THEN
*           1 x 1 diagonal NNB
             work(k,invd) = one / REAL ( A( K, K ) )
             work(k,invd+1) = 0
            k=k-1
         ELSE
*           2 x 2 diagonal NNB
             t = abs( work(k-1,1) )
             ak = REAL ( A( K-1, K-1 ) ) / t
             akp1 = REAL ( A( K, K ) ) / t
             akkp1 = work(k-1,1) / t
             d = t*( ak*akp1-one )
             work(k-1,invd) = akp1 / d
             work(k,invd) = ak / d
             work(k,invd+1) = -akkp1 / d  
             work(k-1,invd+1) = dconjg(work(k,invd+1) )
            k=k-2
         END IF
        END DO
*
*       inv(U**H) = (inv(U))**H
*
*       inv(U**H)*inv(D)*inv(U)
*
        cut=0
        DO WHILE (cut .LT. n)
           nnb=nb
           IF (cut + nnb .GE. n) THEN
              nnb=n-cut
           ELSE
              count = 0
*             count negative elements, 
              DO i=cut+1,cut+nnb
                  IF (ipiv(i) .LT. 0) count=count+1
              END DO
*             need a even number for a clear cut
              IF (mod(count,2) .EQ. 1) nnb=nnb+1
           END IF
*      L21 Block
           DO i=1,n-cut-nnb
             DO j=1,nnb
              work(i,j)=a(cut+nnb+i,cut+j)
             END DO
           END DO
*     L11 Block
           DO i=1,nnb
             work(u11+i,i)=cone
             DO j=i+1,nnb
                work(u11+i,j)=zero
             END DO
             DO j=1,i-1
                work(u11+i,j)=a(cut+i,cut+j)
             END DO
           END DO
*
*          invD*L21
*
           i=n-cut-nnb
           DO WHILE (i .GE. 1)
             IF (ipiv(cut+nnb+i) > 0) THEN
                DO j=1,nnb
                    work(i,j)=work(cut+nnb+i,invd)*work(i,j)
                END DO
                i=i-1
             ELSE
                DO j=1,nnb
                   u01_i_j = work(i,j)
                   u01_ip1_j = work(i-1,j)
                   work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
     $                        work(cut+nnb+i,invd+1)*u01_ip1_j
                   work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
     $                        work(cut+nnb+i-1,invd)*u01_ip1_j
                END DO
                i=i-2
             END IF
           END DO
*
*        invD1*L11
*
           i=nnb
           DO WHILE (i .GE. 1)
             IF (ipiv(cut+i) > 0) THEN
                DO j=1,nnb
                    work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
                END DO
                i=i-1
             ELSE
                DO j=1,nnb
                   u11_i_j = work(u11+i,j)
                   u11_ip1_j = work(u11+i-1,j)
                work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
     $                      work(cut+i,invd+1)*u11_ip1_j
                work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
     $                      work(cut+i-1,invd)*u11_ip1_j
                END DO
                i=i-2
             END IF
           END DO
*    
*       L11**H*invD1*L11->L11
*
        CALL ztrmm('L',uplo,'C','U',nnb, nnb,
     $             cone,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
*
         DO i=1,nnb
            DO j=1,i
              a(cut+i,cut+j)=work(u11+i,j)
            END DO
         END DO
*
        IF ( (cut+nnb) .LT. n ) THEN
*
*          L21**H*invD2*L21->A(CUT+I,CUT+J)
*
         CALL zgemm('C','N',nnb,nnb,n-nnb-cut,cone,a(cut+nnb+1,cut+1)
     $             ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
       
*
*        L11 =  L11**H*invD1*L11 + U01**H*invD*U01
*
         DO i=1,nnb
            DO j=1,i
              a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
            END DO
         END DO
*
*        L01 =  L22**H*invD2*L21
*
         CALL ztrmm('L',uplo,'C','U', n-nnb-cut, nnb,
     $             cone,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)

*      Update L21
         DO i=1,n-cut-nnb
           DO j=1,nnb
              a(cut+nnb+i,cut+j)=work(i,j)
           END DO
         END DO
       ELSE
*
*        L11 =  L11**H*invD1*L11
*
         DO i=1,nnb
            DO j=1,i
              a(cut+i,cut+j)=work(u11+i,j)
            END DO
         END DO
       END IF
*
*      Next Block
*
           cut=cut+nnb
       END DO
*
*        Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
* 
            i=n
            DO WHILE ( i .GE. 1 )
               IF( ipiv(i) .GT. 0 ) THEN
                  ip=ipiv(i)
                 IF (i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip  )
                 IF (i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 IF ( i .LT. ip) CALL zheswapr( uplo, n, a, lda, i ,ip )
                 IF ( i .GT. ip) CALL zheswapr( uplo, n, a, lda, ip ,i )
                 i=i-1
               ENDIF
               i=i-1
            END DO
      END IF
*
      RETURN
*
*     End of ZHETRI2X
*

Here is the call graph for this function:

Here is the caller graph for this function: