Cholesky Factorization and Condition Number of a Symmetric Positive Definite Matrix in Banded Storage

Calling Sequence

CALL DPBCO	(DA, LDA, N, NDIAG, DRCOND, DWORK, INFO)
CALL SPBCO	(SA, LDA, N, NDIAG, SRCOND, SWORK, INFO)
CALL ZPBCO	(ZA, LDA, N, NDIAG, DRCOND, ZWORK, INFO)
CALL CPBCO	(CA, LDA, N, NDIAG, SRCOND, CWORK, INFO)
void dpbco	(double *da, int lda, int n, int ndiag, double *drcond, int *info)
void spbco	(float *sa, int lda, int n, int ndiag, float *srcond, int *info)
void zpbco	(doublecomplex *za, int lda, int n, int ndiag, double *drcond, int *info)
void cpbco	(complex *ca, int lda, int n, int ndiag, float *srcond, int *info)

Arguments

xA	On entry, the upper triangle of the matrix A. On exit, a Cholesky factorization of the matrix A.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA NDIAG + 1.
N	Order of the matrix A. N 0.
NDIAG	Number of diagonals. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.
xRCOND	On exit, an estimate of the reciprocal condition number of A. 0.0 RCOND 1.0. As the value of RCOND gets smaller, operations with A such as solving Ax = b may become less stable. If RCOND satisfies RCOND + 1.0 = 1.0 then A may be singular to working precision.
xWORK	Scratch array with a dimension of N.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k if the leading minor of order k is not positive definite.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, N, NDIAG

PARAMETER (N = 4)

PARAMETER (NDIAG = 1)

PARAMETER (LDA = NDIAG + 1)

C

DOUBLE PRECISION A(LDA,N), B(N), RCOND, WORK(N)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DPBCO, DPBSL

C

C Initialize the array A to store in banded storage mode

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

C store the vector B shown below.

C

C 2 -1 0 0 60

C A = -1 2 -1 0 b = 60

C 0 -1 2 -1 60

C 0 0 -1 2 60

C

DATA A / 8D8, 2.0D0, -1.0D0, 2.0D0, -1.0D0, 2.0D0, -1.0D0,

$ 2.0D0 /

DATA B / 4*6.0D1 /

C

PRINT 1000

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

PRINT 1010, A(3,1), A(2,2), A(1,3)

PRINT 1020, A(3,2), A(2,3), A(1,4)

PRINT 1030, A(3,3), A(2,4)

PRINT 1040

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

PRINT 1050

PRINT 1060, B

CALL DPBCO (A, LDA, N, NDIAG, RCOND, WORK, INFO)

IF (INFO .EQ. 0) THEN

IF ((RCOND + 1.0D0) .EQ. 1.0D0) THEN

PRINT 1100

END IF

CALL DPBSL (A, LDA, N, NDIAG, B)

PRINT 1070, RCOND

PRINT 1080

PRINT 1060, B

ELSE

PRINT 1090

END IF

C

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

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

1020 FORMAT (8X, 3(3X, F5.1))

1030 FORMAT (16X, 3(3X, F5.1))

1040 FORMAT (/1X, 'A in banded form: (* in unused entries)')

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

1060 FORMAT (3X, F5.1)

1070 FORMAT (/1X, 'Reciprocal condition number of A:', F5.1)

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

1090 FORMAT (/1X, 'A is not positive definite.')

1100 FORMAT (1X, 'A may be singular to working precision.')

C

END

Sample Output

A in full form:
2.0 -1.0
-1.0 2.0 -1.0
-1.0 2.0 -1.0
-1.0 2.0
A in banded form: (* in unused entries)
***** -1.0 -1.0 -1.0
2.0 2.0 2.0 2.0
b:
60.0
60.0
60.0
60.0
Reciprocal condition number of A: 0.1
A**(-1) * b:
120.0
180.0
180.0
120.0