      SUBROUTINE VK1 (M, X, F, WORK, IWORK, INFO)
C***BEGIN PROLOGUE  VK1
C***PURPOSE  Computes the hyperbolic Bessel function of the third kind
C            of order one (K1) for a vector of real arguments
C***LIBRARY   VFNLIB
C***CATEGORY  C10B1
C***TYPE      SINGLE PRECISION (VK1-S, DVK1-D)
C***KEYWORDS  BESSEL FUNCTION,THIRD 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   VK1 computes the modified (hyperbolic) Bessel function of the
C   third kind of order one (K1) 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 X(i) so big K1 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 X(i) is zero or negative.
C                     The index of the first offending argument is
C                     returned in IWORK(1).
C            3   No   Error: Some X(i) so small K1 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 BESK1 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  WK1
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C   930203  Minor modifications to prologue.
C***END PROLOGUE  VK1
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, IWTC, IWXC, IWYC, IWZC, JIN
C
C***FIRST EXECUTABLE STATEMENT  VK1
C
C     ... PARTITION WORK ARRAYS
C
      IWTC   = 1
      IWXC  = IWTC + M
      IWYC  = IWXC + M
      IWZC  = IWYC + M
      IWB0  = IWZC + M
      IWB1  = IWB0 + M
      IWB2  = IWB1 + M
C     Total = IB2  + M
C
      JIN   = 1
C     Total = JIN  + M
C
C     ... WK1 DOES ALL THE WORK
C
      CALL WK1(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC),
     +         IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO)
C
      RETURN
      END

      SUBROUTINE WK1 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1,
     +   B2, INFO)
C***BEGIN PROLOGUE  WK1
C***SUBSIDIARY
C***PURPOSE  Computes the hyperbolic Bessel function of the third kind
C            of order one (K1) 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, WNLE, WGTHR, WGTLE, WGT, WLE,
C                    WSCTR, WCS
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C***END PROLOGUE  WK1
C
C  ----------
C  PARAMETERS
C  ----------
C
      INTEGER INFO, INDX, M
      REAL    B0, B1, B2, F, X, TCMP, XCMP, YCMP, ZCMP
C
      DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), TCMP(M),
     +          XCMP(M), YCMP(M), ZCMP(M)
C
C  ---------------
C  LOCAL VARIABLES
C  ---------------
C
      INTEGER LAK1, LAK12, LBI1, LBK1
      PARAMETER ( LAK1=17, LAK12=14, LBI1=11, LBK1=11 )
C
      INTEGER I, IWCS, J, KEY, N, NTI1, NTK1, NTAK1, NTAK12, NTOT
      REAL    AK1CS, AK12CS, BI1CS, BK1CS, EPMACH, EPS, R1MACH,
     +        XMAX, XMIN, XSML
C
      DIMENSION AK1CS(LAK1), AK12CS(LAK12), BI1CS(LBI1), BK1CS(LBK1)
C
      SAVE AK1CS, AK12CS, BI1CS, BK1CS, NTAK1, NTAK12, NTI1, NTK1,
     +     XMAX, XMIN, XSML
C
C----------------------------------------------------------------------
C
C Series for BK1        on the interval  0.          to  4.00000D+00
C                                        with weighted error   7.02E-18
C                                         log weighted error  17.15
C                               significant figures required  16.73
C                                    decimal places required  17.67
C
      DATA BK1 CS( 1) /    .0253002273 389477705E0 /
      DATA BK1 CS( 2) /   -.3531559607 76544876E0 /
      DATA BK1 CS( 3) /   -.1226111808 22657148E0 /
      DATA BK1 CS( 4) /   -.0069757238 596398643E0 /
      DATA BK1 CS( 5) /   -.0001730288 957513052E0 /
      DATA BK1 CS( 6) /   -.0000024334 061415659E0 /
      DATA BK1 CS( 7) /   -.0000000221 338763073E0 /
      DATA BK1 CS( 8) /   -.0000000001 411488392E0 /
      DATA BK1 CS( 9) /   -.0000000000 006666901E0 /
      DATA BK1 CS(10) /   -.0000000000 000024274E0 /
      DATA BK1 CS(11) /   -.0000000000 000000070E0 /
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 AK1        on the interval  1.25000D-01 to  5.00000D-01
C                                        with weighted error   6.06E-17
C                                         log weighted error  16.22
C                               significant figures required  15.41
C                                    decimal places required  16.83
C
      DATA AK1 CS( 1) /    .2744313406 973883E0 /
      DATA AK1 CS( 2) /    .0757198995 3199368E0 /
      DATA AK1 CS( 3) /   -.0014410515 5647540E0 /
      DATA AK1 CS( 4) /    .0000665011 6955125E0 /
      DATA AK1 CS( 5) /   -.0000043699 8470952E0 /
      DATA AK1 CS( 6) /    .0000003540 2774997E0 /
      DATA AK1 CS( 7) /   -.0000000331 1163779E0 /
      DATA AK1 CS( 8) /    .0000000034 4597758E0 /
      DATA AK1 CS( 9) /   -.0000000003 8989323E0 /
      DATA AK1 CS(10) /    .0000000000 4720819E0 /
      DATA AK1 CS(11) /   -.0000000000 0604783E0 /
      DATA AK1 CS(12) /    .0000000000 0081284E0 /
      DATA AK1 CS(13) /   -.0000000000 0011386E0 /
      DATA AK1 CS(14) /    .0000000000 0001654E0 /
      DATA AK1 CS(15) /   -.0000000000 0000248E0 /
      DATA AK1 CS(16) /    .0000000000 0000038E0 /
      DATA AK1 CS(17) /   -.0000000000 0000006E0 /
C
C----------------------------------------------------------------------
C
C Series for AK12       on the interval  0.          to  1.25000D-01
C                                        with weighted error   2.58E-17
C                                         log weighted error  16.59
C                               significant figures required  15.22
C                                    decimal places required  17.16
C
      DATA AK12CS( 1) /    .0637930834 3739001E0 /
      DATA AK12CS( 2) /    .0283288781 3049721E0 /
      DATA AK12CS( 3) /   -.0002475370 6739052E0 /
      DATA AK12CS( 4) /    .0000057719 7245160E0 /
      DATA AK12CS( 5) /   -.0000002068 9392195E0 /
      DATA AK12CS( 6) /    .0000000097 3998344E0 /
      DATA AK12CS( 7) /   -.0000000005 5853361E0 /
      DATA AK12CS( 8) /    .0000000000 3732996E0 /
      DATA AK12CS( 9) /   -.0000000000 0282505E0 /
      DATA AK12CS(10) /    .0000000000 0023720E0 /
      DATA AK12CS(11) /   -.0000000000 0002176E0 /
      DATA AK12CS(12) /    .0000000000 0000215E0 /
      DATA AK12CS(13) /   -.0000000000 0000022E0 /
      DATA AK12CS(14) /    .0000000000 0000002E0 /
C
C----------------------------------------------------------------------
C
      DATA NTK1 / 0 /
C
C***FIRST EXECUTABLE STATEMENT  WK1
C
      IF (M .LE. 0)  GO TO 910
C
      IF (NTK1 .EQ. 0) THEN
         EPMACH = R1MACH(3)
         EPS = 0.10E0*EPMACH
         NTK1 = IWCS(BK1CS, LBK1, EPS)
         NTI1 = IWCS(BI1CS, LBI1, EPS)
         NTAK1 = IWCS(AK1CS, LAK1, EPS)
         NTAK12 = IWCS(AK12CS, LAK12, EPS)
         XMIN = EXP(MAX( LOG(R1MACH(1)), -LOG(R1MACH(2))) + 0.010E0)
         XSML = SQRT(4.0E0*EPMACH)
         XMAX = -LOG(R1MACH(1))
         XMAX = XMAX - 0.50E0*XMAX*LOG(XMAX)/(XMAX + 0.50E0)
      ENDIF
C
      NTOT = 0
C
      CALL WNLE(M,X,0.0E0,KEY)
      IF (KEY .NE. 0)  GO TO 920
C
      CALL WNLE(M,X,XMIN,KEY)
      IF (KEY .NE. 0)  GO TO 930
C
C  ------------------
C  CASE   X .GT. XMAX
C  ------------------
C
C     NOTE -- K0 UNDERFLOWS FOR X .GT. XMAX
C
      DO 5 I=1,M
         F(I) = 0.0E0
   5  CONTINUE
C
C  ----------------
C  CASE  X .LE. 2.0
C  ----------------
C
      CALL WLE(M,X,2.0E0,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,X,INDX,XCMP)
C
C        ... COMPUTE I1(X) ... RESULT IN ZCMP
C
         DO 20 J=1,N
            TCMP(J) = XCMP(J)**2/4.50E0 - 1.0E0
   20    CONTINUE
         CALL WCS(N,TCMP,BI1CS,NTI1,ZCMP,B0,B1,B2)
         DO 30 J=1,N
            ZCMP(J) = XCMP(J)*(0.8750E0 + ZCMP(J))
   30    CONTINUE
C
         DO 40 J=1,N
            TCMP(J) = 0.50E0*XCMP(J)**2 - 1.0E0
   40    CONTINUE
         CALL WCS(N,TCMP,BK1CS,NTK1,YCMP,B0,B1,B2)
         DO 50 J=1,N
            ZCMP(J) = LOG(0.50E0*XCMP(J))*ZCMP(J) +
     +                (0.750E0 + YCMP(J))/XCMP(J)
   50    CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
C
C  -------------------------
C  CASE  2.0 .LT. X .LE. 8.0
C  -------------------------
C
      CALL WGTLE(M,X,2.0E0,8.0E0,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,X,INDX,XCMP)
         DO 60 J=1,N
            TCMP(J) = (16.0E0/XCMP(J) - 5.0E0)/3.0E0
   60    CONTINUE
         CALL WCS(N,TCMP,AK1CS,NTAK1,YCMP,B0,B1,B2)
         DO 70 J=1,N
            ZCMP(J) = EXP(-XCMP(J))*(1.250E0 + YCMP(J))/SQRT(XCMP(J))
   70    CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
C
C  --------------------------
C  CASE  8.0 .LT. X .LE. XMAX
C  --------------------------
C
      CALL WGTLE(M,X,8.0E0,XMAX,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,X,INDX,XCMP)
         DO 80 J=1,N
            TCMP(J) = 16.0E0/XCMP(J) - 1.0E0
   80    CONTINUE
         CALL WCS(N,TCMP,AK12CS,NTAK12,YCMP,B0,B1,B2)
         DO 90 J=1,N
            ZCMP(J) = EXP(-XCMP(J))*(1.250E0 + YCMP(J))/SQRT(XCMP(J))
   90    CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
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     ... X IS ZERO OR NEGATIVE
C
  920 CONTINUE
      INFO = 2
      INDX(1) = KEY
      GO TO 999
C
C     ... X SO SMALL K1 OVERFLOWS
C
  930 CONTINUE
      INFO = 3
      INDX(1) = KEY
      GO TO 999
C
  999 CONTINUE
      RETURN
      END