1*DECK BSPVN
2      SUBROUTINE BSPVN (T, JHIGH, K, INDEX, X, ILEFT, VNIKX, WORK,
3     +   IWORK)
4C***BEGIN PROLOGUE  BSPVN
5C***PURPOSE  Calculate the value of all (possibly) nonzero basis
6C            functions at X.
7C***LIBRARY   SLATEC
8C***CATEGORY  E3, K6
9C***TYPE      SINGLE PRECISION (BSPVN-S, DBSPVN-D)
10C***KEYWORDS  EVALUATION OF B-SPLINE
11C***AUTHOR  Amos, D. E., (SNLA)
12C***DESCRIPTION
13C
14C     Written by Carl de Boor and modified by D. E. Amos
15C
16C     Abstract
17C         BSPVN is the BSPLVN routine of the reference.
18C
19C         BSPVN calculates the value of all (possibly) nonzero basis
20C         functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where
21C         T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine
22C         on the first call when INDEX=1.  ILEFT is such that T(ILEFT)
23C         .LE. X .LT. T(ILEFT+1).  A call to INTRV(T,N+1,X,ILO,ILEFT,
24C         MFLAG) produces the proper ILEFT.  BSPVN calculates using the
25C         basic algorithm needed in BSPVD.  If only basis functions are
26C         desired, setting JHIGH=K and INDEX=1 can be faster than
27C         calling BSPVD, but extra coding is required for derivatives
28C         (INDEX=2) and BSPVD is set up for this purpose.
29C
30C         Left limiting values are set up as described in BSPVD.
31C
32C     Description of Arguments
33C         Input
34C          T       - knot vector of length N+K, where
35C                    N = number of B-spline basis functions
36C                    N = sum of knot multiplicities-K
37C          JHIGH   - order of B-spline, 1 .LE. JHIGH .LE. K
38C          K       - highest possible order
39C          INDEX   - INDEX = 1 gives basis functions of order JHIGH
40C                          = 2 denotes previous entry with WORK, IWORK
41C                              values saved for subsequent calls to
42C                              BSPVN.
43C          X       - argument of basis functions,
44C                    T(K) .LE. X .LE. T(N+1)
45C          ILEFT   - largest integer such that
46C                    T(ILEFT) .LE. X .LT. T(ILEFT+1)
47C
48C         Output
49C          VNIKX   - vector of length K for spline values.
50C          WORK    - a work vector of length 2*K
51C          IWORK   - a work parameter.  Both WORK and IWORK contain
52C                    information necessary to continue for INDEX = 2.
53C                    When INDEX = 1 exclusively, these are scratch
54C                    variables and can be used for other purposes.
55C
56C     Error Conditions
57C         Improper input is a fatal error.
58C
59C***REFERENCES  Carl de Boor, Package for calculating with B-splines,
60C                 SIAM Journal on Numerical Analysis 14, 3 (June 1977),
61C                 pp. 441-472.
62C***ROUTINES CALLED  XERMSG
63C***REVISION HISTORY  (YYMMDD)
64C   800901  DATE WRITTEN
65C   890831  Modified array declarations.  (WRB)
66C   890831  REVISION DATE from Version 3.2
67C   891214  Prologue converted to Version 4.0 format.  (BAB)
68C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
69C   900326  Removed duplicate information from DESCRIPTION section.
70C           (WRB)
71C   920501  Reformatted the REFERENCES section.  (WRB)
72C***END PROLOGUE  BSPVN
73C
74      INTEGER ILEFT, IMJP1, INDEX, IPJ, IWORK, JHIGH, JP1, JP1ML, K, L
75      REAL T, VM, VMPREV, VNIKX, WORK, X
76C     DIMENSION T(ILEFT+JHIGH)
77      DIMENSION T(*), VNIKX(*), WORK(*)
78C     CONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS.
79C     WORK(I) = DELTAP(I), WORK(K+I) = DELTAM(I), I = 1,K
80C***FIRST EXECUTABLE STATEMENT  BSPVN
81      IF(K.LT.1) GO TO 90
82      IF(JHIGH.GT.K .OR. JHIGH.LT.1) GO TO 100
83      IF(INDEX.LT.1 .OR. INDEX.GT.2) GO TO 105
84      IF(X.LT.T(ILEFT) .OR. X.GT.T(ILEFT+1)) GO TO 110
85      GO TO (10, 20), INDEX
86   10 IWORK = 1
87      VNIKX(1) = 1.0E0
88      IF (IWORK.GE.JHIGH) GO TO 40
89C
90   20 IPJ = ILEFT + IWORK
91      WORK(IWORK) = T(IPJ) - X
92      IMJP1 = ILEFT - IWORK + 1
93      WORK(K+IWORK) = X - T(IMJP1)
94      VMPREV = 0.0E0
95      JP1 = IWORK + 1
96      DO 30 L=1,IWORK
97        JP1ML = JP1 - L
98        VM = VNIKX(L)/(WORK(L)+WORK(K+JP1ML))
99        VNIKX(L) = VM*WORK(L) + VMPREV
100        VMPREV = VM*WORK(K+JP1ML)
101   30 CONTINUE
102      VNIKX(JP1) = VMPREV
103      IWORK = JP1
104      IF (IWORK.LT.JHIGH) GO TO 20
105C
106   40 RETURN
107C
108C
109   90 CONTINUE
110      CALL XERMSG ('SLATEC', 'BSPVN', 'K DOES NOT SATISFY K.GE.1', 2,
111     +   1)
112      RETURN
113  100 CONTINUE
114      CALL XERMSG ('SLATEC', 'BSPVN',
115     +   'JHIGH DOES NOT SATISFY 1.LE.JHIGH.LE.K', 2, 1)
116      RETURN
117  105 CONTINUE
118      CALL XERMSG ('SLATEC', 'BSPVN', 'INDEX IS NOT 1 OR 2', 2, 1)
119      RETURN
120  110 CONTINUE
121      CALL XERMSG ('SLATEC', 'BSPVN',
122     +   'X DOES NOT SATISFY T(ILEFT).LE.X.LE.T(ILEFT+1)', 2, 1)
123      RETURN
124      END
125