Inverse of a UDU- or LDL-Factored Symmetric Matrix in Packed Storage

Calling Sequence

CALL DSPTRI	(UPLO, N, DA, IPIVOT, DWORK, INFO)
CALL SSPTRI	(UPLO, N, SA, IPIVOT, SWORK, INFO)
CALL ZSPTRI	(UPLO, N, ZA, IPIVOT, ZWORK, INFO)
CALL CSPTRI	(UPLO, N, CA, IPIVOT, CWORK, INFO)
void dsptri	(char uplo, int n, double *da, int *ipivot, int *info)
void ssptri	(char uplo, int n, float *sa, int *ipivot, int *info)
void zsptri	(char uplo, int n, doublecomplex *za, int *ipivot, int *info)
void csptri	(char uplo, int n, complex *ca, 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 xSPTRF. The dimension of xA is (N × N + N) / 2. On exit, the symmetric inverse of the matrix A.
IPIVOT	Pivot indices as computed by xSPTRF.
xWORK	Scratch array with a dimension of N.
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, N

PARAMETER (N = 4)

PARAMETER (LDA = (N * (N + 1)) / 2)

C

DOUBLE PRECISION A(LDA), WORK(N)

INTEGER INFO, IPIVOT(N)

C

EXTERNAL DSPTRF, DSPTRI

INTRINSIC ABS

C

C Initialize the array A to store in symmetric packed form the

C 4x4 symmetric 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, 1.0D0, 2.0D0, 0.0D0, 1.0D0, 2.0D0,

$ 0.0D0, 0.0D0, 1.0D0, 1.0D0 /

C

C Print the initial values of the arrays.

C

PRINT 1000

PRINT 1010, A(1), A(2), A(4), A(7)

PRINT 1010, A(2), A(3), A(5), A(8)

PRINT 1010, A(4), A(5), A(6), A(9)

PRINT 1010, A(7), A(8), A(9), A(10)

C

C LDL factor A.

C

CALL DSPTRF ('UPPER TRIANGLE OF A STORED', N, A, IPIVOT,

$ INFO)

IF (INFO .NE. 0) THEN

PRINT 1020, 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 DSPTRI ('UPPER TRIANGLE OF A STORED', N, A, IPIVOT,

$ WORK, INFO)

IF (INFO .NE. 0) THEN

PRINT 1030, ABS(INFO)

STOP 2

END IF

PRINT 1040

PRINT 1010, A(1), A(2), A(4), A(7)

PRINT 1010, A(2), A(3), A(5), A(8)

PRINT 1010, A(4), A(5), A(6), A(9)

PRINT 1010, A(7), A(8), A(9), A(10)

C

1000 FORMAT (1X, 'A:')

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

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

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

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

C

END

Sample Output

A:
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**(-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