C C --------------------------------------- C Compare scalar fnlib routines to vfnlib C --------------------------------------- C C Ron Boisvert and Bonita Saunders C Computing and Applied Mathematics Laboratory C National Institute of Standards and Technology C Gaithersburg, MD 20899 C C boisvert@cam.nist.gov C saunders@cam.nist.gov C C 21 Feb 91 C C---------------------------------------------------------------------- C INTEGER M PARAMETER ( M=2000 ) C REAL XLO, XHI, X, F, FX, WORK INTEGER I, IWORK, J C DIMENSION XLO(3,8), XHI(3,8), X(M), F(M), FX(M), WORK(7*M), + IWORK(M) C SAVE XLO, XHI C EXTERNAL VI0, BESI0, VI1, BESI1, VJ0, BESJ0, VJ1, BESJ1, + VK0, BESK0, VK1, BESK1, VY0, BESY0, VY1, BESY1 C DATA ((XLO(I,J),XHI(I,J),I=1,3),J=1,8) / + 0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0, + 0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0, + 0.0E0, 4.0E0 , -8.0E0, -4.0E0, 8.0E0, 50.0E0, + 0.0E0, 4.0E0 , -8.0E0, -4.0E0, 8.0E0, 50.0E0, + 0.0E0, 2.0E0 , 2.0E0, 8.0E0, 10.0E0, 50.0E0, + 0.0E0, 2.0E0 , 2.0E0, 8.0E0, 10.0E0, 50.0E0, + 0.0E0, 1.0E-3, 1.0E0, 4.0E0, 4.0E0, 50.0E0, + 0.0E0, 1.0E-3, 1.0E0, 4.0E0, 4.0E0, 50.0E0 / C C---------------------------------------------------------------------- C C ... TEST EACH USER-CALLABLE ROUTINE C CALL TEST('I0',BESI0,VI0,XLO(1,1),XHI(1,1),M,X,F,FX,WORK,IWORK) CALL TEST('I1',BESI1,VI1,XLO(1,2),XHI(1,2),M,X,F,FX,WORK,IWORK) CALL TEST('J0',BESJ0,VJ0,XLO(1,3),XHI(1,3),M,X,F,FX,WORK,IWORK) CALL TEST('J1',BESJ1,VJ1,XLO(1,4),XHI(1,4),M,X,F,FX,WORK,IWORK) CALL TEST('K0',BESK0,VK0,XLO(1,5),XHI(1,5),M,X,F,FX,WORK,IWORK) CALL TEST('K1',BESK1,VK1,XLO(1,6),XHI(1,6),M,X,F,FX,WORK,IWORK) CALL TEST('Y0',BESY0,VY0,XLO(1,7),XHI(1,7),M,X,F,FX,WORK,IWORK) CALL TEST('Y1',BESY1,VY1,XLO(1,8),XHI(1,8),M,X,F,FX,WORK,IWORK) C END SUBROUTINE TEST (NAME, SFUN, VSUB, XLOW, XHIGH, M, X, F, FX, + WORK, IWORK) C C---------------------------------------------------------------------- C C compare scalar function sfun to vectorized subroutine vsub for m C arguments distributed in various ways in the intervals C (xlo(i),xhi(i)), i=1,2,3. C C---------------------------------------------------------------------- C C ---------- C PARAMETERS C ---------- C REAL SFUN, XLOW, XHIGH, X, F, FX, WORK INTEGER IWORK, M CHARACTER*(*) NAME C DIMENSION XLOW(3), XHIGH(3), X(M), F(M), FX(M), WORK(6*M), + IWORK(M) C EXTERNAL SFUN, VSUB C C --------------- C LOCAL VARIABLES C --------------- C LOGICAL DEBUG INTEGER NREP, NTRY PARAMETER ( DEBUG=.FALSE., NREP=15, NTRY=19 ) C INTEGER I, INFO, IREP, IOFF, IPC(NTRY), JPC(NTRY), K, KPC, + MA, MB, MC, NFAIL, NREPV REAL DIFF, DMAX, DX, EPS, R1MACH REAL SECOND, SU, SUMAX, SUMIN, SUAVE, T0, T1, TS, TSAVE, TV, + TVAVE C SAVE IPC, JPC C DATA IPC / 0, 0, 0, 0, 0, 25, 25, 25, 25, 33, 50, 50, 50, + 75, 75,100, 1, 1, 98/ DATA JPC / 0, 25, 50, 75,100, 0, 25, 50, 75, 33, 0, 25, 50, + 0, 25, 0, 1, 98, 1/ C C---------------------------------------------------------------------- C C ntry = number of trials (each with different argument distribution) C nrep = number of times test repeated to get reliable timings C C the distribution of arguments in trial k is C C ipc(k)% in (xlow(1),xhigh(1)) C jpc(k)% in (xlow(2),xhigh(2)) C 100-ipc(k)-jpc(k)% in (xlow(3),xhigh(3)) C C scalar results (fx) and vector results (f) are compared and accepted C if abs((f - fx)/fx) < eps for fx > 1 (relative difference) C or abs(f - fx) < eps for fx < 1 (absolute difference) C C---------------------------------------------------------------------- C EPS = 5.0E0*R1MACH(4) C ... TRIGGER INITIALIZATIONS C X(1) = 1.0E0 CALL VSUB(1,X,F,WORK,IWORK,INFO) FX(1) = SFUN(X(1)) C PRINT * PRINT * PRINT *,' ----------' PRINT *,' TEST OF ',NAME PRINT *,' ----------' PRINT * PRINT *,' Vector Length = ',M PRINT *,' Repetitions = ',NREP PRINT *,' Eps = ',EPS PRINT * PRINT *,' Range A = (',XLOW(1),',',XHIGH(1),')' PRINT *,' Range B = (',XLOW(2),',',XHIGH(2),')' PRINT *,' Range C = (',XLOW(3),',',XHIGH(3),')' PRINT * PRINT *,' %A = percent arguments in range A' PRINT *,' %B = percent arguments in range B' PRINT *,' %C = percent arguments in range C' PRINT *,' Scalar = time for scalar code (FNLIB)' PRINT *,' Vector = time for vector code (VFNLIB)' PRINT *,' SU = speedup (Scalar/Vector)' PRINT *,' Nfail = number of scalar-vector differences .gt. eps' PRINT *,' Dmax = maximum scalar-vector difference' PRINT * PRINT *,' All times in seconds/argument.' PRINT * PRINT *,' Scalar-vector differences measured by' PRINT *,' abs(s-v)/max(abs(s),1)' PRINT *,' where s,v are scalar,vector results.' PRINT * PRINT *,' %A %B %C Scalar Vector SU Nfail Dmax' PRINT *,' -------------------------------------------------------' C SUMAX = 0.0 SUMIN = 1000.0 SUAVE = 0.0 TSAVE = 0.0 TVAVE = 0.0 C C ... PERFORM EACH TEST C DO 100 K=1,NTRY C MA = IPC(K)*M/100 MB = JPC(K)*M/100 KPC = 100 - IPC(K) - JPC(K) MC = M - MA - MB C C ... SELECT X VECTOR C IOFF = 0 C IF (MA .GT. 0) THEN DX = (XHIGH(1) - XLOW(1))/MA DO 10 I=1,MA X(I) = XLOW(1) + I*DX 10 CONTINUE IOFF = IOFF + MA ENDIF C IF (MB .GT. 0) THEN DX = (XHIGH(2) - XLOW(2))/MB DO 20 I=1,MB X(I+IOFF) = XLOW(2) + I*DX 20 CONTINUE IOFF = IOFF + MB ENDIF C IF (MC .GT. 0) THEN DX = (XHIGH(3) - XLOW(3))/MC DO 30 I=1,MC X(I+IOFF) = XLOW(3) + I*DX 30 CONTINUE ENDIF C C ... RUN VECTOR CODE C NREPV = 5*NREP T0 = SECOND() DO 40 IREP=1,NREPV CALL VSUB(M,X,F,WORK,IWORK,INFO) 40 CONTINUE T1 = SECOND() TV = (T1 - T0)/NREPV IF (INFO .LT. 0) PRINT *,'Warning in ',NAME,' : Info = ',INFO IF (INFO .GT. 0) PRINT *,'Failure in ',NAME,' : Info = ',INFO, + ' Position = ',IWORK(1) C C ... RUN SCALAR CODE C T0 = SECOND() DO 60 IREP=1,NREP DO 50 I=1,M FX(I) = SFUN(X(I)) 50 CONTINUE 60 CONTINUE T1 = SECOND() TS = (T1 - T0)/NREP C C ... COMPUTE TIMINGS C SU = TS/TV SUMIN = MIN(SU,SUMIN) SUMAX = MAX(SU,SUMAX) SUAVE = SUAVE + SU TSAVE = TSAVE + TS TVAVE = TVAVE + TV C C ... EVALUATE RESULTS C DMAX = 0.0E0 NFAIL = 0 DO 70 I=1,M DIFF = ABS(FX(I)-F(I))/MAX(ABS(FX(I)),1.0E0) DMAX = MAX(DIFF,DMAX) IF (DIFF .GT. EPS) THEN NFAIL = NFAIL + 1 IF (DEBUG) PRINT *,'Failure ',I,X(I),F(I),FX(I),DIFF ENDIF 70 CONTINUE C C ... PRINT RESULTS C PRINT 1000, IPC(K),JPC(K),KPC,TS/M,TV/M,SU,NFAIL,DMAX C 100 CONTINUE C SUAVE = SUAVE/NTRY TSAVE = TSAVE/M/NTRY TVAVE = TVAVE/M/NTRY C PRINT * PRINT 1002,'Ave Scalar = ',TSAVE,' sec/argument' PRINT 1002,'Ave Vector = ',TVAVE,' sec/argument' PRINT * PRINT 1001,'Max speed up = ',SUMAX PRINT 1001,'Min speed up = ',SUMIN PRINT * PRINT 1001,'Ave speed up = ',SUAVE C RETURN 1000 FORMAT(3(1X,I3),1P,2(3X,E7.1),0P,1X,F5.1,1X,I4,3X,1P,E7.1) 1001 FORMAT(1X,A,F5.1,A) 1002 FORMAT(1X,A,1P,E7.1,A) END function besi0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bi0cs(12) c external alog, besi0e, csevl, exp, inits, r1mach, sqrt c c series for bi0 on the interval 0. to 9.00000d+00 c with weighted error 2.46e-18 c log weighted error 17.61 c significant figures required 17.90 c decimal places required 18.15 c data bi0 cs( 1) / -.0766054725 2839144951e0 / data bi0 cs( 2) / 1.9273379539 93808270e0 / data bi0 cs( 3) / .2282644586 920301339e0 / data bi0 cs( 4) / .0130489146 6707290428e0 / data bi0 cs( 5) / .0004344270 9008164874e0 / data bi0 cs( 6) / .0000094226 5768600193e0 / data bi0 cs( 7) / .0000001434 0062895106e0 / data bi0 cs( 8) / .0000000016 1384906966e0 / data bi0 cs( 9) / .0000000000 1396650044e0 / data bi0 cs(10) / .0000000000 0009579451e0 / data bi0 cs(11) / .0000000000 0000053339e0 / data bi0 cs(12) / .0000000000 0000000245e0 / c data nti0, xsml, xmax / 0, 0., 0. / c if (nti0.ne.0) go to 10 nti0 = inits (bi0cs, 12, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) xmax = alog (r1mach(2)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi0 = 1.0 if (y.gt.xsml) besi0 = 2.75 + csevl (y*y/4.5-1.0, bi0cs, nti0) return c 20 if (y.gt.xmax) call seteru (34hbesi0 abs(x) so big i0 overflows, 1 34, 1, 2) c besi0 = exp(y) * besi0e(x) c return end function besi1 (x) c oct 1983 version. w. fullerton, c3, los alamos scientific lab. dimension bi1cs(11) c external alog, besi1e, csevl, exp, inits, r1mach, sqrt c c series for bi1 on the interval 0. to 9.00000d+00 c with weighted error 2.40e-17 c log weighted error 16.62 c significant figures required 16.23 c decimal places required 17.14 c data bi1 cs( 1) / -.0019717132 61099859e0 / data bi1 cs( 2) / .4073488766 7546481e0 / data bi1 cs( 3) / .0348389942 99959456e0 / data bi1 cs( 4) / .0015453945 56300123e0 / data bi1 cs( 5) / .0000418885 21098377e0 / data bi1 cs( 6) / .0000007649 02676483e0 / data bi1 cs( 7) / .0000000100 42493924e0 / data bi1 cs( 8) / .0000000000 99322077e0 / data bi1 cs( 9) / .0000000000 00766380e0 / data bi1 cs(10) / .0000000000 00004741e0 / data bi1 cs(11) / .0000000000 00000024e0 / c c data nti1, xmin, xsml, xmax / 0, 3*0. / c if (nti1.ne.0) go to 10 nti1 = inits (bi1cs, 11, 0.1*r1mach(3)) xmin = 2.0*r1mach(1) xsml = sqrt (8.0*r1mach(3)) xmax = alog (r1mach(2)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi1 = 0.0 if (y.eq.0.0) return c if (y.le.xmin) call seteru ( 1 37hbesi1 abs(x) so small i1 underflows, 37, 1, 0) if (y.gt.xmin) besi1 = 0.5*x if (y.gt.xsml) besi1 = x * (.875 + csevl(y*y/4.5-1., bi1cs, nti1)) return c 20 if (y.gt.xmax) call seteru ( 1 34hbesi1 abs(x) so big i1 overflows, 34, 2, 2) c besi1 = exp(y) * besi1e(x) c return end function besj0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bj0cs(13), bm0cs(21), bth0cs(24) c external cos, csevl, inits, r1mach, sqrt c c series for bj0 on the interval 0. to 1.60000d+01 c with weighted error 7.47e-18 c log weighted error 17.13 c significant figures required 16.98 c decimal places required 17.68 c data bj0 cs( 1) / .1002541619 68939137e0 / data bj0 cs( 2) / -.6652230077 64405132e0 / data bj0 cs( 3) / .2489837034 98281314e0 / data bj0 cs( 4) / -.0332527231 700357697e0 / data bj0 cs( 5) / .0023114179 304694015e0 / data bj0 cs( 6) / -.0000991127 741995080e0 / data bj0 cs( 7) / .0000028916 708643998e0 / data bj0 cs( 8) / -.0000000612 108586630e0 / data bj0 cs( 9) / .0000000009 838650793e0 / data bj0 cs(10) / -.0000000000 124235515e0 / data bj0 cs(11) / .0000000000 001265433e0 / data bj0 cs(12) / -.0000000000 000010619e0 / data bj0 cs(13) / .0000000000 000000074e0 / c c series for bm0 on the interval 0. to 6.25000d-02 c with weighted error 4.98e-17 c log weighted error 16.30 c significant figures required 14.97 c decimal places required 16.96 c data bm0 cs( 1) / .0928496163 7381644e0 / data bm0 cs( 2) / -.0014298770 7403484e0 / data bm0 cs( 3) / .0000283057 9271257e0 / data bm0 cs( 4) / -.0000014330 0611424e0 / data bm0 cs( 5) / .0000001202 8628046e0 / data bm0 cs( 6) / -.0000000139 7113013e0 / data bm0 cs( 7) / .0000000020 4076188e0 / data bm0 cs( 8) / -.0000000003 5399669e0 / data bm0 cs( 9) / .0000000000 7024759e0 / data bm0 cs(10) / -.0000000000 1554107e0 / data bm0 cs(11) / .0000000000 0376226e0 / data bm0 cs(12) / -.0000000000 0098282e0 / data bm0 cs(13) / .0000000000 0027408e0 / data bm0 cs(14) / -.0000000000 0008091e0 / data bm0 cs(15) / .0000000000 0002511e0 / data bm0 cs(16) / -.0000000000 0000814e0 / data bm0 cs(17) / .0000000000 0000275e0 / data bm0 cs(18) / -.0000000000 0000096e0 / data bm0 cs(19) / .0000000000 0000034e0 / data bm0 cs(20) / -.0000000000 0000012e0 / data bm0 cs(21) / .0000000000 0000004e0 / c c series for bth0 on the interval 0. to 6.25000d-02 c with weighted error 3.67e-17 c log weighted error 16.44 c significant figures required 15.53 c decimal places required 17.13 c data bth0cs( 1) / -.2463916377 4300119e0 / data bth0cs( 2) / .0017370983 07508963e0 / data bth0cs( 3) / -.0000621836 33402968e0 / data bth0cs( 4) / .0000043680 50165742e0 / data bth0cs( 5) / -.0000004560 93019869e0 / data bth0cs( 6) / .0000000621 97400101e0 / data bth0cs( 7) / -.0000000103 00442889e0 / data bth0cs( 8) / .0000000019 79526776e0 / data bth0cs( 9) / -.0000000004 28198396e0 / data bth0cs(10) / .0000000001 02035840e0 / data bth0cs(11) / -.0000000000 26363898e0 / data bth0cs(12) / .0000000000 07297935e0 / data bth0cs(13) / -.0000000000 02144188e0 / data bth0cs(14) / .0000000000 00663693e0 / data bth0cs(15) / -.0000000000 00215126e0 / data bth0cs(16) / .0000000000 00072659e0 / data bth0cs(17) / -.0000000000 00025465e0 / data bth0cs(18) / .0000000000 00009229e0 / data bth0cs(19) / -.0000000000 00003448e0 / data bth0cs(20) / .0000000000 00001325e0 / data bth0cs(21) / -.0000000000 00000522e0 / data bth0cs(22) / .0000000000 00000210e0 / data bth0cs(23) / -.0000000000 00000087e0 / data bth0cs(24) / .0000000000 00000036e0 / c data pi4 / 0.7853981633 9744831e0 / data ntj0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./ c if (ntj0.ne.0) go to 10 ntj0 = inits (bj0cs, 13, 0.1*r1mach(3)) ntm0 = inits (bm0cs, 21, 0.1*r1mach(3)) ntth0 = inits (bth0cs, 24, 0.1*r1mach(3)) c xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 y = abs(x) if (y.gt.4.0) go to 20 c besj0 = 1.0 if (y.gt.xsml) besj0 = csevl (.125*y*y-1., bj0cs, ntj0) return c 20 if (y.gt.xmax) call seteru ( 1 42hbesj0 no precision because abs(x) is big, 42, 1, 2) c z = 32.0/y**2 - 1.0 ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(y) theta = y - pi4 + csevl (z, bth0cs, ntth0) / y besj0 = ampl * cos (theta) c return end function besj1 (x) c sept 1983 edition. w. fullerton, c3, los alamos scientific lab. dimension bj1cs(12), bm1cs(21), bth1cs(24) c external cos, csevl, inits, r1mach, sqrt c c series for bj1 on the interval 0. to 1.60000d+01 c with weighted error 4.48e-17 c log weighted error 16.35 c significant figures required 15.77 c decimal places required 16.89 c data bj1 cs( 1) / -.1172614151 3332787e0 / data bj1 cs( 2) / -.2536152183 0790640e0 / data bj1 cs( 3) / .0501270809 84469569e0 / data bj1 cs( 4) / -.0046315148 09625081e0 / data bj1 cs( 5) / .0002479962 29415914e0 / data bj1 cs( 6) / -.0000086789 48686278e0 / data bj1 cs( 7) / .0000002142 93917143e0 / data bj1 cs( 8) / -.0000000039 36093079e0 / data bj1 cs( 9) / .0000000000 55911823e0 / data bj1 cs(10) / -.0000000000 00632761e0 / data bj1 cs(11) / .0000000000 00005840e0 / data bj1 cs(12) / -.0000000000 00000044e0 / c c series for bm1 on the interval 0. to 6.25000d-02 c with weighted error 5.61e-17 c log weighted error 16.25 c significant figures required 14.97 c decimal places required 16.91 c data bm1 cs( 1) / .1047362510 931285e0 / data bm1 cs( 2) / .0044244389 3702345e0 / data bm1 cs( 3) / -.0000566163 9504035e0 / data bm1 cs( 4) / .0000023134 9417339e0 / data bm1 cs( 5) / -.0000001737 7182007e0 / data bm1 cs( 6) / .0000000189 3209930e0 / data bm1 cs( 7) / -.0000000026 5416023e0 / data bm1 cs( 8) / .0000000004 4740209e0 / data bm1 cs( 9) / -.0000000000 8691795e0 / data bm1 cs(10) / .0000000000 1891492e0 / data bm1 cs(11) / -.0000000000 0451884e0 / data bm1 cs(12) / .0000000000 0116765e0 / data bm1 cs(13) / -.0000000000 0032265e0 / data bm1 cs(14) / .0000000000 0009450e0 / data bm1 cs(15) / -.0000000000 0002913e0 / data bm1 cs(16) / .0000000000 0000939e0 / data bm1 cs(17) / -.0000000000 0000315e0 / data bm1 cs(18) / .0000000000 0000109e0 / data bm1 cs(19) / -.0000000000 0000039e0 / data bm1 cs(20) / .0000000000 0000014e0 / data bm1 cs(21) / -.0000000000 0000005e0 / c c series for bth1 on the interval 0. to 6.25000d-02 c with weighted error 4.10e-17 c log weighted error 16.39 c significant figures required 15.96 c decimal places required 17.08 c data bth1cs( 1) / .7406014102 6313850e0 / data bth1cs( 2) / -.0045717556 59637690e0 / data bth1cs( 3) / .0001198185 10964326e0 / data bth1cs( 4) / -.0000069645 61891648e0 / data bth1cs( 5) / .0000006554 95621447e0 / data bth1cs( 6) / -.0000000840 66228945e0 / data bth1cs( 7) / .0000000133 76886564e0 / data bth1cs( 8) / -.0000000024 99565654e0 / data bth1cs( 9) / .0000000005 29495100e0 / data bth1cs(10) / -.0000000001 24135944e0 / data bth1cs(11) / .0000000000 31656485e0 / data bth1cs(12) / -.0000000000 08668640e0 / data bth1cs(13) / .0000000000 02523758e0 / data bth1cs(14) / -.0000000000 00775085e0 / data bth1cs(15) / .0000000000 00249527e0 / data bth1cs(16) / -.0000000000 00083773e0 / data bth1cs(17) / .0000000000 00029205e0 / data bth1cs(18) / -.0000000000 00010534e0 / data bth1cs(19) / .0000000000 00003919e0 / data bth1cs(20) / -.0000000000 00001500e0 / data bth1cs(21) / .0000000000 00000589e0 / data bth1cs(22) / -.0000000000 00000237e0 / data bth1cs(23) / .0000000000 00000097e0 / data bth1cs(24) / -.0000000000 00000040e0 / c c data pi4 / 0.7853981633 9744831e0 / data ntj1, ntm1, ntth1, xsml, xmin, xmax / 3*0, 3*0./ c if (ntj1.ne.0) go to 10 ntj1 = inits (bj1cs, 12, 0.1*r1mach(3)) ntm1 = inits (bm1cs, 21, 0.1*r1mach(3)) ntth1 = inits (bth1cs, 24, 0.1*r1mach(3)) c xsml = sqrt (8.0*r1mach(3)) xmin = 2.0*r1mach(1) xmax = 1.0/r1mach(4) c 10 y = abs(x) if (y.gt.4.0) go to 20 c besj1 = 0.0 if (y.eq.0.0) return if (y.le.xmin) call seteru ( 1 37hbesj1 abs(x) so small j1 underflows, 37, 1, 0) if (y.gt.xmin) besj1 = 0.5*x if (y.gt.xsml) besj1 = x * (.25 + csevl(.125*y*y-1., bj1cs, ntj1)) return c 20 if (y.gt.xmax) call seteru ( 1 42hbesj1 no precision because abs(x) is big, 42, 2, 2) z = 32.0/y**2 - 1.0 ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(y) theta = y - 3.0*pi4 + csevl (z, bth1cs, ntth1) / y besj1 = sign (ampl, x) * cos (theta) c return end function besk0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bk0cs(11) c external alog, besi0, besk0e, csevl, exp, inits, r1mach, sqrt c c series for bk0 on the interval 0. to 4.00000d+00 c with weighted error 3.57e-19 c log weighted error 18.45 c significant figures required 17.99 c decimal places required 18.97 c data bk0 cs( 1) / -.0353273932 3390276872e0 / data bk0 cs( 2) / .3442898999 246284869e0 / data bk0 cs( 3) / .0359799365 1536150163e0 / data bk0 cs( 4) / .0012646154 1144692592e0 / data bk0 cs( 5) / .0000228621 2103119451e0 / data bk0 cs( 6) / .0000002534 7910790261e0 / data bk0 cs( 7) / .0000000019 0451637722e0 / data bk0 cs( 8) / .0000000000 1034969525e0 / data bk0 cs( 9) / .0000000000 0004259816e0 / data bk0 cs(10) / .0000000000 0000013744e0 / data bk0 cs(11) / .0000000000 0000000035e0 / c data ntk0, xsml, xmax / 0, 0., 0. / c if (ntk0.ne.0) go to 10 ntk0 = inits (bk0cs, 11, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) xmax = -alog(r1mach(1)) xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01 c 10 if (x.le.0.) call seteru (29hbesk0 x is zero or negative, 29, 1 2, 2) if (x.gt.2.) go to 20 c y = 0. if (x.gt.xsml) y = x*x besk0 = -alog(0.5*x)*besi0(x) - .25 + csevl (.5*y-1., bk0cs, ntk0) return c 20 besk0 = 0. if (x.gt.xmax) call seteru (30hbesk0 x so big k0 underflows, 30, 1 1, 0) if (x.gt.xmax) return c besk0 = exp(-x) * besk0e(x) c return end function besk1 (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk1cs(11) c external alog, besi1, besk1e, csevl, exp, inits, r1mach, sqrt c c series for bk1 on the interval 0. to 4.00000d+00 c with weighted error 7.02e-18 c log weighted error 17.15 c significant figures required 16.73 c decimal places required 17.67 c data bk1 cs( 1) / .0253002273 389477705e0 / data bk1 cs( 2) / -.3531559607 76544876e0 / data bk1 cs( 3) / -.1226111808 22657148e0 / data bk1 cs( 4) / -.0069757238 596398643e0 / data bk1 cs( 5) / -.0001730288 957513052e0 / data bk1 cs( 6) / -.0000024334 061415659e0 / data bk1 cs( 7) / -.0000000221 338763073e0 / data bk1 cs( 8) / -.0000000001 411488392e0 / data bk1 cs( 9) / -.0000000000 006666901e0 / data bk1 cs(10) / -.0000000000 000024274e0 / data bk1 cs(11) / -.0000000000 000000070e0 / c data ntk1, xmin, xsml, xmax /0, 0., 0., 0. / c if (ntk1.ne.0) go to 10 ntk1 = inits (bk1cs, 11, 0.1*r1mach(3)) xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01) xsml = sqrt (4.0*r1mach(3)) xmax = -alog(r1mach(1)) xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01 c 10 if (x.le.0.) call seteru (29hbesk1 x is zero or negative, 29, 1 2, 2) if (x.gt.2.0) go to 20 c if (x.lt.xmin) call seteru (31hbesk1 x so small k1 overflows, 1 31, 3, 2) y = 0. if (x.gt.xsml) y = x*x besk1 = alog(0.5*x)*besi1(x) + 1 (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x return c 20 besk1 = 0. if (x.gt.xmax) call seteru (30hbesk1 x so big k1 underflows, 1 30, 1, 0) if (x.gt.xmax) return c besk1 = exp(-x) * besk1e(x) c return end function besy0 (x) c august 1980 version. w. fullerton, c3, los alamos scientific lab. dimension by0cs(13), bm0cs(21), bth0cs(24) c external alog, besj0, csevl, inits, r1mach, sin, sqrt c c series for by0 on the interval 0. to 1.60000d+01 c with weighted error 1.20e-17 c log weighted error 16.92 c significant figures required 16.15 c decimal places required 17.48 c data by0 cs( 1) / -.0112778393 92865573e0 / data by0 cs( 2) / -.1283452375 6042035e0 / data by0 cs( 3) / -.1043788479 9794249e0 / data by0 cs( 4) / .0236627491 83969695e0 / data by0 cs( 5) / -.0020903916 47700486e0 / data by0 cs( 6) / .0001039754 53939057e0 / data by0 cs( 7) / -.0000033697 47162423e0 / data by0 cs( 8) / .0000000772 93842676e0 / data by0 cs( 9) / -.0000000013 24976772e0 / data by0 cs(10) / .0000000000 17648232e0 / data by0 cs(11) / -.0000000000 00188105e0 / data by0 cs(12) / .0000000000 00001641e0 / data by0 cs(13) / -.0000000000 00000011e0 / c c series for bm0 on the interval 0. to 6.25000d-02 c with weighted error 4.98e-17 c log weighted error 16.30 c significant figures required 14.97 c decimal places required 16.96 c data bm0 cs( 1) / .0928496163 7381644e0 / data bm0 cs( 2) / -.0014298770 7403484e0 / data bm0 cs( 3) / .0000283057 9271257e0 / data bm0 cs( 4) / -.0000014330 0611424e0 / data bm0 cs( 5) / .0000001202 8628046e0 / data bm0 cs( 6) / -.0000000139 7113013e0 / data bm0 cs( 7) / .0000000020 4076188e0 / data bm0 cs( 8) / -.0000000003 5399669e0 / data bm0 cs( 9) / .0000000000 7024759e0 / data bm0 cs(10) / -.0000000000 1554107e0 / data bm0 cs(11) / .0000000000 0376226e0 / data bm0 cs(12) / -.0000000000 0098282e0 / data bm0 cs(13) / .0000000000 0027408e0 / data bm0 cs(14) / -.0000000000 0008091e0 / data bm0 cs(15) / .0000000000 0002511e0 / data bm0 cs(16) / -.0000000000 0000814e0 / data bm0 cs(17) / .0000000000 0000275e0 / data bm0 cs(18) / -.0000000000 0000096e0 / data bm0 cs(19) / .0000000000 0000034e0 / data bm0 cs(20) / -.0000000000 0000012e0 / data bm0 cs(21) / .0000000000 0000004e0 / c c series for bth0 on the interval 0. to 6.25000d-02 c with weighted error 3.67e-17 c log weighted error 16.44 c significant figures required 15.53 c decimal places required 17.13 c data bth0cs( 1) / -.2463916377 4300119e0 / data bth0cs( 2) / .0017370983 07508963e0 / data bth0cs( 3) / -.0000621836 33402968e0 / data bth0cs( 4) / .0000043680 50165742e0 / data bth0cs( 5) / -.0000004560 93019869e0 / data bth0cs( 6) / .0000000621 97400101e0 / data bth0cs( 7) / -.0000000103 00442889e0 / data bth0cs( 8) / .0000000019 79526776e0 / data bth0cs( 9) / -.0000000004 28198396e0 / data bth0cs(10) / .0000000001 02035840e0 / data bth0cs(11) / -.0000000000 26363898e0 / data bth0cs(12) / .0000000000 07297935e0 / data bth0cs(13) / -.0000000000 02144188e0 / data bth0cs(14) / .0000000000 00663693e0 / data bth0cs(15) / -.0000000000 00215126e0 / data bth0cs(16) / .0000000000 00072659e0 / data bth0cs(17) / -.0000000000 00025465e0 / data bth0cs(18) / .0000000000 00009229e0 / data bth0cs(19) / -.0000000000 00003448e0 / data bth0cs(20) / .0000000000 00001325e0 / data bth0cs(21) / -.0000000000 00000522e0 / data bth0cs(22) / .0000000000 00000210e0 / data bth0cs(23) / -.0000000000 00000087e0 / data bth0cs(24) / .0000000000 00000036e0 / c c twodpi = 2.0/pi data twodpi / 0.6366197723 6758134e0 / data pi4 / 0.7853981633 9744831e0 / data alnhaf / -0.6931471805 59945309e0 / data nty0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./ c if (nty0.ne.0) go to 10 nty0 = inits (by0cs, 13, 0.1*r1mach(3)) ntm0 = inits (bm0cs, 21, 0.1*r1mach(3)) ntth0 = inits (bth0cs, 24, 0.1*r1mach(3)) c xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 if (x.le.0.) call seteru (29hbesy0 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0) go to 20 c y = 0. if (x.gt.xsml) y = x*x besy0 = twodpi*(alnhaf+alog(x))*besj0(x) + .375 + 1 csevl (.125*y-1., by0cs, nty0) return c 20 if (x.gt.xmax) call seteru ( 1 37hbesy0 no precision because x is big, 37, 2, 2) c z = 32.0/x**2 - 1.0 ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(x) theta = x - pi4 + csevl (z, bth0cs, ntth0) / x besy0 = ampl * sin (theta) c return end function besy1 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension by1cs(14), bm1cs(21), bth1cs(24) c external alog, besj1, csevl, exp, inits, r1mach, sin, sqrt c c series for by1 on the interval 0. to 1.60000d+01 c with weighted error 1.87e-18 c log weighted error 17.73 c significant figures required 17.83 c decimal places required 18.30 c data by1 cs( 1) / .0320804710 0611908629e0 / data by1 cs( 2) / 1.2627078974 33500450e0 / data by1 cs( 3) / .0064999618 9992317500e0 / data by1 cs( 4) / -.0893616452 8860504117e0 / data by1 cs( 5) / .0132508812 2175709545e0 / data by1 cs( 6) / -.0008979059 1196483523e0 / data by1 cs( 7) / .0000364736 1487958306e0 / data by1 cs( 8) / -.0000010013 7438166600e0 / data by1 cs( 9) / .0000000199 4539657390e0 / data by1 cs(10) / -.0000000003 0230656018e0 / data by1 cs(11) / .0000000000 0360987815e0 / data by1 cs(12) / -.0000000000 0003487488e0 / data by1 cs(13) / .0000000000 0000027838e0 / data by1 cs(14) / -.0000000000 0000000186e0 / c c series for bm1 on the interval 0. to 6.25000d-02 c with weighted error 5.61e-17 c log weighted error 16.25 c significant figures required 14.97 c decimal places required 16.91 c data bm1 cs( 1) / .1047362510 931285e0 / data bm1 cs( 2) / .0044244389 3702345e0 / data bm1 cs( 3) / -.0000566163 9504035e0 / data bm1 cs( 4) / .0000023134 9417339e0 / data bm1 cs( 5) / -.0000001737 7182007e0 / data bm1 cs( 6) / .0000000189 3209930e0 / data bm1 cs( 7) / -.0000000026 5416023e0 / data bm1 cs( 8) / .0000000004 4740209e0 / data bm1 cs( 9) / -.0000000000 8691795e0 / data bm1 cs(10) / .0000000000 1891492e0 / data bm1 cs(11) / -.0000000000 0451884e0 / data bm1 cs(12) / .0000000000 0116765e0 / data bm1 cs(13) / -.0000000000 0032265e0 / data bm1 cs(14) / .0000000000 0009450e0 / data bm1 cs(15) / -.0000000000 0002913e0 / data bm1 cs(16) / .0000000000 0000939e0 / data bm1 cs(17) / -.0000000000 0000315e0 / data bm1 cs(18) / .0000000000 0000109e0 / data bm1 cs(19) / -.0000000000 0000039e0 / data bm1 cs(20) / .0000000000 0000014e0 / data bm1 cs(21) / -.0000000000 0000005e0 / c c series for bth1 on the interval 0. to 6.25000d-02 c with weighted error 4.10e-17 c log weighted error 16.39 c significant figures required 15.96 c decimal places required 17.08 c data bth1cs( 1) / .7406014102 6313850e0 / data bth1cs( 2) / -.0045717556 59637690e0 / data bth1cs( 3) / .0001198185 10964326e0 / data bth1cs( 4) / -.0000069645 61891648e0 / data bth1cs( 5) / .0000006554 95621447e0 / data bth1cs( 6) / -.0000000840 66228945e0 / data bth1cs( 7) / .0000000133 76886564e0 / data bth1cs( 8) / -.0000000024 99565654e0 / data bth1cs( 9) / .0000000005 29495100e0 / data bth1cs(10) / -.0000000001 24135944e0 / data bth1cs(11) / .0000000000 31656485e0 / data bth1cs(12) / -.0000000000 08668640e0 / data bth1cs(13) / .0000000000 02523758e0 / data bth1cs(14) / -.0000000000 00775085e0 / data bth1cs(15) / .0000000000 00249527e0 / data bth1cs(16) / -.0000000000 00083773e0 / data bth1cs(17) / .0000000000 00029205e0 / data bth1cs(18) / -.0000000000 00010534e0 / data bth1cs(19) / .0000000000 00003919e0 / data bth1cs(20) / -.0000000000 00001500e0 / data bth1cs(21) / .0000000000 00000589e0 / data bth1cs(22) / -.0000000000 00000237e0 / data bth1cs(23) / .0000000000 00000097e0 / data bth1cs(24) / -.0000000000 00000040e0 / c c twodpi = 2.0/pi data twodpi / 0.6366197723 6758134e0 / data pi4 / 0.7853981633 9744831e0 / data nty1, ntm1, ntth1, xmin, xsml, xmax / 3*0, 3*0./ c if (nty1.ne.0) go to 10 nty1 = inits (by1cs, 14, 0.1*r1mach(3)) ntm1 = inits (bm1cs, 21, 0.1*r1mach(3)) ntth1 = inits (bth1cs, 24, 0.1*r1mach(3)) c xmin = 1.571*exp ( amax1(alog(r1mach(1)), -alog(r1mach(2)))+.01) xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 if (x.le.0.) call seteru (29hbesy1 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0) go to 20 c if (x.lt.xmin) call seteru (31hbesy1 x so small y1 overflows, 1 31, 3, 2) y = 0. if (x.gt.xsml) y = x*x besy1 = twodpi*alog(0.5*x)*besj1(x) + 1 (0.5 + csevl (.125*y-1., by1cs, nty1))/x return c 20 if (x.gt.xmax) call seteru ( 1 37hbesy1 no precision because x is big, 37, 2, 2) c z = 32.0/x**2 - 1.0 ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(x) theta = x - 3.0*pi4 + csevl (z, bth1cs, ntth1) / x besy1 = ampl * sin (theta) c return end function inits (os, nos, eta) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c initialize the orthogonal series so that inits is the number of terms c needed to insure the error is no larger than eta. ordinarily, eta c will be chosen to be one-tenth machine precision. c c input arguments -- c os array of nos coefficients in an orthogonal series. c nos number of coefficients in os. c eta requested accuracy of series. c dimension os(nos) c if (nos.lt.1) call seteru ( 1 35hinits number of coefficients lt 1, 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(os(i)) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call seteru (28hinits eta may be too small, 28, 1 1, 2) inits = i c return end function csevl (x, cs, n) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c evaluate the n-term chebyshev series cs at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox c and parker, chebyshev polys in numerical analysis, oxford press, p.56. c c input arguments -- c x value at which the series is to be evaluated. c cs array of n terms of a chebyshev series. in eval- c uating cs, only half the first coef is summed. c n number of terms in array cs. c dimension cs(1) c if (n.lt.1) call seteru (28hcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hcsevl number of terms gt 1000, 1 31, 3, 2) if (x.lt.(-1.1) .or. x.gt.1.1) call seteru ( 1 25hcsevl x outside (-1,+1), 25, 1, 1) c b1 = 0. b0 = 0. twox = 2.*x do 10 i=1,n b2 = b1 b1 = b0 ni = n + 1 - i b0 = twox*b1 - b2 + cs(ni) 10 continue c csevl = 0.5 * (b0-b2) c return end subroutine seteru (messg, nmessg, nerr, iopt) common /cseter/ iunflo integer messg(1) data iunflo / 0 / c if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt) if (iopt.ne.0) return c if (iunflo.le.0) return call seterr (messg, nmessg, nerr, 1) c return end function besi0e (x) c july 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bi0cs(12), ai0cs(21), ai02cs(22) c external csevl, exp, inits, r1mach, sqrt c c series for bi0 on the interval 0. to 9.00000d+00 c with weighted error 2.46e-18 c log weighted error 17.61 c significant figures required 17.90 c decimal places required 18.15 c data bi0 cs( 1) / -.0766054725 2839144951e0 / data bi0 cs( 2) / 1.9273379539 93808270e0 / data bi0 cs( 3) / .2282644586 920301339e0 / data bi0 cs( 4) / .0130489146 6707290428e0 / data bi0 cs( 5) / .0004344270 9008164874e0 / data bi0 cs( 6) / .0000094226 5768600193e0 / data bi0 cs( 7) / .0000001434 0062895106e0 / data bi0 cs( 8) / .0000000016 1384906966e0 / data bi0 cs( 9) / .0000000000 1396650044e0 / data bi0 cs(10) / .0000000000 0009579451e0 / data bi0 cs(11) / .0000000000 0000053339e0 / data bi0 cs(12) / .0000000000 0000000245e0 / c c series for ai0 on the interval 1.25000d-01 to 3.33333d-01 c with weighted error 7.87e-17 c log weighted error 16.10 c significant figures required 14.69 c decimal places required 16.76 c data ai0 cs( 1) / .0757599449 4023796e0 / data ai0 cs( 2) / .0075913808 1082334e0 / data ai0 cs( 3) / .0004153131 3389237e0 / data ai0 cs( 4) / .0000107007 6463439e0 / data ai0 cs( 5) / -.0000079011 7997921e0 / data ai0 cs( 6) / -.0000007826 1435014e0 / data ai0 cs( 7) / .0000002783 8499429e0 / data ai0 cs( 8) / .0000000082 5247260e0 / data ai0 cs( 9) / -.0000000120 4463945e0 / data ai0 cs(10) / .0000000015 5964859e0 / data ai0 cs(11) / .0000000002 2925563e0 / data ai0 cs(12) / -.0000000001 1916228e0 / data ai0 cs(13) / .0000000000 1757854e0 / data ai0 cs(14) / .0000000000 0112822e0 / data ai0 cs(15) / -.0000000000 0114684e0 / data ai0 cs(16) / .0000000000 0027155e0 / data ai0 cs(17) / -.0000000000 0002415e0 / data ai0 cs(18) / -.0000000000 0000608e0 / data ai0 cs(19) / .0000000000 0000314e0 / data ai0 cs(20) / -.0000000000 0000071e0 / data ai0 cs(21) / .0000000000 0000007e0 / c c series for ai02 on the interval 0. to 1.25000d-01 c with weighted error 3.79e-17 c log weighted error 16.42 c significant figures required 14.86 c decimal places required 17.09 c data ai02cs( 1) / .0544904110 1410882e0 / data ai02cs( 2) / .0033691164 7825569e0 / data ai02cs( 3) / .0000688975 8346918e0 / data ai02cs( 4) / .0000028913 7052082e0 / data ai02cs( 5) / .0000002048 9185893e0 / data ai02cs( 6) / .0000000226 6668991e0 / data ai02cs( 7) / .0000000033 9623203e0 / data ai02cs( 8) / .0000000004 9406022e0 / data ai02cs( 9) / .0000000000 1188914e0 / data ai02cs(10) / -.0000000000 3149915e0 / data ai02cs(11) / -.0000000000 1321580e0 / data ai02cs(12) / -.0000000000 0179419e0 / data ai02cs(13) / .0000000000 0071801e0 / data ai02cs(14) / .0000000000 0038529e0 / data ai02cs(15) / .0000000000 0001539e0 / data ai02cs(16) / -.0000000000 0004151e0 / data ai02cs(17) / -.0000000000 0000954e0 / data ai02cs(18) / .0000000000 0000382e0 / data ai02cs(19) / .0000000000 0000176e0 / data ai02cs(20) / -.0000000000 0000034e0 / data ai02cs(21) / -.0000000000 0000027e0 / data ai02cs(22) / .0000000000 0000003e0 / c data nti0, ntai0, ntai02, xsml / 3*0, 0. / c if (nti0.ne.0) go to 10 nti0 = inits (bi0cs, 12, 0.1*r1mach(3)) ntai0 = inits (ai0cs, 21, 0.1*r1mach(3)) ntai02 = inits (ai02cs, 22, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi0e = 1.0 if (y.gt.xsml) besi0e = exp(-y) * ( 2.75 + 1 csevl (y*y/4.5-1.0, bi0cs, nti0) ) return c 20 if (y.le.8.) besi0e = (.375 + csevl ((48./y-11.)/5., ai0cs, ntai0) 1 ) / sqrt(y) if (y.gt.8.) besi0e = (.375 + csevl (16./y-1., ai02cs, ntai02)) 1 / sqrt(y) c return end function besi1e (x) c oct 1983 version. w. fullerton, c3, los alamos scientific lab. dimension bi1cs(11), ai1cs(21), ai12cs(22) c external csevl, exp, inits, r1mach, sqrt c c series for bi1 on the interval 0. to 9.00000d+00 c with weighted error 2.40e-17 c log weighted error 16.62 c significant figures required 16.23 c decimal places required 17.14 c data bi1 cs( 1) / -.0019717132 61099859e0 / data bi1 cs( 2) / .4073488766 7546481e0 / data bi1 cs( 3) / .0348389942 99959456e0 / data bi1 cs( 4) / .0015453945 56300123e0 / data bi1 cs( 5) / .0000418885 21098377e0 / data bi1 cs( 6) / .0000007649 02676483e0 / data bi1 cs( 7) / .0000000100 42493924e0 / data bi1 cs( 8) / .0000000000 99322077e0 / data bi1 cs( 9) / .0000000000 00766380e0 / data bi1 cs(10) / .0000000000 00004741e0 / data bi1 cs(11) / .0000000000 00000024e0 / c c series for ai1 on the interval 1.25000d-01 to 3.33333d-01 c with weighted error 6.98e-17 c log weighted error 16.16 c significant figures required 14.53 c decimal places required 16.82 c data ai1 cs( 1) / -.0284674418 1881479e0 / data ai1 cs( 2) / -.0192295323 1443221e0 / data ai1 cs( 3) / -.0006115185 8579437e0 / data ai1 cs( 4) / -.0000206997 1253350e0 / data ai1 cs( 5) / .0000085856 1914581e0 / data ai1 cs( 6) / .0000010494 9824671e0 / data ai1 cs( 7) / -.0000002918 3389184e0 / data ai1 cs( 8) / -.0000000155 9378146e0 / data ai1 cs( 9) / .0000000131 8012367e0 / data ai1 cs(10) / -.0000000014 4842341e0 / data ai1 cs(11) / -.0000000002 9085122e0 / data ai1 cs(12) / .0000000001 2663889e0 / data ai1 cs(13) / -.0000000000 1664947e0 / data ai1 cs(14) / -.0000000000 0166665e0 / data ai1 cs(15) / .0000000000 0124260e0 / data ai1 cs(16) / -.0000000000 0027315e0 / data ai1 cs(17) / .0000000000 0002023e0 / data ai1 cs(18) / .0000000000 0000730e0 / data ai1 cs(19) / -.0000000000 0000333e0 / data ai1 cs(20) / .0000000000 0000071e0 / data ai1 cs(21) / -.0000000000 0000006e0 / c c series for ai12 on the interval 0. to 1.25000d-01 c with weighted error 3.55e-17 c log weighted error 16.45 c significant figures required 14.69 c decimal places required 17.12 c data ai12cs( 1) / .0285762350 1828014e0 / data ai12cs( 2) / -.0097610974 9136147e0 / data ai12cs( 3) / -.0001105889 3876263e0 / data ai12cs( 4) / -.0000038825 6480887e0 / data ai12cs( 5) / -.0000002512 2362377e0 / data ai12cs( 6) / -.0000000263 1468847e0 / data ai12cs( 7) / -.0000000038 3538039e0 / data ai12cs( 8) / -.0000000005 5897433e0 / data ai12cs( 9) / -.0000000000 1897495e0 / data ai12cs(10) / .0000000000 3252602e0 / data ai12cs(11) / .0000000000 1412580e0 / data ai12cs(12) / .0000000000 0203564e0 / data ai12cs(13) / -.0000000000 0071985e0 / data ai12cs(14) / -.0000000000 0040836e0 / data ai12cs(15) / -.0000000000 0002101e0 / data ai12cs(16) / .0000000000 0004273e0 / data ai12cs(17) / .0000000000 0001041e0 / data ai12cs(18) / -.0000000000 0000382e0 / data ai12cs(19) / -.0000000000 0000186e0 / data ai12cs(20) / .0000000000 0000033e0 / data ai12cs(21) / .0000000000 0000028e0 / data ai12cs(22) / -.0000000000 0000003e0 / c data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0. / c if (nti1.ne.0) go to 10 nti1 = inits (bi1cs, 11, 0.1*r1mach(3)) ntai1 = inits (ai1cs, 21, 0.1*r1mach(3)) ntai12 = inits (ai12cs, 22, 0.1*r1mach(3)) c xmin = 2.0*r1mach(1) xsml = sqrt (8.0*r1mach(3)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi1e = 0.0 if (y.eq.0.0) return c if (y.le.xmin) call seteru ( 1 37hbesi1e abs(x) so small i1 underflows, 37, 1, 0) besi1e = 0.0 if (y.gt.xmin) besi1e = 0.5*x if (y.gt.xsml) besi1e = x * (.875 + csevl(y*y/4.5-1., bi1cs,nti1)) besi1e = exp(-y) * besi1e return c 20 if (y.le.8.) besi1e = (.375 + csevl ((48./y-11.)/5., ai1cs, ntai1) 1 ) / sqrt(y) if (y.gt.8.) besi1e = (.375 + csevl (16./y-1.0, ai12cs, ntai12)) 1 / sqrt(y) besi1e = sign (besi1e, x) c return end function besk0e (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk0cs(11), ak0cs(17), ak02cs(14) c external alog, besi0, csevl, exp, inits, r1mach, sqrt c c series for bk0 on the interval 0. to 4.00000d+00 c with weighted error 3.57e-19 c log weighted error 18.45 c significant figures required 17.99 c decimal places required 18.97 c data bk0 cs( 1) / -.0353273932 3390276872e0 / data bk0 cs( 2) / .3442898999 246284869e0 / data bk0 cs( 3) / .0359799365 1536150163e0 / data bk0 cs( 4) / .0012646154 1144692592e0 / data bk0 cs( 5) / .0000228621 2103119451e0 / data bk0 cs( 6) / .0000002534 7910790261e0 / data bk0 cs( 7) / .0000000019 0451637722e0 / data bk0 cs( 8) / .0000000000 1034969525e0 / data bk0 cs( 9) / .0000000000 0004259816e0 / data bk0 cs(10) / .0000000000 0000013744e0 / data bk0 cs(11) / .0000000000 0000000035e0 / c c series for ak0 on the interval 1.25000d-01 to 5.00000d-01 c with weighted error 5.34e-17 c log weighted error 16.27 c significant figures required 14.92 c decimal places required 16.89 c data ak0 cs( 1) / -.0764394790 3327941e0 / data ak0 cs( 2) / -.0223565260 5699819e0 / data ak0 cs( 3) / .0007734181 1546938e0 / data ak0 cs( 4) / -.0000428100 6688886e0 / data ak0 cs( 5) / .0000030817 0017386e0 / data ak0 cs( 6) / -.0000002639 3672220e0 / data ak0 cs( 7) / .0000000256 3713036e0 / data ak0 cs( 8) / -.0000000027 4270554e0 / data ak0 cs( 9) / .0000000003 1694296e0 / data ak0 cs(10) / -.0000000000 3902353e0 / data ak0 cs(11) / .0000000000 0506804e0 / data ak0 cs(12) / -.0000000000 0068895e0 / data ak0 cs(13) / .0000000000 0009744e0 / data ak0 cs(14) / -.0000000000 0001427e0 / data ak0 cs(15) / .0000000000 0000215e0 / data ak0 cs(16) / -.0000000000 0000033e0 / data ak0 cs(17) / .0000000000 0000005e0 / c c series for ak02 on the interval 0. to 1.25000d-01 c with weighted error 2.34e-17 c log weighted error 16.63 c significant figures required 14.67 c decimal places required 17.20 c data ak02cs( 1) / -.0120186982 6307592e0 / data ak02cs( 2) / -.0091748526 9102569e0 / data ak02cs( 3) / .0001444550 9317750e0 / data ak02cs( 4) / -.0000040136 1417543e0 / data ak02cs( 5) / .0000001567 8318108e0 / data ak02cs( 6) / -.0000000077 7011043e0 / data ak02cs( 7) / .0000000004 6111825e0 / data ak02cs( 8) / -.0000000000 3158592e0 / data ak02cs( 9) / .0000000000 0243501e0 / data ak02cs(10) / -.0000000000 0020743e0 / data ak02cs(11) / .0000000000 0001925e0 / data ak02cs(12) / -.0000000000 0000192e0 / data ak02cs(13) / .0000000000 0000020e0 / data ak02cs(14) / -.0000000000 0000002e0 / c data ntk0, ntak0, ntak02, xsml / 3*0, 0. / c if (ntk0.ne.0) go to 10 ntk0 = inits (bk0cs, 11, 0.1*r1mach(3)) ntak0 = inits (ak0cs, 17, 0.1*r1mach(3)) ntak02 = inits (ak02cs, 14, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) c 10 if (x.le.0.) call seteru (29hbesk0e x is zero or negative, 29, 1 1, 2) if (x.gt.2.) go to 20 c y = 0. if (x.gt.xsml) y = x*x besk0e = exp(x) * (-alog(0.5*x)*besi0(x) 1 - .25 + csevl (.5*y-1., bk0cs, ntk0) ) return c 20 if (x.le.8.) besk0e = (1.25 + csevl ((16./x-5.)/3., ak0cs, ntak0)) 1 / sqrt(x) if (x.gt.8.) besk0e = (1.25 + csevl (16./x-1., ak02cs, ntak02)) 1 / sqrt(x) c return end function besk1e (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk1cs(11), ak1cs(17), ak12cs(14) c external alog, besi1, csevl, exp, inits, r1mach, sqrt c c series for bk1 on the interval 0. to 4.00000d+00 c with weighted error 7.02e-18 c log weighted error 17.15 c significant figures required 16.73 c decimal places required 17.67 c data bk1 cs( 1) / .0253002273 389477705e0 / data bk1 cs( 2) / -.3531559607 76544876e0 / data bk1 cs( 3) / -.1226111808 22657148e0 / data bk1 cs( 4) / -.0069757238 596398643e0 / data bk1 cs( 5) / -.0001730288 957513052e0 / data bk1 cs( 6) / -.0000024334 061415659e0 / data bk1 cs( 7) / -.0000000221 338763073e0 / data bk1 cs( 8) / -.0000000001 411488392e0 / data bk1 cs( 9) / -.0000000000 006666901e0 / data bk1 cs(10) / -.0000000000 000024274e0 / data bk1 cs(11) / -.0000000000 000000070e0 / c c series for ak1 on the interval 1.25000d-01 to 5.00000d-01 c with weighted error 6.06e-17 c log weighted error 16.22 c significant figures required 15.41 c decimal places required 16.83 c data ak1 cs( 1) / .2744313406 973883e0 / data ak1 cs( 2) / .0757198995 3199368e0 / data ak1 cs( 3) / -.0014410515 5647540e0 / data ak1 cs( 4) / .0000665011 6955125e0 / data ak1 cs( 5) / -.0000043699 8470952e0 / data ak1 cs( 6) / .0000003540 2774997e0 / data ak1 cs( 7) / -.0000000331 1163779e0 / data ak1 cs( 8) / .0000000034 4597758e0 / data ak1 cs( 9) / -.0000000003 8989323e0 / data ak1 cs(10) / .0000000000 4720819e0 / data ak1 cs(11) / -.0000000000 0604783e0 / data ak1 cs(12) / .0000000000 0081284e0 / data ak1 cs(13) / -.0000000000 0011386e0 / data ak1 cs(14) / .0000000000 0001654e0 / data ak1 cs(15) / -.0000000000 0000248e0 / data ak1 cs(16) / .0000000000 0000038e0 / data ak1 cs(17) / -.0000000000 0000006e0 / c c series for ak12 on the interval 0. to 1.25000d-01 c with weighted error 2.58e-17 c log weighted error 16.59 c significant figures required 15.22 c decimal places required 17.16 c data ak12cs( 1) / .0637930834 3739001e0 / data ak12cs( 2) / .0283288781 3049721e0 / data ak12cs( 3) / -.0002475370 6739052e0 / data ak12cs( 4) / .0000057719 7245160e0 / data ak12cs( 5) / -.0000002068 9392195e0 / data ak12cs( 6) / .0000000097 3998344e0 / data ak12cs( 7) / -.0000000005 5853361e0 / data ak12cs( 8) / .0000000000 3732996e0 / data ak12cs( 9) / -.0000000000 0282505e0 / data ak12cs(10) / .0000000000 0023720e0 / data ak12cs(11) / -.0000000000 0002176e0 / data ak12cs(12) / .0000000000 0000215e0 / data ak12cs(13) / -.0000000000 0000022e0 / data ak12cs(14) / .0000000000 0000002e0 / c data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0. / c if (ntk1.ne.0) go to 10 ntk1 = inits (bk1cs, 11, 0.1*r1mach(3)) ntak1 = inits (ak1cs, 17, 0.1*r1mach(3)) ntak12 = inits (ak12cs, 14, 0.1*r1mach(3)) c xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01) xsml = sqrt (4.0*r1mach(3)) c 10 if (x.le.0.) call seteru (29hbesk1e x is zero or negative, 29, 1 1, 2) if (x.gt.2.0) go to 20 c if (x.lt.xmin) call seteru (31hbesk1e x so small k1 overflows, 1 31, 2, 2) y = 0. if (x.gt.xsml) y = x*x besk1e = exp(x) * (alog(0.5*x)*besi1(x) + 1 (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x ) return c 20 if (x.le.8.) besk1e = (1.25 + csevl ((16./x-5.)/3., ak1cs, ntak1)) 1 / sqrt(x) if (x.gt.8.) besk1e = (1.25 + csevl (16./x-1., ak12cs, ntak12)) 1 / sqrt(x) c return end integer function i8save(isw,ivalue,set) c c if (isw = 1) i8save returns the current error number and c sets it to ivalue if set = .true. . c c if (isw = 2) i8save returns the current recovery switch and c sets it to ivalue if set = .true. . c logical set c integer iparam(2) c iparam(1) is the error number and iparam(2) is the recovery switch. c c start execution error free and with recovery turned off. c data iparam(1) /0/, iparam(2) /2/ c i8save=iparam(isw) if (set) iparam(isw)=ivalue c return c end subroutine eprint c c this subroutine prints the last error message, if any. c integer messg(1) c call e9rint(messg,1,1,.false.) return c end subroutine s88fmt( n, w, ifmt ) c c s88fmt replaces ifmt(1), ... , ifmt(n) with c the characters corresponding to the n least significant c digits of w. c integer n,w,ifmt(n) c integer nt,wt,digits(10) c data digits( 1) / 1h0 / data digits( 2) / 1h1 / data digits( 3) / 1h2 / data digits( 4) / 1h3 / data digits( 5) / 1h4 / data digits( 6) / 1h5 / data digits( 7) / 1h6 / data digits( 8) / 1h7 / data digits( 9) / 1h8 / data digits(10) / 1h9 / c nt = n wt = w c 10 if (nt .le. 0) return idigit = mod( wt, 10 ) ifmt(nt) = digits(idigit+1) wt = wt/10 nt = nt - 1 go to 10 c end subroutine seterr (messg, nmessg, nerr, iopt) c c this version modified by w. fullerton to dump if iopt = 1 and c not recovering. c seterr sets lerror = nerr, optionally prints the message and dumps c according to the following rules... c c if iopt = 1 and recovering - just remember the error. c if iopt = 1 and not recovering - print, dump and stop. c if iopt = 2 - print, dump and stop. c c input c c messg - the error message. c nmessg - the length of the message, in characters. c nerr - the error number. must have nerr non-zero. c iopt - the option. must have iopt=1 or 2. c c error states - c c 1 - message length not positive. c 2 - cannot have nerr=0. c 3 - an unrecovered error followed by another error. c 4 - bad value for iopt. c c only the first 72 characters of the message are printed. c c the error handler calls a subroutine named fdump to produce a c symbolic dump. to complete the package, a dummy version of fdump c is supplied, but it should be replaced by a locally written version c which at least gives a trace-back. c integer messg(1) external i1mach, i8save c c the unit for error messages. c iwunit=i1mach(4) c if (nmessg.ge.1) go to 10 c c a message of non-positive length is fatal. c write(iwunit,9000) 9000 format(52h1error 1 in seterr - message length not positive.) go to 60 c c nw is the number of words the message occupies. c 10 nw=(min0(nmessg,72)-1)/i1mach(6)+1 c if (nerr.ne.0) go to 20 c c cannot turn the error state off using seterr. c write(iwunit,9001) 9001 format(42h1error 2 in seterr - cannot have nerr=0// 1 34h the current error message follows///) call e9rint(messg,nw,nerr,.true.) itemp=i8save(1,1,.true.) go to 50 c c set lerror and test for a previous unrecovered error. c 20 if (i8save(1,nerr,.true.).eq.0) go to 30 c write(iwunit,9002) 9002 format(23h1error 3 in seterr -, 1 48h an unrecovered error followed by another error.// 2 48h the previous and current error messages follow.///) call eprint call e9rint(messg,nw,nerr,.true.) go to 50 c c save this message in case it is not recovered from properly. c 30 call e9rint(messg,nw,nerr,.true.) c if (iopt.eq.1 .or. iopt.eq.2) go to 40 c c must have iopt = 1 or 2. c write(iwunit,9003) 9003 format(42h1error 4 in seterr - bad value for iopt// 1 34h the current error message follows///) go to 50 c c test for recovery. c 40 if (iopt.eq.2) go to 50 c if (i8save(2,0,.false.).eq.1) return c c call eprint c stop c 50 call eprint 60 call fdump stop c end subroutine e9rint(messg,nw,nerr,save) c c this routine stores the current error message or prints the old one, c if any, depending on whether or not save = .true. . c integer messg(nw) logical save external i1mach, i8save c c messgp stores at least the first 72 characters of the previous c message. its length is machine dependent and must be at least c c 1 + 71/(the number of characters stored per integer word). c integer messgp(36),fmt(14),ccplus c c start with no previous message. c data messgp(1)/1h1/, nwp/0/, nerrp/0/ c c set up the format for printing the error message. c the format is simply (a1,14x,72axx) where xx=i1mach(6) is the c number of characters stored per integer word. c data ccplus / 1h+ / c data fmt( 1) / 1h( / data fmt( 2) / 1ha / data fmt( 3) / 1h1 / data fmt( 4) / 1h, / data fmt( 5) / 1h1 / data fmt( 6) / 1h4 / data fmt( 7) / 1hx / data fmt( 8) / 1h, / data fmt( 9) / 1h7 / data fmt(10) / 1h2 / data fmt(11) / 1ha / data fmt(12) / 1hx / data fmt(13) / 1hx / data fmt(14) / 1h) / c if (.not.save) go to 20 c c save the message. c nwp=nw nerrp=nerr do 10 i=1,nw 10 messgp(i)=messg(i) c go to 30 c 20 if (i8save(1,0,.false.).eq.0) go to 30 c c print the message. c iwunit=i1mach(4) write(iwunit,9000) nerrp 9000 format(7h error ,i4,4h in ) c call s88fmt(2,i1mach(6),fmt(12)) write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp) c 30 return c end subroutine fdump return end REAL FUNCTION SECOND () C C OBTAIN ELAPSED CPU TIME SINCE START OF RUN. C C ON CRAY -- THIS ROUTINE NOT NEEDED. C C ON CONVEX -- MORE ACCURATE TIMINGS ARE OBTAINED BY CALLING C THE VECLIB ROUTINE CPUTIME AS FOLLOWS C C SECOND = CPUTIME(0.0E0) C C ON UNIX -- (SUN, FOR EXAMPLE) THE FOLLOWING MAY BE USED C REAL TARRAY(2), ETIME, TOTAL TOTAL = ETIME(TARRAY) SECOND = TARRAY(1) + TARRAY(2) C RETURN END