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.