Inverse of a UDU- or LDL-Factored Symmetric Matrix

Calling Sequence

CALL DSYTRI	(UPLO, N, DA, LDA, IPIVOT, DWORK, INFO)
CALL SSYTRI	(UPLO, N, SA, LDA, IPIVOT, SWORK, INFO)
CALL ZSYTRI	(UPLO, N, ZA, LDA, IPIVOT, ZWORK, INFO)
CALL CSYTRI	(UPLO, N, CA, LDA, IPIVOT, CWORK, INFO)
void dsytri	(char uplo, int n, double *da, int lda, int *ipivot, int *info)
void ssytri	(char uplo, int n, float *sa, int lda, int *ipivot, int *info)
void zsytri	(char uplo, int n, doublecomplex *za, int lda, int *ipivot, int *info)
void csytri	(char uplo, int n, complex *ca, int lda, int *ipivot, int *info)

Arguments

UPLO	Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.
	'U' or 'u'	xA contains the upper triangle.
	'L' or 'l'	xA contains the lower triangle.
N	Order of the matrix A. N 0.
xA	On entry, the UDU or LDL factorization of the matrix A, as computed by xSYTRF. On exit, the symmetric inverse of the matrix A. If UPLO = 'U' or 'u', the upper triangular part of the inverse is formed and the part of A below the diagonal is not used. If UPLO = 'L' or 'l' the lower triangular part of the inverse is formed and the part of A above the diagonal is not used.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, N).
IPIVOT	Pivot indices as computed by xSYTRF.
xWORK	Scratch array with a dimension of N for real subroutines or 2 × N for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	INFO > 0	D(i,i), where i = INFO, is zero. The matrix is therefore singular and its inverse could not be computed.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDWORK, N

PARAMETER (N = 4)

PARAMETER (LDA = N)

PARAMETER (LDWORK = N)

C

DOUBLE PRECISION A(LDA,N), WORK(LDWORK)

INTEGER ICOL, INFO, IPIVOT(N), IROW

C

EXTERNAL DSYTRF, DSYTRI

INTRINSIC ABS

C

C Initialize the array A to store in symmetric form the 4x4

C symmetric coefficient matrix A shown below.

C

C 2 1 0 0

C A = 1 2 1 0

C 0 1 2 1

C 0 0 1 1

C

DATA A / 2.0D0, 3*8D8, 1.0D0, 2.0D0, 2*8D8, 0.0D0, 1.0D0,

$ 2.0D0, 8D8, 0.0D0, 0.0D0, 1.0D0, 1.0D0 /

C

C Print the initial values of the arrays.

C

PRINT 1000

DO 100, IROW = 1, N

PRINT 1010, (A(ICOL,IROW), ICOL = 1, IROW - 1),

$ (A(IROW,ICOL), ICOL = IROW, N)

100 CONTINUE

PRINT 1020

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)

C

C UDU factor A.

C

CALL DSYTRF ('UPPER TRIANGLE OF A STORED', N, A, LDA,

$ IPIVOT, WORK, LDWORK, INFO)

IF (INFO .NE. 0) THEN

PRINT 1030, INFO

STOP 1

END IF

C

C Use the factored form of A to compute the inverse of A and

C print the inverse thus computed.

C

CALL DSYTRI ('UPPER TRIANGULAR FACTOR', N, A, LDA, IPIVOT,

$ WORK, INFO)

IF (INFO .NE. 0) THEN

PRINT 1040, ABS(INFO)

STOP 2

END IF

PRINT 1050

DO 200, IROW = 1, N

PRINT 1010, (A(ICOL,IROW), ICOL = 1, IROW - 1),

$ (A(IROW,ICOL), ICOL = IROW, N)

200 CONTINUE

C

1000 FORMAT (1X, 'A in full form:')

1010 FORMAT (4(3X, F6.3))

1020 FORMAT (/1X, 'A in symmetric form: (* in unused elements)')

1030 FORMAT (1X, 'Error factoring A, INFO = ', I5)

1040 FORMAT (1X, 'Error computing inverse of A, INFO = ', I5)

1050 FORMAT (/1X, 'A**(-1):')

C

END

Sample Output

A in full form:
2.000 1.000 0.000 0.000
1.000 2.000 1.000 0.000
0.000 1.000 2.000 1.000
0.000 0.000 1.000 1.000
A in symmetric form: (* in unused elements)
2.000 1.000 0.000 0.000
****** 2.000 1.000 0.000
****** ****** 2.000 1.000
****** ****** ****** 1.000
A**(-1):
1.000 -1.000 1.000 -1.000
-1.000 2.000 -2.000 2.000
1.000 -2.000 3.000 -3.000
-1.000 2.000 -3.000 4.000