1      SUBROUTINE DBSPVD(T,K,NDERIV,X,ILEFT,LDVNIK,VNIKX,WORK)
2C***BEGIN PROLOGUE  DBSPVD
3C***DATE WRITTEN   800901   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***REVISION HISTORY  (YYMMDD)
6C   000330  Modified array declarations.  (JEC)
7C
8C***CATEGORY NO.  E3,K6
9C***KEYWORDS  B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION,
10C             SPLINE
11C***AUTHOR  AMOS, D. E., (SNLA)
12C***PURPOSE  Calculates the value and all derivatives of order less than
13C            NDERIV of all basis functions which do not vanish at X.
14C***DESCRIPTION
15C
16C     Written by Carl de Boor and modified by D. E. Amos
17C
18C     Reference
19C         SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472.
20C
21C     Abstract    **** a double precision routine ****
22C
23C         DBSPVD is the BSPLVD routine of the reference.
24C
25C         DBSPVD calculates the value and all derivatives of order
26C         less than NDERIV of all basis functions which do not
27C         (possibly) vanish at X.  ILEFT is input such that
28C         T(ILEFT) .LE. X .LT. T(ILEFT+1).  A call to INTRV(T,N+1,X,
29C         ILO,ILEFT,MFLAG) will produce the proper ILEFT.  The output of
30C         DBSPVD is a matrix VNIKX(I,J) of dimension at least (K,NDERIV)
31C         whose columns contain the K nonzero basis functions and
32C         their NDERIV-1 right derivatives at X, I=1,K, J=1,NDERIV.
33C         These basis functions have indices ILEFT-K+I, I=1,K,
34C         K .LE. ILEFT .LE. N.  The nonzero part of the I-th basis
35C         function lies in (T(I),T(I+K)), I=1,N).
36C
37C         If X=T(ILEFT+1) then VNIKX contains left limiting values
38C         (left derivatives) at T(ILEFT+1).  In particular, ILEFT = N
39C         produces left limiting values at the right end point
40C         X=T(N+1).  To obtain left limiting values at T(I), I=K+1,N+1,
41C         set X= next lower distinct knot, call INTRV to get ILEFT,
42C         set X=T(I), and then call DBSPVD.
43C
44C         DBSPVD calls DBSPVN
45C
46C     Description of Arguments
47C         Input      T,X are double precision
48C          T       - knot vector of length N+K, where
49C                    N = number of B-spline basis functions
50C                    N = sum of knot multiplicities-K
51C          K       - order of the B-spline, K .GE. 1
52C          NDERIV  - number of derivatives = NDERIV-1,
53C                    1 .LE. NDERIV .LE. K
54C          X       - argument of basis functions,
55C                    T(K) .LE. X .LE. T(N+1)
56C          ILEFT   - largest integer such that
57C                    T(ILEFT) .LE. X .LT.  T(ILEFT+1)
58C          LDVNIK  - leading dimension of matrix VNIKX
59C
60C         Output     VNIKX,WORK are double precision
61C          VNIKX   - matrix of dimension at least (K,NDERIV) contain-
62C                    ing the nonzero basis functions at X and their
63C                    derivatives columnwise.
64C          WORK    - a work vector of length (K+1)*(K+2)/2
65C
66C     Error Conditions
67C         Improper input is a fatal error
68C***REFERENCES  C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*,
69C                 SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3,
70C                 JUNE 1977, PP. 441-472.
71C***ROUTINES CALLED  DBSPVN,XERROR
72C***END PROLOGUE  DBSPVD
73C
74C
75      INTEGER I,IDERIV,ILEFT,IPKMD,J,JJ,JLOW,JM,JP1MID,K,KMD, KP1, L,
76     1 LDUMMY, M, MHIGH, NDERIV
77      DOUBLE PRECISION FACTOR, FKMD, T, V, VNIKX, WORK, X
78C     DIMENSION T(ILEFT+K), WORK((K+1)*(K+2)/2)
79C     A(I,J) = WORK(I+J*(J+1)/2),  I=1,J+1  J=1,K-1
80C     A(I,K) = W0RK(I+K*(K-1)/2)  I=1.K
81C     WORK(1) AND WORK((K+1)*(K+2)/2) ARE NOT USED.
82      DIMENSION T(*), VNIKX(LDVNIK,NDERIV), WORK(*)
83C***FIRST EXECUTABLE STATEMENT  DBSPVD
84      IF(K.LT.1) GO TO 200
85      IF(NDERIV.LT.1 .OR. NDERIV.GT.K) GO TO 205
86      IF(LDVNIK.LT.K) GO TO 210
87      IDERIV = NDERIV
88      KP1 = K + 1
89      JJ = KP1 - IDERIV
90      CALL DBSPVN(T, JJ, K, 1, X, ILEFT, VNIKX, WORK, IWORK)
91      IF (IDERIV.EQ.1) GO TO 100
92      MHIGH = IDERIV
93      DO 20 M=2,MHIGH
94        JP1MID = 1
95        DO 10 J=IDERIV,K
96          VNIKX(J,IDERIV) = VNIKX(JP1MID,1)
97          JP1MID = JP1MID + 1
98   10   CONTINUE
99        IDERIV = IDERIV - 1
100        JJ = KP1 - IDERIV
101        CALL DBSPVN(T, JJ, K, 2, X, ILEFT, VNIKX, WORK, IWORK)
102   20 CONTINUE
103C
104      JM = KP1*(KP1+1)/2
105      DO 30 L = 1,JM
106        WORK(L) = 0.0D0
107   30 CONTINUE
108C     A(I,I) = WORK(I*(I+3)/2) = 1.0       I = 1,K
109      L = 2
110      J = 0
111      DO 40 I = 1,K
112        J = J + L
113        WORK(J) = 1.0D0
114        L = L + 1
115   40 CONTINUE
116      KMD = K
117      DO 90 M=2,MHIGH
118        KMD = KMD - 1
119        FKMD = FLOAT(KMD)
120        I = ILEFT
121        J = K
122        JJ = J*(J+1)/2
123        JM = JJ - J
124        DO 60 LDUMMY=1,KMD
125          IPKMD = I + KMD
126          FACTOR = FKMD/(T(IPKMD)-T(I))
127          DO 50 L=1,J
128            WORK(L+JJ) = (WORK(L+JJ)-WORK(L+JM))*FACTOR
129   50     CONTINUE
130          I = I - 1
131          J = J - 1
132          JJ = JM
133          JM = JM - J
134   60   CONTINUE
135C
136        DO 80 I=1,K
137          V = 0.0D0
138          JLOW = MAX0(I,M)
139          JJ = JLOW*(JLOW+1)/2
140          DO 70 J=JLOW,K
141            V = WORK(I+JJ)*VNIKX(J,M) + V
142            JJ = JJ + J + 1
143   70     CONTINUE
144          VNIKX(I,M) = V
145   80   CONTINUE
146   90 CONTINUE
147  100 RETURN
148C
149C
150  200 CONTINUE
151      CALL XERROR( ' DBSPVD,  K DOES NOT SATISFY K.GE.1', 35, 2, 1)
152      RETURN
153  205 CONTINUE
154      CALL XERROR( ' DBSPVD,  NDERIV DOES NOT SATISFY 1.LE.NDERIV.LE.K',
155     1 50, 2, 1)
156      RETURN
157  210 CONTINUE
158      CALL XERROR( ' DBSPVD,  LDVNIK DOES NOT SATISFY LDVNIK.GE.K',45,
159     1 2, 1)
160      RETURN
161      END
162