Cholesky Factorization of a Symmetric Positive Definite Matrix in Packed Storage

Calling Sequence

CALL DPPFA	(DA, N, INFO)
CALL SPPFA	(SA, N, INFO)
CALL ZPPFA	(ZA, N, INFO)
CALL CPPFA	(CA, N, INFO)
void dppfa	(double *da, int n, int *info)
void sppfa	(float *sa, int n, int *info)
void zppfa	(doublecomplex *za, int n, int *info)
void cppfa	(complex *ca, int n, int *info)

Arguments

xA	On entry, the upper triangle of the matrix A. On exit, a Cholesky factorization of the matrix A.
N	Order of the matrix A. N 0.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k if U(k,k) = 0 to indicate that xGESL will divide by zero if called.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LENGTA, N

PARAMETER (N = 4)

PARAMETER (LENGTA = (N * N + N) / 2)

C

DOUBLE PRECISION A(LENGTA), B(N)

INTEGER INFO

C

EXTERNAL DPPFA, DPPSL

C

C Initialize the array A to store in packed symmetric storage

C mode the matrix A shown below. Initialize the array B to

C store the vector B shown below.

C

C 4 3 2 1 30

C A = 3 4 3 2 b = 20

C 2 3 4 3 20

C 1 2 3 4 40

C

DATA A / 4.0D0, 3.0D0, 4.0D0, 2.0D0, 3.0D0, 4.0D0,

$ 1.0D0, 2.0D0, 3.0D0, 4.0D0 /

DATA B / 6.0D1, 2.0D1, 2.0D1, 6.0D1 /

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)

PRINT 1020

PRINT 1030, B

CALL DPPFA (A, N, INFO)

IF (INFO .EQ. 0) THEN

CALL DPPSL (A, N, B)

PRINT 1040

PRINT 1030, B

ELSE

PRINT 1050

END IF

C

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

1010 FORMAT (5(3X, F5.1))

1020 FORMAT (/1X, 'b:')

1030 FORMAT (3X, F5.1)

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

1050 FORMAT (/1X, 'A is too ill-conditioned.')

C

END

Sample Output

A in full form:
4.0 3.0 2.0 1.0
3.0 4.0 3.0 2.0
2.0 3.0 4.0 3.0
1.0 2.0 3.0 4.0
b:
60.0
20.0
20.0
60.0
A**(-1) * b:
32.0
-20.0
-20.0
32.0