Solution to a Linear System in a Cholesky-Factored Symmetric Positive Definite Matrix in Banded Storage

Calling Sequence

CALL DPBSL	(DA, LDA, N, NDIAG, DB)
CALL SPBSL	(SA, LDA, N, NDIAG, SB)
CALL ZPBSL	(ZA, LDA, N, NDIAG, ZB)
CALL CPBSL	(CA, LDA, N, NDIAG, CB)
void dpbsl	(double *da, int lda, int n, int ndiag, double *db)
void spbsl	(float *sa, int lda, int n, int ndiag, float *sb)
void zpbsl	(doublecomplex *za, int lda, int n, int ndiag, doublecomplex *zb)
void cpbsl	(complex *ca, int lda, int n, int ndiag, complex *cb)

Arguments

xA	Cholesky factorization of the matrix A, as computed by xPBCO or xPBFA.
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 of matrix A. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.
xB	On entry, the right-hand side vector b. On exit, the solution vector x.

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)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DPBFA, 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 = 120

C 0 -1 2 -1 120

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 / 6.0D1, 1.2D1, 1.2D1, 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 DPBFA (A, LDA, N, NDIAG, INFO)

IF (INFO .EQ. 0) THEN

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

PRINT 1070

PRINT 1060, B

ELSE

PRINT 1080

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, 'A**(-1) * b:')

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

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
12.0
12.0
60.0
A**(-1) * b:
72.0
84.0
84.0
72.0