       SUBROUTINE VI1 (M, X, F, WORK, IWORK, INFO)
C***BEGIN PROLOGUE  VI1
C***PURPOSE  Computes the hyperbolic Bessel function of the first kind
C            of order one (I1) for a vector of real arguments
C***LIBRARY   VFNLIB
C***CATEGORY  C10B1
C***TYPE      SINGLE PRECISION (VI1-S, DVI1-D)
C***KEYWORDS  BESSEL FUNCTION,FIRST KIND,HYPERBOLIC BESSEL FUNCTION,
C             MODIFIED BESSEL FUNCTION, ORDER ONE, 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   VI1 computes the modified (hyperbolic) Bessel function of the
C   first kind of order one (I1) 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
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           -1   Yes  Warning: Some abs(X(i)) so small I1 underflows.
C                     The corresponding F(i) are set to zero.
C            0   Yes  Successfull execution.
C            1   No   Error: M .LE. 0
C            2   No   Error: Some abs(X(i)) so big I1 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 BESI1 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  WI1
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C   930203  Minor modifications to prologue.
C***END PROLOGUE  VI1
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  VI1
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     ... WI1 DOES ALL THE WORK
C
      CALL WI1(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC),
     +         IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO)
C
      RETURN
      END

      SUBROUTINE WI1(M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2,
     +   INFO)
C***BEGIN PROLOGUE  WI1
C***SUBSIDIARY
C***PURPOSE  Computes the hyperbolic Bessel function of the first kind
C            of order one (I1) 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  WI1
C
C  ----------
C  PARAMETERS
C  ----------
C
      INTEGER LAI1, LAI12, LBI1
      PARAMETER ( LAI1=21, LAI12=22, LBI1=11 )
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 I, IWCS, J, KEY, N, NTI1, NTAI1, NTAI12, NTOT
      REAL    AI1CS, AI12CS, BI1CS, EPMACH, EPS, R1MACH,
     +        XMAX, XMIN, XSML
C
      DIMENSION AI1CS(LAI1), AI12CS(LAI12), BI1CS(LBI1)
C
      SAVE AI1CS, AI12CS, BI1CS, N, NTAI1, NTAI12, NTI1,
     +     XMAX, XMIN, XSML
C
C----------------------------------------------------------------------
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----------------------------------------------------------------------
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----------------------------------------------------------------------
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
C----------------------------------------------------------------------
C
      DATA NTI1 / 0 /
C
C***FIRST EXECUTABLE STATEMENT  WI1
C
      IF (M .LE. 0)  GO TO 910
C
      IF (NTI1.EQ.0) THEN
         EPMACH = R1MACH(3)
         EPS = 0.10E0*EPMACH
         NTI1 = IWCS (BI1CS, LBI1, EPS)
         NTAI1 = IWCS (AI1CS, LAI1, EPS)
         NTAI12 = IWCS (AI12CS, LAI12, EPS)
         XMIN = 2.0E0*R1MACH(1)
         XSML = SQRT(8.0E0*EPMACH)
         XMAX = LOG(R1MACH(2))
      ENDIF
C
      NTOT = 0
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=0  OR  Y TOO SMALL
C  ---------------------------
C
C     NOTE -- I0 UNDERFLOWS FOR X .LE. XMIN
C
      DO 15 I=1,M
         F(I) = 0.0E0
  15  CONTINUE
C
C  ----------------------------
C  CASE   XMIN .LT. Y .LE. XSML
C  ----------------------------
C
      CALL WGTLE(M,Y,XMIN,XSML,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 20 J=1,N
            ZCMP(J) = 0.50E0*YCMP(J)
  20     CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
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
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 30 J=1,N
            TCMP(J) = YCMP(J)**2/4.50E0 - 1.0E0
  30     CONTINUE
         CALL WCS(N,TCMP,BI1CS,NTI1,ZCMP,B0,B1,B2)
         DO 40 J=1,N
            ZCMP(J) = YCMP(J)*(0.8750E0 + ZCMP(J))
  40     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
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 70 J=1,N
            TCMP(J) = (48.0E0/YCMP(J) - 11.0E0)/5.0E0
  70     CONTINUE
         CALL WCS(N,TCMP,AI1CS,NTAI1,ZCMP,B0,B1,B2)
         DO 80 J=1,N
            ZCMP(J) = EXP(YCMP(J))*(0.3750E0 + ZCMP(J))/SQRT(YCMP(J))
  80     CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
C
C  ----------------
C  CASE  8.0 .LT. Y
C  ----------------
C
      CALL WGT(M,Y,8.0E0,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 100 J=1,N
            TCMP(J) = 16.0E0/YCMP(J) - 1.0E0
  100     CONTINUE
         CALL WCS(N,TCMP,AI12CS,NTAI12,ZCMP,B0,B1,B2)
         DO 110 J=1,N
            ZCMP(J) = EXP(YCMP(J))*(0.3750E0 + ZCMP(J))/SQRT(YCMP(J))
  110    CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
C
C  -------------------------------------------
C  REVERSE SIGN FOR NEGATIVE X (RESULT IS ODD)
C  -------------------------------------------
C
      DO 200 I=1,M
         IF (X(I) .LT. 0.0E0)  F(I) = -F(I)
  200 CONTINUE
C
C  -----
C  EXITS
C  -----
C
C     ... NORMAL
C
      INFO = 0
      IF (NTOT .NE. M)  INFO = -1
      GO TO 999
C
C     ... M .LE. 0
C
  910 CONTINUE
      INFO = 1
      GO TO 999
C
C     ... ABS(X) SO LARGE I1 OVERFLOWS
C
  920 CONTINUE
      INFO = 2
      INDX(1) = KEY
      GO TO 999
C
  999 CONTINUE
      RETURN
      END