
PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDAF, LDB, LDIWRK, LDWORK, LDX, N

INTEGER NDIAG, NRHS

PARAMETER (N = 4)

PARAMETER (NDIAG = 1)

PARAMETER (NRHS = 1)

PARAMETER (LDA = NDIAG + 1)

PARAMETER (LDAF = LDA)

PARAMETER (LDB = N)

PARAMETER (LDIWRK = N)

PARAMETER (LDX = N)

PARAMETER (LDWORK = 3 * N)

C

DOUBLE PRECISION A(LDA,N), AF(LDAF,N), B(LDB,NRHS)

DOUBLE PRECISION BERR(NRHS), FERR(NRHS), RCOND, SCALE(N)

DOUBLE PRECISION WORK(LDWORK), X(LDX,NRHS)

INTEGER ICOL, INFO, IROW, IWORK(LDIWRK)

CHARACTER*1 EQUED

C

EXTERNAL DPBSVX

C

C Initialize the array A to store in symmetric banded form

C the 4x4 symmetric positive definite coefficient matrix A

C with one subdiagonal and one superdiagonal shown below.

C Initialize the array B to store the right hand side

C vector b shown below.

C

C 2 1 0 0 6

C A = 1 2 1 0 b = 12

C 0 1 2 1 12

C 0 0 1 2 6

C

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

$ 1.0D0, 2.0D0, 1.0D0, 2.0D0 /

DATA B / 6.0D0, 1.2D1, 1.2D1, 6.0D0 /

C

C Print the initial values of the arrays.

C

PRINT 1000

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

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

PRINT 1010, 0.0D0, A(1,3), A(2,3), A(1,4)

PRINT 1010, 0.0D0, 0.0D0, A(1,4), A(2,4)

PRINT 1020

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

PRINT 1030

PRINT 1040, B

C

C Solve the system and print the results.

C

CALL DPBSVX ('EQUILIBRATE A', 'UPPER TRIANGLE OF A STORED',

$ N, NDIAG, NRHS, A, LDA, AF, LDAF, EQUED, SCALE,

$ B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK,

$ INFO)

IF (INFO .NE. 0) THEN

PRINT 1050, INFO

IF (INFO .EQ. (N + 1)) THEN

PRINT 1060

STOP 1

END IF

STOP 2

END IF

PRINT 1070

PRINT 1040, X

PRINT 1080, 1.0D0 / RCOND

IF (EQUED .EQ. 'N') THEN

PRINT 1090

ELSE

PRINT 1100

END IF

PRINT 1110, (IROW, BERR(IROW), IROW = 1, NRHS)

PRINT 1120, (IROW, FERR(IROW), IROW = 1, NRHS)

C

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

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

1020 FORMAT (/1X, 'A in banded form: (* in unused elements)')

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

1040 FORMAT (1X, F6.2)

1050 FORMAT (1X, 'Error solving Ax=b, INFO = ', I5)

1060 FORMAT (1X, 'Matrix is singular to working precision.')

1070 FORMAT (/1X, 'x:')

1080 FORMAT (/1X, 'Estimated condition number of A: ', F7.2)

1090 FORMAT (1X, 'No equilibration was required for A.')

1100 FORMAT (1X, 'Diagonal equilibration was done on A.')

1110 FORMAT (/1X, 'Backward error bound for system #', I1, ': ',

$ E12.5)

1120 FORMAT (1X, 'Forward error bound for system #', I1, ': ',

$ E12.5)

C

END

