1      subroutine splder(t,n,c,k,nu,x,y,m,e,wrk,ier)
2      implicit none
3c  subroutine splder evaluates in a number of points x(i),i=1,2,...,m
4c  the derivative of order nu of a spline s(x) of degree k,given in
5c  its b-spline representation.
6c
7c  calling sequence:
8c     call splder(t,n,c,k,nu,x,y,m,e,wrk,ier)
9c
10c  input parameters:
11c    t    : array,length n, which contains the position of the knots.
12c    n    : integer, giving the total number of knots of s(x).
13c    c    : array,length n, which contains the b-spline coefficients.
14c    k    : integer, giving the degree of s(x).
15c    nu   : integer, specifying the order of the derivative. 0<=nu<=k
16c    x    : array,length m, which contains the points where the deriv-
17c           ative of s(x) must be evaluated.
18c    m    : integer, giving the number of points where the derivative
19c           of s(x) must be evaluated
20c    e    : integer, if 0 the spline is extrapolated from the end
21c           spans for points not in the support, if 1 the spline
22c           evaluates to zero for those points, and if 2 ier is set to
23c           1 and the subroutine returns.
24c    wrk  : real array of dimension n. used as working space.
25c
26c  output parameters:
27c    y    : array,length m, giving the value of the derivative of s(x)
28c           at the different points.
29c    ier  : error flag
30c      ier = 0 : normal return
31c      ier = 1 : argument out of bounds and e == 2
32c      ier =10 : invalid input data (see restrictions)
33c
34c  restrictions:
35c    0 <= nu <= k
36c    m >= 1
37c    t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
38c
39c  other subroutines required: fpbspl
40c
41c  references :
42c    de boor c : on calculating with b-splines, j. approximation theory
43c                6 (1972) 50-62.
44c    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths
45c                applics 10 (1972) 134-149.
46c   dierckx p. : curve and surface fitting with splines, monographs on
47c                numerical analysis, oxford university press, 1993.
48c
49c  author :
50c    p.dierckx
51c    dept. computer science, k.u.leuven
52c    celestijnenlaan 200a, b-3001 heverlee, belgium.
53c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
54c
55c  latest update : march 1987
56c
57c++ pearu: 13 aug 20003
58c++   - disabled cliping x values to interval [min(t),max(t)]
59c++   - removed the restriction of the orderness of x values
60c++   - fixed initialization of sp to double precision value
61c
62c  ..scalar arguments..
63      integer n,k,nu,m,e,ier
64c  ..array arguments..
65      real*8 t(n),c(n),x(m),y(m),wrk(n)
66c  ..local scalars..
67      integer i,j,kk,k1,k2,l,ll,l1,l2,nk1,nk2,nn
68      real*8 ak,arg,fac,sp,tb,te
69c++..
70      integer k3
71c..++
72c  ..local arrays ..
73      real*8 h(6)
74c  before starting computations a data check is made. if the input data
75c  are invalid control is immediately repassed to the calling program.
76      ier = 10
77      if(nu.lt.0 .or. nu.gt.k) go to 200
78c--      if(m-1) 200,30,10
79c++..
80      if(m.lt.1) go to 200
81c..++
82c--  10  do 20 i=2,m
83c--        if(x(i).lt.x(i-1)) go to 200
84c--  20  continue
85      ier = 0
86c  fetch tb and te, the boundaries of the approximation interval.
87      k1 = k+1
88      k3 = k1+1
89      nk1 = n-k1
90      tb = t(k1)
91      te = t(nk1+1)
92c  the derivative of order nu of a spline of degree k is a spline of
93c  degree k-nu,the b-spline coefficients wrk(i) of which can be found
94c  using the recurrence scheme of de boor.
95      l = 1
96      kk = k
97      nn = n
98      do 40 i=1,nk1
99         wrk(i) = c(i)
100  40  continue
101      if(nu.eq.0) go to 100
102      nk2 = nk1
103      do 60 j=1,nu
104         ak = kk
105         nk2 = nk2-1
106         l1 = l
107         do 50 i=1,nk2
108            l1 = l1+1
109            l2 = l1+kk
110            fac = t(l2)-t(l1)
111            if(fac.le.0.) go to 50
112            wrk(i) = ak*(wrk(i+1)-wrk(i))/fac
113  50     continue
114         l = l+1
115         kk = kk-1
116  60  continue
117      if(kk.ne.0) go to 100
118c  if nu=k the derivative is a piecewise constant function
119      j = 1
120      do 90 i=1,m
121        arg = x(i)
122c++..
123c  check if arg is in the support
124        if (arg .lt. tb .or. arg .gt. te) then
125            if (e .eq. 0) then
126                goto 65
127            else if (e .eq. 1) then
128                y(i) = 0
129                goto 90
130            else if (e .eq. 2) then
131                ier = 1
132                goto 200
133            endif
134        endif
135c  search for knot interval t(l) <= arg < t(l+1)
136 65     if(arg.ge.t(l) .or. l+1.eq.k3) go to 70
137        l1 = l
138        l = l-1
139        j = j-1
140        go to 65
141c..++
142  70    if(arg.lt.t(l+1) .or. l.eq.nk1) go to 80
143        l = l+1
144        j = j+1
145        go to 70
146  80    y(i) = wrk(j)
147  90  continue
148      go to 200
149
150 100  l = k1
151      l1 = l+1
152      k2 = k1-nu
153c  main loop for the different points.
154      do 180 i=1,m
155c  fetch a new x-value arg.
156        arg = x(i)
157c  check if arg is in the support
158        if (arg .lt. tb .or. arg .gt. te) then
159            if (e .eq. 0) then
160                goto 135
161            else if (e .eq. 1) then
162                y(i) = 0
163                goto 180
164            else if (e .eq. 2) then
165                ier = 1
166                goto 200
167            endif
168        endif
169c  search for knot interval t(l) <= arg < t(l+1)
170 135    if(arg.ge.t(l) .or. l1.eq.k3) go to 140
171        l1 = l
172        l = l-1
173        go to 135
174c..++
175 140    if(arg.lt.t(l1) .or. l.eq.nk1) go to 150
176        l = l1
177        l1 = l+1
178        go to 140
179c  evaluate the non-zero b-splines of degree k-nu at arg.
180 150    call fpbspl(t,n,kk,arg,l,h)
181c  find the value of the derivative at x=arg.
182        sp = 0.0d0
183        ll = l-k1
184        do 160 j=1,k2
185          ll = ll+1
186          sp = sp+wrk(ll)*h(j)
187 160    continue
188        y(i) = sp
189 180  continue
190 200  return
191      end
192