|
PROGRAM TEST
|
IMPLICIT NONE
|
C
|
INTEGER LDAB, LDAFB, LDB, LDX, N, NRHS, NSUB
|
INTEGER NSUPER, SIZEAB, SIZEB, SIZEX
|
PARAMETER (NSUB = 1)
|
PARAMETER (NSUPER = 1)
|
PARAMETER (N = 4)
|
PARAMETER (NRHS = 1)
|
PARAMETER (LDAB = NSUB + NSUPER + 1)
|
PARAMETER (LDAFB = 2*NSUB + NSUPER + 1)
|
PARAMETER (LDB = N)
|
PARAMETER (LDX = N)
|
PARAMETER (SIZEAB = LDAB * N)
|
PARAMETER (SIZEB = LDB * NRHS)
|
PARAMETER (SIZEX = LDX * NRHS)
|
C
|
DOUBLE PRECISION AB(LDAB,N), AFB(LDAFB,N), B(LDB,NRHS)
|
DOUBLE PRECISION BERR(NRHS), COLSCA(N), FERR(NRHS), RCOND
|
DOUBLE PRECISION ROWSCA(N), WORK(3*N), X(LDX,NRHS)
|
INTEGER ICOL, INFO, IPIVOT(N), IROW, IWORK(N)
|
CHARACTER*1 EQUED
|
C
|
EXTERNAL DGBSVX
|
INTRINSIC MAX0, MIN0
|
C
|
C Initialize the array AB to store in banded storage mode the
|
C 4x4 matrix A with one subdiagonal and one superdiagonal
|
C shown below. Initialize the array B to store the vector b
|
C shown below.
|
C
|
C 2 -1 5
|
C A = -1 2 -1 b = 5
|
C -1 2 -1 5
|
C -1 2 5
|
C
|
DATA AB / 8.8D8, 2.0D0, -1.0D0, -1.0D0, 2.0D0, -1.0D0,
|
$ -1.0D0, 2.0D0, -1.0D0, -1.0D0, 2.0D0, 8.8D8 /
|
DATA B / SIZEB*5.0D0 /
|
DATA X / SIZEX*9.9D9 /
|
C
|
C Print the initial values of the arrays.
|
C
|
PRINT 1000
|
PRINT 1010, AB(2,1), AB(1,2), 0.0D0, 0.0D0
|
PRINT 1010, AB(3,1), AB(2,2), AB(1,3), 0.0D0
|
PRINT 1010, 0.0D0, AB(3,2), AB(2,3), AB(1,4)
|
PRINT 1010, 0.0D0, 0.0D0, AB(3,3), AB(2,4)
|
PRINT 1020
|
PRINT 1010, ((AB(IROW,ICOL), ICOL = 1, N), IROW = 1, LDAB)
|
PRINT 1030
|
PRINT 1040, B
|
C
|
C Solve the system Ax=b.
|
C
|
CALL DGBSVX ('EQUILIBRATE A', 'NO TRANSPOSE A', N, NSUB,
|
$ NSUPER, NRHS, AB, LDAB, AFB, LDAFB, IPIVOT,
|
$ EQUED, ROWSCA, COLSCA, B, LDB, X, LDX, RCOND,
|
$ FERR, BERR, WORK, IWORK, INFO)
|
IF (INFO .EQ. 0) THEN
|
PRINT 1050
|
PRINT 1040, X
|
PRINT 1060, 1.0D0 / RCOND
|
IF (EQUED .EQ. 'N') PRINT 1070
|
IF (EQUED .EQ. 'R') PRINT 1080
|
IF (EQUED .EQ. 'C') PRINT 1090
|
IF (EQUED .EQ. 'B') PRINT 1100
|
PRINT 1110, (IROW, BERR(IROW), IROW = 1, NRHS)
|
PRINT 1120, (IROW, FERR(IROW), IROW = 1, NRHS)
|
ELSE
|
PRINT 1130, INFO
|
IF (INFO .EQ. (N + 1)) THEN
|
PRINT 1140
|
END IF
|
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, 'Estimated condition number of A: ', F7.2)
|
1070 FORMAT (1X, 'No equilibration was required for A.')
|
1080 FORMAT (1X, 'Row equilibration was done on A.')
|
1090 FORMAT (1X, 'Column equilibration was done on A.')
|
1100 FORMAT (1X, 'Row and column equilibration was done on A.')
|
1110 FORMAT (/1X, 'Backward error for system #', I1, ': ', E12.6)
|
1120 FORMAT (1X, 'Forward error for system #', I1, ': ', E12.6)
|
1130 FORMAT (1X, 'Error in solving system, INFO = ', I3)
|
1140 FORMAT (1X, 'Matrix is singular to working precision.')
|
C
|
END
|
|