
PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDAB, LDAG, LDB, M, N, NRHS, NSUB, NSUPER

PARAMETER (NSUB = 1)

PARAMETER (NSUPER = 1)

PARAMETER (M = 4)

PARAMETER (N = 4)

PARAMETER (NRHS = 1)

PARAMETER (LDAB = 2*NSUB + NSUPER + 1)

PARAMETER (LDAG = M)

PARAMETER (LDB = N)

C

DOUBLE PRECISION AB(LDAB,N), AG(LDAG,N), B(LDB,NRHS)

INTEGER ICOL, INFO, IPIVOT(N), IROW, IROWB, I1, I2

INTEGER M1

C

EXTERNAL DGBTRF, DGBTRS

INTRINSIC MAX0, MIN0

C

C Initialize the array AG to store the 4x4 matrix A with one

C subdiagonal and one superdiagonal shown below. Initialize

C the array B to store the vector b shown below.

C

C 2 1 5

C AG = 1 2 1 b = 5

C 1 2 1 5

C 1 2 5

C

DATA AB / 16 * 8.0D8 /

DATA AG / 2.0D0, 1.0D0, 2*0.0, 1.0D0, 2.0D0, 1.0D0,

$ 0.0D0, 0.0D0, 1.0D0, 2.0D0, 1.0D0, 2*0D0,

$ 1.0D0, 2.0D0 /

DATA B / N*5.0D0 /

C

C Copy the matrix A from the array AG to the array AB. The

C matrix is stored in general storage mode in AG and it will

C be stored in banded storage mode in AB. The code to copy

C from general to banded storage mode was written for LINPACK

C by Cleve Moler.

C

M1 = NSUB + NSUPER + 1

DO 20, ICOL = 1, N

I1 = MAX0 (1, ICOL  NSUPER)

I2 = MIN0 (N, ICOL + NSUB)

DO 10, IROW = I1, I2

IROWB = IROW  ICOL + M1

AB(IROWB,ICOL) = AG(IROW,ICOL)

10 CONTINUE

20 CONTINUE

C

C Print the initial values of the arrays.

C

PRINT 1000

PRINT 1010, ((AG(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

PRINT 1020

PRINT 1010, ((AB(IROW,ICOL), ICOL = 1, N), IROW = 1, LDAB)

PRINT 1030

PRINT 1040, B

C

C Factor the matrix in banded form and print the results.

C

CALL DGBTRF (M, N, NSUB, NSUPER, AB, LDAB, IPIVOT, INFO)

IF (INFO .EQ. 0) THEN

CALL DGBTRS ('NO TRANSPOSE A', N, NSUB, NSUPER, NRHS,

$ AB, LDAB, IPIVOT, B, LDB, INFO)

IF (INFO .EQ. 0) THEN

PRINT 1050

PRINT 1040, B

ELSE

PRINT 1060, INFO

END IF

ELSE

PRINT 1070, INFO

END IF

C

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

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

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

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

1040 FORMAT (1X, 2X, F4.1)

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

1060 FORMAT (1X, 'Error in solving system, INFO = ', I4)

1070 FORMAT (1X, 'Error in factoring system, INFO = ', I4)

C

END

