|
PROGRAM TEST
|
IMPLICIT NONE
|
C
|
DOUBLE PRECISION ZERO
|
INTEGER LDB, LDB1, LDWORK, LDWRK2, LDX, LENGTA
|
INTEGER N, NRHS
|
PARAMETER (N = 3)
|
PARAMETER (LDB = N)
|
PARAMETER (LDB1 = N)
|
PARAMETER (LDWORK = 2 * N)
|
PARAMETER (LDWRK2 = N)
|
PARAMETER (LDX = LDB)
|
PARAMETER (LENGTA = (N * N + N) / 2)
|
PARAMETER (NRHS = 1)
|
PARAMETER (ZERO = 0.0D0)
|
C
|
DOUBLE PRECISION ANORM, BERR(N), FERR(N), RCOND
|
DOUBLE PRECISION WORK2(LDWRK2)
|
COMPLEX*16 A(LENGTA), AF(LENGTA), B(LDB,NRHS)
|
COMPLEX*16 B1(LDB,NRHS), WORK(LDWORK), X(LDX,NRHS)
|
INTEGER ICOL, INFO, IPIVOT(N), IROW
|
C
|
EXTERNAL ZCOPY, ZHPCON, ZHPMV, ZHPRFS, ZHPTRF
|
EXTERNAL ZHPTRS
|
INTRINSIC ABS, CMPLX, CONJG, DBLE, MAX
|
C
|
C Initialize the array A to store the coefficient array A
|
C shown below. Initialize the array B to store the right
|
C hand side vector b shown below.
|
C
|
C (1e+16, 0) (1e+16, 1) (3e+16, 3e-16)
|
C A = (1e+16, -1) (1e+16, 0) (3e+26, 3e-16)
|
C (3e+16, 3e-16) (3e+26, -3e-16) (7e-16, 0)
|
C
|
C (2e-26, 2e+16)
|
C b = (2e+26, 2e-16)
|
C (2, 2)
|
C
|
DATA A / (1.0D16,8D8), (1.0D16,1.0D0), (1.0D16,8D8),
|
$ (3.0D16,3.0D-16), (3.0D26,3.0D-16), (7.0D-16,8D8) /
|
DATA B / (2.0D-26,2.0D16), (2.0D26,2.0D-16), (2.0D0,2.0D0) /
|
C
|
C Print the initial value of the arrays.
|
C
|
PRINT 1000
|
PRINT 1010, CMPLX(DBLE(A(1)),ZERO), A(2), A(4)
|
PRINT 1010, CONJG(A(2)), CMPLX(DBLE(A(3)),ZERO), A(5)
|
PRINT 1010, CONJG(A(4)), CONJG(A(5)), CMPLX(DBLE(A(6)),ZERO)
|
DO 110, ICOL = 1, NRHS
|
PRINT 1020, ICOL
|
DO 100, IROW = 1, N
|
PRINT 1100, B(IROW,ICOL)
|
100 CONTINUE
|
110 CONTINUE
|
C
|
C UDU factor A.
|
C
|
CALL ZCOPY (LENGTA, A, 1, AF, 1)
|
CALL ZHPTRF ('UPPER TRIANGLE OF A STORED', N, AF, IPIVOT,
|
$ INFO)
|
IF (INFO .LT. 0) THEN
|
PRINT 1030, ABS(INFO)
|
STOP 1
|
ELSE IF (INFO .GT. 0) THEN
|
PRINT 1040
|
STOP 2
|
END IF
|
C
|
C Estimate the condition number. Print a warning if it is
|
C unreasonably high.
|
C
|
ANORM = 3.0D26
|
CALL ZHPCON ('UPPER TRIANGULAR FACTOR OF A', N, AF, IPIVOT,
|
$ ANORM, RCOND, WORK, INFO)
|
IF (INFO .LT. 0) THEN
|
PRINT 1050, ABS(INFO)
|
STOP 3
|
END IF
|
PRINT 1060, RCOND
|
IF ((1.0D0 + RCOND) .EQ. 1.0D0) THEN
|
PRINT 1070
|
END IF
|
C
|
C Compute and print the solution.
|
C
|
CALL ZCOPY (LDB * NRHS, B, 1, X, 1)
|
CALL ZHPTRS ('UPPER TRIANGULAR FACTOR OF A', N, NRHS, AF,
|
$ IPIVOT, X, LDX, INFO)
|
IF (INFO .LT. 0) THEN
|
PRINT 1080, ABS(INFO)
|
STOP 4
|
END IF
|
DO 210, ICOL = 1, NRHS
|
PRINT 1090, ICOL
|
DO 200, IROW = 1, N
|
PRINT 1100, X(IROW,ICOL)
|
200 CONTINUE
|
210 CONTINUE
|
DO 220, ICOL = 1, NRHS
|
CALL ZHPMV ('UPPER TRIANGLE OF A STORED', N,
|
$ (1.0D0,0.0D0), A, X(1,ICOL), 1, (0.0D0,0.0D0),
|
$ B1(1,ICOL), 1)
|
220 CONTINUE
|
DO 240, ICOL = 1, NRHS
|
PRINT 1110, ICOL
|
DO 230, IROW = 1, N
|
PRINT 1100, B1(IROW,ICOL)
|
230 CONTINUE
|
240 CONTINUE
|
C
|
C Refine the solution and print the refined solution with the
|
C error bounds.
|
C
|
CALL ZHPRFS ('UPPER TRIANGLE OF A STORED', N, NRHS, A, AF,
|
$ IPIVOT, B, LDB, X, LDX, FERR, BERR, WORK,
|
$ WORK2, INFO)
|
IF (INFO .LT. 0) THEN
|
PRINT 1120, ABS(INFO)
|
STOP 5
|
END IF
|
DO 310, ICOL = 1, NRHS
|
PRINT 1130, ICOL
|
DO 300, IROW = 1, N
|
PRINT 1100, X(IROW,ICOL)
|
300 CONTINUE
|
310 CONTINUE
|
DO 320, ICOL = 1, NRHS
|
CALL ZHPMV ('UPPER TRIANGLE OF A STORED', N,
|
$ (1.0D0,0.0D0), A, X(1,ICOL), 1, (0.0D0,0.0D0),
|
$ B1(1,ICOL), 1)
|
320 CONTINUE
|
DO 340, ICOL = 1, NRHS
|
PRINT 1140, ICOL
|
DO 330, IROW = 1, N
|
PRINT 1100, B1(IROW,ICOL)
|
330 CONTINUE
|
340 CONTINUE
|
DO 350, IROW = 1, NRHS
|
PRINT 1150, IROW, BERR(IROW)
|
PRINT 1160, IROW, FERR(IROW)
|
350 CONTINUE
|
C
|
1000 FORMAT (1X, 'A:')
|
1010 FORMAT (1X, 10(: 1X, '(', E8.2, ', ', E8.2, ')'))
|
1020 FORMAT (/1X, 'b for system #', I1)
|
1030 FORMAT (1X, 'Illegal argument to ZHPTRF, argument #', I2)
|
1040 FORMAT (1X, 'A is singular to working precision.')
|
1050 FORMAT (1X, 'Illegal argument to ZHPCON, argument #', I2)
|
1060 FORMAT (/1X,
|
$ 'Estimated reciprocal condition number of A: ',
|
$ E9.3)
|
1070 FORMAT (1X, 'A may be singular to working precision.')
|
1080 FORMAT (1X, 'Illegal argument to ZHPTRS, argument #', I2)
|
1090 FORMAT (/1X, 'Initial xb for system #', I1)
|
1100 FORMAT (3X, '(', E20.14, ', ', E20.14, ')')
|
1110 FORMAT (/1X, 'Ax with initial x for system #', I1)
|
1120 FORMAT (1X, 'Illegal argument to ZHPRFS, argument #', I2)
|
1130 FORMAT (/1X, 'Refined xb for system #', I1)
|
1140 FORMAT (/1X, 'Ax with refined x for system #', I1)
|
1150 FORMAT (/1X, 'Backward error for system #', I1, ': ', E12.6)
|
1160 FORMAT (1X, 'Forward error for system #', I1, ': ', E12.6)
|
C
|
END
|
|