1      subroutine DSBASD(KORDER, LEFT, T, X, IDERIV, BDERIV)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 1994-10-20 DSBASD Krogh  Changes to use M77CON
6c>> 1993-01-12 DSBASD CLL  Bringing all computation inline.
7c>> 1992-11-02 C. L. Lawson, JPL
8c>> 1988-03-22 C. L. Lawson, JPL
9c
10c     Given the knot-set T(i), i = 1, ..., NT, this subr computes
11c     values at X of the derivative of order IDERIV of each of the
12c     the KORDER B-spline basis functions of order KORDER
13c     that are nonzero on the open interval (T(LEFT), T(LEFT+1).
14c     Require T(LEFT) < T(LEFT+1) and KORDER .le. LEFT .le.
15c     NT-KORDER.
16c     For proper evaluation, X must lie in the closed interval
17c     [T(LEFT), T(LEFT+1], otherwise the computation constitutes
18c     extrapolation.
19c     The basis functions whose derivatives are returned will be those
20c     indexed from LEFT+1-KORDER through LEFT, and their values will be
21c     stored in BDERIV(I), I = 1, ..., KORDER.
22c     ------------------------------------------------------------------
23c  Method:
24c
25c  In general there are two stages: First compute values of basis
26c  functions of order KORDER-IDERIV.  Then, if IDERIV > 0, transform
27c  these values to produce values of higher order derivatives of
28c  higher order basis functions, ultimately obtaining values of
29c  the IDERIV derivative of basis functions of order KORDER.
30c
31c  The first stage uses the recursion:
32C
33C                       X - T(I)              T(I+J+1) - X
34C     B(I,J+1)(X)  =  -----------B(I,J)(X) + ---------------B(I+1,J)(X)
35C                     T(I+J)-T(I)            T(I+J+1)-T(I+1)
36C
37c  where B(I,J) denotes the basis function of order J and index I.
38c  For order J, the only basis functions that can be nonzero on
39c  [T(LEFT), T(LEFT+1] are those indexed from LEFT+1-J to LEFT.
40c  For order 1, the only basis function nonzero on this interval is
41c  is B(LEFT,1), whose value is 1.  The organization of the calculation
42c  using the above recursion follows Algorithm (8) in Chapter X of the
43c  book by DeBoor.
44c
45c  For the second stage, let B(ID, K, J) denote the value at X of the
46c  IDth derivative of the basis function of order K, with index J.
47c  From the first stage we will have values of B(0, KORDER-IDERIV, j)
48c  for j = LEFT+1-KORDER+IDERIV), ..., LEFT, stored in BDERIV(i), for
49c  i = 1+IDERIV, ..., KORDER.
50c
51c     Loop for ID = 1, ..., IDERIV
52c        Replace the contents of BDERIV() by values of
53c        B(ID, KORDER-IDERIV+ID, j) for
54c        j = LEFT+1-KORDER+IDERIV-ID), ..., LEFT, storing these in
55c        BDERIV(i), i = 1+IDERIV-ID, ..., KORDER.
56c     End loop
57c
58c  The above loop uses formula (10), p. 138, of
59c  A PRACTICAL GUIDE TO SPLINES by Carl DeBoor, Springer-Verlag,
60c  1978, when ID = 1, and successive derivatives of that formula
61c  when ID > 1.  Note that we are using Formula (10) in the special
62c  case in which (Alpha sub i) = 1, and all other Alpha's are zero.
63c
64c  This approach and implementation by C. L. Lawson, JPL, March 1988.
65C     ------------------------------------------------------------------
66C  KORDER [in]  Gives both the order of the B-spline basis functions
67c        whose derivatives are to be evaluated,
68c        and the number of such functions to be evaluated.
69c        Note that the polynomial degree of these functions is one less
70c        than the order.  Example:  For cubic splines the order is 4.
71c        Require 1 .le. KORDER .le. KMAX, where KMAX is an internal
72c        parameter.
73c  LEFT  [in]  Index identifying the left end of the interval
74c        [T(LEFT), T(LEFT+1)] relative to which the B-spline basis
75c        functions will be evaluated.
76c        Require        KORDER .le. LEFT .le. NT-KORDER
77c        and            T(LEFT) < T(LEFT+1)
78c        The evaluation is proper if T(LEFT) .le. X .le. T(LEFT+1), and
79c        otherwise constitutes extrapolation.
80C        DIVISION BY ZERO  will result if T(LEFT) = T(LEFT+1)
81C  T()   [in]  Knot sequence, indexed from 1 to NT, with
82c        NT .ge. LEFT+KORDER.  Knot values must be nonincreasing.
83c        Repetition of values is permitted and has the effect of
84c        decreasing the order of contimuity at the repeated knot.
85c        Proper function representation by splines of order K is
86c        supported on the interval from T(K) to T(NT+1-K).
87c        Extrapolation can be done outside this interval.
88C  X     [in]  The abcissa at which the B-spline basis functions are to
89c        be evaluated.  The evaluation is proper if T(KORDER) .le. X .le
90c        T(NT+1-KORDER), and constitutes extrapolation otherwise.
91c  IDERIV [in]  Order of derivative requested.  Require IDERIV .ge. 0.
92c        Derivatives of order .ge. KORDER are zero.
93C  BDERIV()  [out]  On normal return, this array will contain in
94c        BDERIV(i), i = 1, ...,KORDER, the values computed by this subr.
95c        These are values at X of the derivative of order IDERIV of each
96c        of the basis functions indexed LEFT+1-KORDER through LEFT.
97C     ------------------------------------------------------------------
98c--D replaces "?": ?SBASD
99c     Both versions use IERM1, IERV1
100C     ------------------------------------------------------------------
101      integer KMAX
102      parameter(KMAX = 20)
103      integer I, I1, ID, IDERIV, ISTART, IT1, IT2
104      integer J, JP1, K1, KORDER, LEFT
105      double precision BDERIV(KORDER), DELTAL(KMAX), DELTAR(KMAX), FLK1
106      double precision SAVED, T(LEFT+KORDER), TERM, X
107C     ------------------------------------------------------------------
108      ID = IDERIV
109      if(ID .lt. 0 ) then
110         call IERM1('DSBASD',2,2,'Require IDERIV .ge. 0',
111     *   'IDERIV',ID,'.')
112         return
113      endif
114      if(ID .ge. KORDER) then
115         do 5 I = 1, KORDER
116            BDERIV(I) = 0.0d0
117    5    continue
118         return
119      endif
120      if(KORDER .gt. KMAX) then
121         call IERM1('DSBASD',1,2,'Require KORDER .le. KMAX.',
122     *   'KORDER',KORDER,',')
123         call IERV1('KMAX',KMAX,'.')
124         return
125      endif
126      ISTART = 1+ID
127      K1 = KORDER - ID
128c
129c        Evaluate K1 basis functions of order K1.  Store values in
130c        BDERIV(ISTART:ISTART+K1-1)
131c
132      BDERIV(ISTART) = 1.0d0
133      do 30 J = 1, K1-1
134         JP1 = J + 1
135         DELTAR(J) = T(LEFT+J) - X
136         DELTAL(J) = X - T(LEFT+1-J)
137         SAVED = 0.0d0
138         do 26 I=1,J
139            TERM = BDERIV(ID+I)/(DELTAR(I) + DELTAL(JP1-I))
140            BDERIV(ID+I) = SAVED + DELTAR(I)*TERM
141            SAVED = DELTAL(JP1-I)*TERM
142   26    continue
143         BDERIV(ID+JP1) = SAVED
144   30 continue
145c
146c     Loop IDERIV times, each time advancing the order of the
147c     derivative by one, so final results are values of derivatives
148c     of order IDERIV.
149c
150      do 70 I1 = ISTART, 2, -1
151         FLK1 = dble(K1)
152         IT1 = LEFT+1-I1
153         IT2 = IT1 - K1
154         do 50 I = I1, KORDER
155            BDERIV(I) = BDERIV(I) * FLK1 /(T(IT1+I) - T(IT2+I))
156   50    continue
157c
158         BDERIV(I1-1) = -BDERIV(I1)
159         do 60 I = I1, KORDER-1
160            BDERIV(I) = BDERIV(I) - BDERIV(I+1)
161   60    continue
162         K1 = K1 + 1
163   70 continue
164      return
165      end
166