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.