1      FUNCTION BVALU(T,A,N,K,IDERIV,X,INBV,WORK)
2C***BEGIN PROLOGUE  BVALU
3C***DATE WRITTEN   800901   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***CATEGORY NO.  E3,K6
6C***KEYWORDS  B-SPLINE,DATA FITTING,INTERPOLATION,SPLINE
7C***AUTHOR  AMOS, D. E., (SNLA)
8C***PURPOSE  Evaluates the B-representation of a B-spline at X for the
9C            function value or any of its derivatives.
10C***DESCRIPTION
11C
12C     Written by Carl de Boor and modified by D. E. Amos
13C
14C     Reference
15C         SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472.
16C
17C     Abstract
18C         BVALU is the BVALUE function of the reference.
19C
20C         BVALU evaluates the B-representation (T,A,N,K) of a B-spline
21C         at X for the function value on IDERIV = 0 or any of its
22C         derivatives on IDERIV = 1,2,...,K-1.  Right limiting values
23C         (right derivatives) are returned except at the right end
24C         point X=T(N+1) where left limiting values are computed.  The
25C         spline is defined on T(K) .LE. X .LE. T(N+1).  BVALU returns
26C         a fatal error message when X is outside of this interval.
27C
28C         To compute left derivatives or left limiting values at a
29C         knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1.
30C
31C         BVALU calls INTRV
32C
33C     Description of Arguments
34C         Input
35C          T       - knot vector of length N+K
36C          A       - B-spline coefficient vector of length N
37C          N       - number of B-spline coefficients
38C                    N = sum of knot multiplicities-K
39C          K       - order of the B-spline, K .GE. 1
40C          IDERIV  - order of the derivative, 0 .LE. IDERIV .LE. K-1
41C                    IDERIV=0 returns the B-spline value
42C          X       - argument, T(K) .LE. X .LE. T(N+1)
43C          INBV    - an initialization parameter which must be set
44C                    to 1 the first time BVALU is called.
45C
46C         Output
47C          INBV    - INBV contains information for efficient process-
48C                    ing after the initial call and INBV must not
49C                    be changed by the user.  Distinct splines require
50C                    distinct INBV parameters.
51C          WORK    - work vector of length 3*K.
52C          BVALU   - value of the IDERIV-th derivative at X
53C
54C     Error Conditions
55C         An improper input is a fatal error
56C***REFERENCES  C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*,
57C                 SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3,
58C                 JUNE 1977, PP. 441-472.
59C***ROUTINES CALLED  INTRV,XERROR
60C***END PROLOGUE  BVALU
61C
62C
63      INTEGER I,IDERIV,IDERP1,IHI,IHMKMJ,ILO,IMK,IMKPJ, INBV, IPJ,
64     1 IP1, IP1MJ, J, JJ, J1, J2, K, KMIDER, KMJ, KM1, KPK, MFLAG, N
65      REAL A, FKMJ, T, WORK, X
66C     DIMENSION T(N+K), WORK(3*K)
67      DIMENSION T(1), A(N), WORK(1)
68C***FIRST EXECUTABLE STATEMENT  BVALU
69      BVALU = 0.0E0
70      IF(K.LT.1) GO TO 102
71      IF(N.LT.K) GO TO 101
72      IF(IDERIV.LT.0 .OR. IDERIV.GE.K) GO TO 110
73      KMIDER = K - IDERIV
74C
75C *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1)
76C     (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)).
77      KM1 = K - 1
78      CALL INTRV(T, N+1, X, INBV, I, MFLAG)
79      IF (X.LT.T(K)) GO TO 120
80      IF (MFLAG.EQ.0) GO TO 20
81      IF (X.GT.T(I)) GO TO 130
82   10 IF (I.EQ.K) GO TO 140
83      I = I - 1
84      IF (X.EQ.T(I)) GO TO 10
85C
86C *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES
87C     WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K
88C
89   20 IMK = I - K
90      DO 30 J=1,K
91        IMKPJ = IMK + J
92        WORK(J) = A(IMKPJ)
93   30 CONTINUE
94      IF (IDERIV.EQ.0) GO TO 60
95      DO 50 J=1,IDERIV
96        KMJ = K - J
97        FKMJ = FLOAT(KMJ)
98        DO 40 JJ=1,KMJ
99          IHI = I + JJ
100          IHMKMJ = IHI - KMJ
101          WORK(JJ) = (WORK(JJ+1)-WORK(JJ))/(T(IHI)-T(IHMKMJ))*FKMJ
102   40   CONTINUE
103   50 CONTINUE
104C
105C *** COMPUTE VALUE AT *X* IN (T(I),(T(I+1)) OF IDERIV-TH DERIVATIVE,
106C     GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV).
107   60 IF (IDERIV.EQ.KM1) GO TO 100
108      IP1 = I + 1
109      KPK = K + K
110      J1 = K + 1
111      J2 = KPK + 1
112      DO 70 J=1,KMIDER
113        IPJ = I + J
114        WORK(J1) = T(IPJ) - X
115        IP1MJ = IP1 - J
116        WORK(J2) = X - T(IP1MJ)
117        J1 = J1 + 1
118        J2 = J2 + 1
119   70 CONTINUE
120      IDERP1 = IDERIV + 1
121      DO 90 J=IDERP1,KM1
122        KMJ = K - J
123        ILO = KMJ
124        DO 80 JJ=1,KMJ
125          WORK(JJ) = (WORK(JJ+1)*WORK(KPK+ILO)+WORK(JJ)
126     1              *WORK(K+JJ))/(WORK(KPK+ILO)+WORK(K+JJ))
127          ILO = ILO - 1
128   80   CONTINUE
129   90 CONTINUE
130  100 BVALU = WORK(1)
131      RETURN
132C
133C
134  101 CONTINUE
135      CALL XERROR( ' BVALU,  N DOES NOT SATISFY N.GE.K',34,2,1)
136      RETURN
137  102 CONTINUE
138      CALL XERROR( ' BVALU,  K DOES NOT SATISFY K.GE.1',34,2,1)
139      RETURN
140  110 CONTINUE
141      CALL XERROR( ' BVALU,  IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K',
142     1 49, 2, 1)
143      RETURN
144  120 CONTINUE
145      CALL XERROR( ' BVALU,  X IS N0T GREATER THAN OR EQUAL TO T(K)',
146     1 47, 2, 1)
147      RETURN
148  130 CONTINUE
149      CALL XERROR( ' BVALU,  X IS NOT LESS THAN OR EQUAL TO T(N+1)',
150     1 46, 2, 1)
151      RETURN
152  140 CONTINUE
153      CALL XERROR( ' BVALU,  A LEFT LIMITING VALUE CANN0T BE OBTAINED AT
154     1 T(K)',     57, 2, 1)
155      RETURN
156      END
157