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

CSYTRI2X

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

Purpose:

 CSYTRI2X computes the inverse of a real symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 CSYTRF.

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**T; = 'L': Lower triangular, form is A = L*D*L**T.
[in]	N	N is INTEGER The order of the matrix A. N >= 0.
[in,out]	A	A is COMPLEX 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 CSYTRF. 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 CSYTRF.
[out]	WORK	WORK is COMPLEX array, dimension (N+NNB+1,NNB+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 2011

Definition at line 122 of file csytri2x.f.

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

      COMPLEX   ak, akkp1, akp1, d, t
      COMPLEX   u01_i_j, u01_ip1_j
      COMPLEX   u11_i_j, u11_ip1_j
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           csyconv, xerbla, ctrtri
      EXTERNAL           cgemm, ctrmm, csyswapr
*     ..
*     .. 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( 'CSYTRI2X', -info )
         RETURN
      END IF
      IF( n.EQ.0 )
     $   RETURN
*
*     Convert A
*     Workspace got Non-diag elements of D
*
      CALL csyconv( 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**T)*inv(D)*inv(U)*P**T.
*
        CALL ctrtri( 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 /  a( k, k )
             work(k,invd+1) = 0
            k=k+1
         ELSE
*           2 x 2 diagonal NNB
             t = work(k+1,1)
             ak = a( k, k ) / t
             akp1 = 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) = -akkp1 / d  
            k=k+2
         END IF
        END DO
*
*       inv(U**T) = (inv(U))**T
*
*       inv(U**T)*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)=one
             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**T*invD1*U11->U11
*
        CALL ctrmm('L','U','T','U',nnb, nnb,
     $             one,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**T*invD*U01->A(CUT+I,CUT+J)
*
         CALL cgemm('T','N',nnb,nnb,cut,one,a(1,cut+1),lda,
     $              work,n+nb+1, zero, work(u11+1,1), n+nb+1)
*
*        U11 =  U11**T*invD1*U11 + U01**T*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**T*invD0*U01
*
         CALL ctrmm('L',uplo,'T','U',cut, nnb,
     $             one,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**T: P * inv(U**T)*inv(D)*inv(U) *P**T
*  
            i=1
            DO WHILE ( i .LE. n )
               IF( ipiv(i) .GT. 0 ) THEN
                  ip=ipiv(i)
                 IF (i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip )
                 IF (i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 i=i+1
                 IF ( (i-1) .LT. ip) 
     $                  CALL csyswapr( uplo, n, a, lda, i-1 ,ip )
                 IF ( (i-1) .GT. ip) 
     $                  CALL csyswapr( uplo, n, a, lda, ip ,i-1 )
              ENDIF
               i=i+1
            END DO
      ELSE
*
*        LOWER...
*
*        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
*
         CALL ctrtri( 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 /  a( k, k )
             work(k,invd+1) = 0
            k=k-1
         ELSE
*           2 x 2 diagonal NNB
             t = work(k-1,1)
             ak = a( k-1, k-1 ) / t
             akp1 = 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) = -akkp1 / d  
            k=k-2
         END IF
        END DO
*
*       inv(U**T) = (inv(U))**T
*
*       inv(U**T)*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)=one
             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**T*invD1*L11->L11
*
        CALL ctrmm('L',uplo,'T','U',nnb, nnb,
     $             one,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**T*invD2*L21->A(CUT+I,CUT+J)
*
         CALL cgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
     $             ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
       
*
*        L11 =  L11**T*invD1*L11 + U01**T*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**T*invD2*L21
*
         CALL ctrmm('L',uplo,'T','U', n-nnb-cut, nnb,
     $             one,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**T*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**T: P * inv(U**T)*inv(D)*inv(U) *P**T
* 
            i=n
            DO WHILE ( i .GE. 1 )
               IF( ipiv(i) .GT. 0 ) THEN
                  ip=ipiv(i)
                 IF (i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip  )
                 IF (i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 IF ( i .LT. ip) CALL csyswapr( uplo, n, a, lda, i ,ip )
                 IF ( i .GT. ip) CALL csyswapr( uplo, n, a, lda, ip ,i )
                 i=i-1
               ENDIF
               i=i-1
            END DO
      END IF
*
      RETURN
*
*     End of CSYTRI2X
*

Here is the call graph for this function:

Here is the caller graph for this function: