Portable Vectorized Software for Bessel Function Evaluation Ronald F. Boisvert and Bonita V. Saunders Computing and Applied Mathematics Laboratory National Institute of Standards and Technology Gaithersburg, MD 20899 boisvert@cam.nist.gov saunders@cam.nist.gov June 19, 1991 Revised February 3, 1993 VFNLIB is a suite of Fortran subprograms for the evaluation of Bessel functions and modified Bessel functions of orders zero and one for a vector of real arguments. Distinguishing characteristics of these subprograms are that (a) they are portable across a wide range of machines, and (b) they are vectorized in the case when multiple function evaluation are to be performed. Single precision subprograms : VI0 : hyperbolic Bessel function of the first kind of order zero (I0) for a vector of real arguments VI1 : hyperbolic Bessel function of the first kind of order one (I1) for a vector of real arguments VJ0 : Bessel function of the first kind of order zero (J0) for a vector of real arguments VJ1 : Bessel function of the first kind of order one (J1) for a vector of real arguments VK0 : hyperbolic Bessel function of the third kind of order zero (K0) for a vector of real arguments VK1 : hyperbolic Bessel function of the third kind of order one (K1) for a vector of real arguments VY0 : Bessel function of the second kind of order zero (Y0) for a vector of real arguments VY1 : Bessel function of the second kind of order one (Y1) for a vector of real arguments Double precision subprograms : DVI0 : hyperbolic Bessel function of the first kind of order zero (I0) for a vector of real arguments DVI1 : hyperbolic Bessel function of the first kind of order one (I1) for a vector of real arguments DVJ0 : Bessel function of the first kind of order zero (J0) for a vector of real arguments DVJ1 : Bessel function of the first kind of order one (J1) for a vector of real arguments DVK0 : hyperbolic Bessel function of the third kind of order zero (K0) for a vector of real arguments DVK1 : hyperbolic Bessel function of the third kind of order one (K1) for a vector of real arguments DVY0 : Bessel function of the second kind of order zero (Y0) for a vector of real arguments DVY1 : Bessel function of the second kind of order one (Y1) for a vector of real arguments Calling sequences : All VFNLIB routines have calling sequences of the form CALL NAME (M, X, F, WORK, IWORK, INFO) where M is the number of arguments, X is the vector of arguments, F is the vector of function values, WORK is a workspace vector of length 7M, IWORK is a workspace vector of length M, and INFO is a return code. Test programs : Single and double precision test programs are included which compare the accuracy and timing of each VFNLIB Bessel function with its FNLIB counterpart. Portability notes : These routine are (mostly) coded in ANSI Fortran 77. The following should be noted. (a) The PORT library machine constants routines I1MACH, R1MACH, and D1MACH are used, but are not included. Use the netlib request "send i1mach r1mach d1mach from port" to obtain them. These routines need to be modified when moving VFNLIB to a new machine type. (b) The accuracy of these routines is limited by the precision with which Chebyshev series coefficients are represented. The coefficients in VFNLIB are given to about 16 decimal digits in single precision and 31 digits in double precision. This is sufficient for 64-bit/128-bit arithmetic as found on Cray Research computers, for example. A problem occurs when compiler automatic precision doubling switches are used on machines with 32-bit/64-bit arithmetic (IEEE arithmetic is one such case). When auto-doubling is used the effective precision can be increased to about 16 decimal digits in single precision and 34 decimal digits in double precision. In this case VFNLIB routines will return a fatal error because the precision requested in evaluating the Chebyshev series will be too great. If accuracy comparable to the Cray is sufficient in these cases, then simply comment out the IF statement immediately following statement number 20 in the routines IWCS and IDWCS. (c) The test routines call the real function SECOND to obtain the elapsed user CP time from the start of the run. This routine is distributed with calls to Unix time utilities. This may need changing on other systems. See the comments in the code for suggestions.