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

SSYTRI2X

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

Purpose:

 SSYTRI2X 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
 SSYTRF.

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 REAL 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 SSYTRF. 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 SSYTRF.
[out]	WORK	WORK is REAL 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 ssytri2x.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( * )
      REAL               a( lda, * ), work( n+nb+1,* )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               one, zero
      parameter                ( one = 1.0e+0, zero = 0.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            upper
      INTEGER            i, iinfo, ip, k, cut, nnb
      INTEGER            count
      INTEGER            j, u11, invd

      REAL               ak, akkp1, akp1, d, t
      REAL               u01_i_j, u01_ip1_j
      REAL               u11_i_j, u11_ip1_j
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           ssyconv, xerbla, strtri
      EXTERNAL           sgemm, strmm, ssyswapr
*     ..
*     .. 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( 'SSYTRI2X', -info )
         RETURN
      END IF
      IF( n.EQ.0 )
     $   RETURN
*
*     Convert A
*     Workspace got Non-diag elements of D
*
      CALL ssyconv( 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 strtri( 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 strmm('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 sgemm('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 strmm('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 ssyswapr( uplo, n, a, lda, i ,ip )
                 IF (i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 i=i+1
                 IF ( (i-1) .LT. ip) 
     $                  CALL ssyswapr( uplo, n, a, lda, i-1 ,ip )
                 IF ( (i-1) .GT. ip) 
     $                  CALL ssyswapr( 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 strtri( 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 .GT. 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 strmm('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 sgemm('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 strmm('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 ssyswapr( uplo, n, a, lda, i ,ip  )
                 IF (i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
               ELSE
                 ip=-ipiv(i)
                 IF ( i .LT. ip) CALL ssyswapr( uplo, n, a, lda, i ,ip )
                 IF ( i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
                 i=i-1
               ENDIF
               i=i-1
            END DO
      END IF
*
      RETURN
*
*     End of SSYTRI2X
*

Here is the call graph for this function:

Here is the caller graph for this function: