SUBROUTINE VI0 (M, X, F, WORK, IWORK, INFO)
C***BEGIN PROLOGUE VI0
C***PURPOSE Computes the hyperbolic Bessel function of the first kind
C of order zero (I0) for a vector of real arguments
C***LIBRARY VFNLIB
C***CATEGORY C10B1
C***TYPE SINGLE PRECISION (VI0-S, DVI0-D)
C***KEYWORDS BESSEL FUNCTION,FIRST KIND,HYPERBOLIC BESSEL FUNCTION,
C MODIFIED BESSEL FUNCTION, ORDER ZERO, VECTORIZED
C***AUTHOR SAUNDERS, B. V., (NIST)
C COMPUTING AND APPLIED MATHEMATICS LABORATORY
C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C GAITHERSBURG, MD 20899
C BOISVERT, R. F., (NIST)
C COMPUTING AND APPLIED MATHEMATICS LABORATORY
C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C GAITHERSBURG, MD 20899
C***DESCRIPTION
C
C VI0 computes the modified (hyperbolic) Bessel function of the
C first kind of order zero (I0) for real arguments using uniform
C approximation by Chebyshev polynomials.
C
C
C P A R A M E T E R S
C
C M (Input) Integer (M .GT. 0)
C The number of arguments at which the function is to be
C evaluated.
C
C X (Input) Real array of length M
C The arguments at which the function is to be evaluated are
C stored in X(1) to X(M) in any order.
C
C F (Output) Real array of length M
C F(i) contains the value of the function at X(i), i=1,..,M.
C
C WORK (Work) Real vector of length 7*M
C
C IWORK (Work) Integer vector of length M
C
C INFO (Output) Integer
C Indicates status of computed result. The following table
C lists possible values and their meanings. If OK=Yes then
C all F(i) have been set, otherwise none have been set.
C
C INFO OK Description
C ------------------------------------------------------------
C 0 Yes Successfull execution.
C 1 No Error: M .LE. 0
C 2 No Error: Some abs(X(i)) so big I0 overflows.
C The index of the first offending argument is
C returned in IWORK(1).
C
C
C *********************************************************************
C This routine is a modification of the function BESI0 developed by
C W. Fullerton of LANL.
C *********************************************************************
C
C***REFERENCES Ronald F. Boisvert and Bonita V. Saunders, Portable
C Vectorized Software for Bessel Function Evaluation,
C ACM Transactions on Mathematical Software 18 (1992),
C pp 456-469.
C***ROUTINES CALLED WI0
C***REVISION HISTORY (YYMMDD)
C 910226 DATE WRITTEN
C 930203 Minor modifications to prologue.
C***END PROLOGUE VI0
C
C ----------
C PARAMETERS
C ----------
C
INTEGER INFO, IWORK, M
REAL F, X, WORK
C
DIMENSION X(M), F(M), WORK(7*M), IWORK(M)
C
C ---------------
C LOCAL VARIABLES
C ---------------
C
INTEGER IWB0, IWB1, IWB2, IWY, IWTC, IWYC, IWZC, JIN
C
C***FIRST EXECUTABLE STATEMENT VI0
C
C ... PARTITION WORK ARRAYS
C
IWY = 1
IWTC = IWY + M
IWYC = IWTC + M
IWZC = IWYC + M
IWB0 = IWZC + M
IWB1 = IWB0 + M
IWB2 = IWB1 + M
C Total = IWB2 + M
C
JIN = 1
C Total = JIN + M
C
C ... WI0 DOES ALL THE WORK
C
CALL WI0(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC),
+ IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO)
C
RETURN
END
SUBROUTINE WI0 (M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2,
+ INFO)
C***BEGIN PROLOGUE WI0
C***SUBSIDIARY
C***PURPOSE Computes the hyperbolic Bessel function of the first kind
C of order zero (I0) for a vector of arguments
C***LIBRARY VFNLIB
C***AUTHOR SAUNDERS, B. V., (NIST)
C COMPUTING AND APPLIED MATHEMATICS LABORATORY
C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C GAITHERSBURG, MD 20899
C BOISVERT, R. F., (NIST)
C COMPUTING AND APPLIED MATHEMATICS LABORATORY
C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C GAITHERSBURG, MD 20899
C***ROUTINES CALLED IWCS, R1MACH, WNGT, WGTHR, WGTLE, WGT, WLE,
C WSCTR, WCS
C***REVISION HISTORY (YYMMDD)
C 910226 DATE WRITTEN
C***END PROLOGUE WI0
C
C ----------
C PARAMETERS
C ----------
C
INTEGER INFO, INDX, M
REAL B0, B1, B2, F, X, Y, TCMP, YCMP, ZCMP
C
DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), Y(M),
+ TCMP(M), YCMP(M), ZCMP(M)
C
C ---------------
C LOCAL VARIABLES
C ---------------
C
INTEGER LAI0, LAI02, LBI0
PARAMETER ( LAI0=21, LAI02=22, LBI0=12 )
C
INTEGER I, IWCS, J, KEY, N, NTI0, NTAI0, NTAI02
REAL AI0CS, AI02CS, BI0CS, EPMACH, EPS, R1MACH, XSML,
+ XMAX
C
DIMENSION AI0CS(LAI0), AI02CS(LAI02), BI0CS(LBI0)
C
SAVE AI0CS, AI02CS, BI0CS, N, NTAI0, NTAI02, NTI0, XSML, XMAX
C
C----------------------------------------------------------------------
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----------------------------------------------------------------------
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----------------------------------------------------------------------
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
C----------------------------------------------------------------------
C
DATA NTI0 / 0 /
C
C***FIRST EXECUTABLE STATEMENT WI0
C
IF (M .LE. 0) GO TO 910
C
IF (NTI0.EQ.0) THEN
EPMACH = R1MACH(3)
EPS = 0.10E0*EPMACH
NTI0 = IWCS(BI0CS, LBI0, EPS)
NTAI0 = IWCS(AI0CS, LAI0, EPS)
NTAI02 = IWCS(AI02CS, LAI02, EPS)
XSML = SQRT(4.0E0*EPMACH)
XMAX = LOG(R1MACH(2))
ENDIF
C
DO 10 I=1,M
Y(I) = ABS(X(I))
10 CONTINUE
C
CALL WNGT(M,Y,XMAX,KEY)
IF (KEY .NE. 0) GO TO 920
C
C ----------------
C CASE Y .LE. XSML
C ----------------
C
DO 15 I=1,M
F(I) = 1.0E0
15 CONTINUE
C
C --------------------------
C CASE XSML .LT. Y .LE. 3.0
C --------------------------
C
CALL WGTLE(M,Y,XSML,3.0E0,N,INDX)
IF (N .GT. 0) THEN
CALL WGTHR(N,Y,INDX,YCMP)
DO 20 J=1,N
TCMP(J) = YCMP(J)**2/4.50E0 - 1.0E0
20 CONTINUE
CALL WCS(N,TCMP,BI0CS,NTI0,ZCMP,B0,B1,B2)
DO 30 J=1,N
ZCMP(J) = 2.750E0 + ZCMP(J)
30 CONTINUE
CALL WSCTR(N,ZCMP,INDX,F)
ENDIF
C
C -------------------------
C CASE 3.0 .LT. Y .LE. 8.0
C -------------------------
C
CALL WGTLE(M,Y,3.0E0,8.0E0,N,INDX)
IF (N .GT. 0) THEN
CALL WGTHR(N,Y,INDX,YCMP)
DO 50 J=1,N
TCMP(J) = (48.0E0/YCMP(J) - 11.0E0)/5.0E0
50 CONTINUE
CALL WCS(N,TCMP,AI0CS,NTAI0,ZCMP,B0,B1,B2)
DO 60 J=1,N
ZCMP(J) = EXP(YCMP(J))*(0.3750E0+ZCMP(J))/SQRT(YCMP(J))
60 CONTINUE
CALL WSCTR(N,ZCMP,INDX,F)
ENDIF
C
C ----------------
C CASE Y .GT. 8.0
C ----------------
C
CALL WGT(M,Y,8.0E0,N,INDX)
IF (N .GT. 0) THEN
CALL WGTHR(N,Y,INDX,YCMP)
DO 80 J=1,N
TCMP(J) = 16.0E0/YCMP(J) - 1.0E0
80 CONTINUE
CALL WCS(N,TCMP,AI02CS,NTAI02,ZCMP,B0,B1,B2)
DO 90 J=1,N
ZCMP(J) = EXP(YCMP(J))*(0.3750E0+ZCMP(J))/SQRT(YCMP(J))
90 CONTINUE
CALL WSCTR(N,ZCMP,INDX,F)
ENDIF
C
C -----
C EXITS
C -----
C
C ... NORMAL
C
INFO = 0
GO TO 999
C
C ... M .LE. 0
C
910 CONTINUE
INFO = 1
GO TO 999
C
C ... ABS(X) SO LARGE I0 OVERFLOWS
C
920 CONTINUE
INFO = 2
INDX(1) = KEY
GO TO 999
C
999 CONTINUE
RETURN
END
.