Inverse of a Cholesky-Factored Symmetric Positive Definite Matrix in Packed Storage

Calling Sequence

CALL DPPTRI	(UPLO, N, DA, INFO)
CALL SPPTRI	(UPLO, N, SA, INFO)
CALL ZPPTRI	(UPLO, N, ZA, INFO)
CALL CPPTRI	(UPLO, N, CA, INFO)
void dpptri	(char uplo, int n, double *da, int *info)
void spptri	(char uplo, int n, float *sa, int *info)
void zpptri	(char uplo, int n, doublecomplex *za, int *info)
void cpptri	(char uplo, int n, complex *ca, 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 Cholesky factorization of the matrix A as computed by xPPTRF. The dimension of xA is (N × N + N) / 2. On exit, the symmetric (or Hermitian) inverse of A.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	INFO > 0	Element (i,i) of the factor U or L, 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)

INTEGER INFO

C

EXTERNAL DPPTRF, DPPTRI

INTRINSIC ABS

C

C Initialize the array A to store in symmetric form the

C 4x4 symmetric positive definite coefficient matrix A

C shown below. Initialize the array B to store the right

C hand side vector b 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 Cholesky factor A.

C

CALL DPPTRF ('UPPER TRIANGLE OF A STORED', N, A, 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 DPPTRI ('UPPER TRIANGLE OF A STORED', N, A, 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