1      subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,
2     * wrk,lwrk,iwrk,kwrk,ier)
3c  subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
4c  ,my the partial derivative ( order nux,nuy) of a bivariate spline
5c  s(x,y) of degrees kx and ky, given in the b-spline representation.
6c
7c  calling sequence:
8c     call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
9c    * iwrk,kwrk,ier)
10c
11c  input parameters:
12c   tx    : real array, length nx, which contains the position of the
13c           knots in the x-direction.
14c   nx    : integer, giving the total number of knots in the x-direction
15c   ty    : real array, length ny, which contains the position of the
16c           knots in the y-direction.
17c   ny    : integer, giving the total number of knots in the y-direction
18c   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
19c           b-spline coefficients.
20c   kx,ky : integer values, giving the degrees of the spline.
21c   nux   : integer values, specifying the order of the partial
22c   nuy     derivative. 0<=nux<kx, 0<=nuy<ky.
23c   x     : real array of dimension (mx).
24c           before entry x(i) must be set to the x co-ordinate of the
25c           i-th grid point along the x-axis.
26c           tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
27c   mx    : on entry mx must specify the number of grid points along
28c           the x-axis. mx >=1.
29c   y     : real array of dimension (my).
30c           before entry y(j) must be set to the y co-ordinate of the
31c           j-th grid point along the y-axis.
32c           ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
33c   my    : on entry my must specify the number of grid points along
34c           the y-axis. my >=1.
35c   wrk   : real array of dimension lwrk. used as workspace.
36c   lwrk  : integer, specifying the dimension of wrk.
37c           lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
38c   iwrk  : integer array of dimension kwrk. used as workspace.
39c   kwrk  : integer, specifying the dimension of iwrk. kwrk >= mx+my.
40c
41c  output parameters:
42c   z     : real array of dimension (mx*my).
43c           on successful exit z(my*(i-1)+j) contains the value of the
44c           specified partial derivative of s(x,y) at the point
45c           (x(i),y(j)),i=1,...,mx;j=1,...,my.
46c   ier   : integer error flag
47c    ier=0 : normal return
48c    ier=10: invalid input data (see restrictions)
49c
50c  restrictions:
51c   mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my
52c   lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
53c   tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
54c   ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
55c
56c  other subroutines required:
57c    fpbisp,fpbspl
58c
59c  references :
60c    de boor c : on calculating with b-splines, j. approximation theory
61c                6 (1972) 50-62.
62c   dierckx p. : curve and surface fitting with splines, monographs on
63c                numerical analysis, oxford university press, 1993.
64c
65c  author :
66c    p.dierckx
67c    dept. computer science, k.u.leuven
68c    celestijnenlaan 200a, b-3001 heverlee, belgium.
69c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
70c
71c  latest update : march 1989
72c
73c  ..scalar arguments..
74      integer nx,ny,kx,ky,nux,nuy,mx,my,lwrk,kwrk,ier
75c  ..array arguments..
76      integer iwrk(kwrk)
77      real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my),
78     * wrk(lwrk)
79c  ..local scalars..
80      integer i,iwx,iwy,j,kkx,kky,kx1,ky1,lx,ly,lwest,l1,l2,m,m0,m1,
81     * nc,nkx1,nky1,nxx,nyy
82      real*8 ak,fac
83c  ..
84c  before starting computations a data check is made. if the input data
85c  are invalid control is immediately repassed to the calling program.
86      ier = 10
87      kx1 = kx+1
88      ky1 = ky+1
89      nkx1 = nx-kx1
90      nky1 = ny-ky1
91      nc = nkx1*nky1
92      if(nux.lt.0 .or. nux.ge.kx) go to 400
93      if(nuy.lt.0 .or. nuy.ge.ky) go to 400
94      lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my
95      if(lwrk.lt.lwest) go to 400
96      if(kwrk.lt.(mx+my)) go to 400
97      if (mx.lt.1) go to 400
98      if (mx.eq.1) go to 30
99      go to 10
100  10  do 20 i=2,mx
101        if(x(i).lt.x(i-1)) go to 400
102  20  continue
103  30  if (my.lt.1) go to 400
104      if (my.eq.1) go to 60
105      go to 40
106  40  do 50 i=2,my
107        if(y(i).lt.y(i-1)) go to 400
108  50  continue
109  60  ier = 0
110      nxx = nkx1
111      nyy = nky1
112      kkx = kx
113      kky = ky
114c  the partial derivative of order (nux,nuy) of a bivariate spline of
115c  degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
116c  we calculate the b-spline coefficients of this spline
117      do 70 i=1,nc
118        wrk(i) = c(i)
119  70  continue
120      if(nux.eq.0) go to 200
121      lx = 1
122      do 100 j=1,nux
123        ak = kkx
124        nxx = nxx-1
125        l1 = lx
126        m0 = 1
127        do 90 i=1,nxx
128          l1 = l1+1
129          l2 = l1+kkx
130          fac = tx(l2)-tx(l1)
131          if(fac.le.0.) go to 90
132          do 80 m=1,nyy
133            m1 = m0+nyy
134            wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
135            m0  = m0+1
136  80      continue
137  90    continue
138        lx = lx+1
139        kkx = kkx-1
140 100  continue
141 200  if(nuy.eq.0) go to 300
142      ly = 1
143      do 230 j=1,nuy
144        ak = kky
145        nyy = nyy-1
146        l1 = ly
147        do 220 i=1,nyy
148          l1 = l1+1
149          l2 = l1+kky
150          fac = ty(l2)-ty(l1)
151          if(fac.le.0.) go to 220
152          m0 = i
153          do 210 m=1,nxx
154            m1 = m0+1
155            wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
156            m0  = m0+nky1
157 210      continue
158 220    continue
159        ly = ly+1
160        kky = kky-1
161 230  continue
162      m0 = nyy
163      m1 = nky1
164      do 250 m=2,nxx
165        do 240 i=1,nyy
166          m0 = m0+1
167          m1 = m1+1
168          wrk(m0) = wrk(m1)
169 240    continue
170        m1 = m1+nuy
171 250  continue
172c  we partition the working space and evaluate the partial derivative
173 300  iwx = 1+nxx*nyy
174      iwy = iwx+mx*(kx1-nux)
175      call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky,
176     * x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1))
177 400  return
178      end
179
180