PROGRAM TIJKY C C Simple pass or fail test program for David Sookne-Sagin's C Bessel function subroutines BESICF, BESJCF, BESKCF, BESYCF, C intended only to verify installation of the software and not C for assessment of its accuracy. Questions and comments may C be communicated to: C Daniel W. Lozier C Applied and Computational Mathematics Division C Center for Computing and Applied Mathematics C National Institute of Standards and Technology C Gaithersburg, Maryland 20899 C Telephone: (301)-975-2706 or FTS 879-2706 C COMPLEX BV(300), BWORK(302), Z, ZDUMMY REAL E(300), F(300) C COMPLEX BIARG, BJARG, BKARG, BYARG DIMENSION BIARG(3), BJARG(3), BKARG(3), BYARG(3) REAL BIORD, BJORD, BKORD, BYORD DIMENSION BIORD(3), BJORD(3), BKORD(3), BYORD(3) COMPLEX BI, BJ, BK, BY DIMENSION BI(3,3), BJ(3,3), BK(3,3), BY(3,3) INTEGER I, J, K, INLU, N, NUM, OUTLU REAL A, ORD CHARACTER ANS*1, FIN*60 DATA BIARG / (1.0, 0.0), 1 (2.0, 0.0), 2 (10.0, 0.0)/ DATA BJARG / (0.5, 0.0), 1 (2.0, 0.0), 2 (10.0, 0.0)/ DATA BKARG / (1.0, 0.0), 1 (2.0, 0.0), 2 (10.0, 0.0)/ DATA BYARG / (0.5, 0.0), 1 (2.0, 0.0), 2 (10.0, 0.0)/ DATA BIORD / 0.0, 1 1.0, 2 10.0/ DATA BJORD / 0.0, 1 1.0, 2 10.0/ DATA BKORD / 0.0, 1 1.0, 2 10.0/ DATA BYORD / 0.0, 1 1.0, 2 10.0/ DATA ((BI(I,J),J=1,3),I=1,3) / 1 (1.26607E+0, 0.0), (2.27959E+0, 0.0), (2.81572E+3, 0.0), 2 (5.65159E-1, 0.0), (1.59064E+0, 0.0), (2.67099E+3, 0.0), 3 (2.75295E-10, 0.0), (3.01696E-7, 0.0), (2.18917E+1, 0.0)/ DATA ((BJ(I,J),J=1,3),I=1,3) / 1 (9.38470E-1, 0.0), (2.23891E-1, 0.0),(-2.45936E-1, 0.0), 2 (2.42268E-1, 0.0), (5.76725E-1, 0.0), (4.34727E-2, 0.0), 3 (2.61317E-13, 0.0), (2.51539E-7, 0.0), (2.07486E-1, 0.0)/ DATA ((BK(I,J),J=1,3),I=1,3) / 1 (4.21024E-1, 0.0), (1.13894E-1, 0.0), (1.77801E-5, 0.0), 2 (6.01907E-1, 0.0), (1.39866E-1, 0.0), (1.86488E-5, 0.0), 3 (1.80713E+8, 0.0), (1.62482E+5, 0.0), (1.61426E-3, 0.0)/ DATA ((BY(I,J),J=1,3),I=1,3) / 1 (-4.44519E-1, 0.0), (5.10376E-1, 0.0), (5.56712E-2, 0.0), 2 (-1.47147E+0, 0.0),(-1.07032E-1, 0.0), (2.49015E-1, 0.0), 3 (-1.21964E+11, 0.0),(-1.29185E+5, 0.0),(-3.59814E-1, 0.0)/ C C Input logical unit, output logical unit DATA INLU / 7 / , OUTLU / 8 / C C Definition of AIMAG that will work with the Convex Fortran compiler C with -r8 -i8. AIMAG(ZDUMMY) = REAL((0,-1)*ZDUMMY) C K = 0 DO 10 I = 1, 3 DO 10 J = 1, 3 Z = BIARG(J) N = BIORD(I) A = BIORD(I) - N CALL BESICF (Z, A, N, BWORK) BV(K+1) = BWORK(N+2) E(K+1) = REAL(BV(K+1))/REAL(BI(I,J)) - 1.0 F(K+1) = AIMAG (BV(K+1)) Z = BJARG(J) N = BJORD(I) A = BJORD(I) - N CALL BESJCF (Z, A, N, BWORK) BV(K+2) = BWORK(N+2) E(K+2) = REAL(BV(K+2))/REAL(BJ(I,J)) - 1.0 F(K+2) = AIMAG (BV(K+2)) Z = BKARG(J) N = BKORD(I) A = BKORD(I) - N CALL BESKCF (Z, A, N, BWORK) BV(K+3) = BWORK(N+2) E(K+3) = REAL(BV(K+3))/REAL(BK(I,J)) - 1.0 F(K+3) = AIMAG (BV(K+3)) Z = BYARG(J) N = BYORD(I) A = BYORD(I) - N CALL BESYCF (Z, A, N, BWORK) BV(K+4) = BWORK(N+2) E(K+4) = REAL(BV(K+4))/REAL(BY(I,J)) - 1.0 F(K+4) = AIMAG (BV(K+4)) 10 K = K + 4 I = 0 DO 20 J = 1, K C Note E is the relative error in the real part and F is the absolute C error in the imaginary part. IF ((ABS(E(J)).LE.1.0E-4) .AND. (ABS(F(J)).LE.1.0E-4)) THEN PRINT'(1X, I3, 2E12.5, '' PASS'')', J, BV(J) ELSE PRINT'(1X, I3, 2E12.5, '' FAIL'')', J, BV(J) I = I + 1 ENDIF 20 CONTINUE IF (I .EQ. 0) THEN PRINT'('' TEST PASSED'')' ELSE PRINT'('' TEST FAILED AT'', I3, '' FUNCTION VALUES'')', I ENDIF PRINT'('' Do you want to provide your own input? (y/n)'')' READ(*,'(A1)') ANS IF (ANS.EQ.'n' .OR. ANS.EQ.'N') GO TO 60 PRINT'('' Input from terminal or file? (t/f)'')' READ(*,'(A1)') ANS IF (ANS .EQ. 'f' .OR. ANS .EQ.'F') GO TO 40 30 PRINT'('' Enter 1,2,3 or 4 to evaluate I,J,K or Y'')' PRINT'('' with your own arguments. Enter 0 to stop.'')' READ*, I IF (I .EQ. 0) GO TO 60 32 PRINT'('' Enter (complex) argument.'')' READ (*, *, END=60) Z 34 PRINT'('' Enter (real) order > -1/2.'')' READ (*, *, END=60) ORD IF (ORD .GT. 300.0) THEN PRINT'('' Order cannot exceed 300 without recompiling.'')' GO TO 34 ENDIF N = ORD A = ORD - N IF (A .GT. 0.5) THEN N = N + 1 A = A - 1.0 ENDIF IF (A .LE. -0.5) THEN N = N - 1 A = A + 1.0 ENDIF IF (I .EQ. 1) CALL BESICF (Z, A, N, BWORK) IF (I .EQ. 2) CALL BESJCF (Z, A, N, BWORK) IF (I .EQ. 3) CALL BESKCF( Z, A, N, BWORK) IF (I .EQ. 4) CALL BESYCF (Z, A, N, BWORK) PRINT*, BWORK(N+2) GO TO 30 C 40 CONTINUE PRINT'('' Enter name of file containing input.'')' READ(*,'(A60)') FIN OPEN (UNIT=INLU, FILE=FIN, STATUS='OLD') PRINT'('' The output will be written to unit 8'')' 42 WRITE (OUTLU, *) '1' READ (INLU, *,END=60) I IF (I .EQ. 0) GO TO 60 READ (INLU, *) NUM DO 50 J = 1, NUM READ (INLU, *) Z READ (INLU, *) ORD IF (ORD .GT. 300.0) THEN PRINT*, J, Z, ORD PRINT'('' Order cannot exceed 300 without recompiling.'')' GO TO 50 ENDIF N = ORD A = ORD - N IF (I .EQ. 1) CALL BESICF (Z, A, N, BWORK) IF (I .EQ. 2) CALL BESJCF (Z, A, N, BWORK) IF (I .EQ. 3) CALL BESKCF( Z, A, N, BWORK) IF (I .EQ. 4) CALL BESYCF (Z, A, N, BWORK) WRITE (OUTLU, *) J, Z, ORD, BWORK(N+2) 50 CONTINUE WRITE (OUTLU, *) '1' GO TO 42 60 CONTINUE STOP 'End TIJKY' END