← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/parder.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/interpolate/fitpack/parder.f
1
 
      subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,
2
 
     * wrk,lwrk,iwrk,kwrk,ier)
3
 
c  subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
4
 
c  ,my the partial derivative ( order nux,nuy) of a bivariate spline
5
 
c  s(x,y) of degrees kx and ky, given in the b-spline representation.
6
 
c
7
 
c  calling sequence:
8
 
c     call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
9
 
c    * iwrk,kwrk,ier)
10
 
c
11
 
c  input parameters:
12
 
c   tx    : real array, length nx, which contains the position of the
13
 
c           knots in the x-direction.
14
 
c   nx    : integer, giving the total number of knots in the x-direction
15
 
c   ty    : real array, length ny, which contains the position of the
16
 
c           knots in the y-direction.
17
 
c   ny    : integer, giving the total number of knots in the y-direction
18
 
c   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
19
 
c           b-spline coefficients.
20
 
c   kx,ky : integer values, giving the degrees of the spline.
21
 
c   nux   : integer values, specifying the order of the partial
22
 
c   nuy     derivative. 0<=nux<kx, 0<=nuy<ky.
23
 
c   x     : real array of dimension (mx).
24
 
c           before entry x(i) must be set to the x co-ordinate of the
25
 
c           i-th grid point along the x-axis.
26
 
c           tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
27
 
c   mx    : on entry mx must specify the number of grid points along
28
 
c           the x-axis. mx >=1.
29
 
c   y     : real array of dimension (my).
30
 
c           before entry y(j) must be set to the y co-ordinate of the
31
 
c           j-th grid point along the y-axis.
32
 
c           ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
33
 
c   my    : on entry my must specify the number of grid points along
34
 
c           the y-axis. my >=1.
35
 
c   wrk   : real array of dimension lwrk. used as workspace.
36
 
c   lwrk  : integer, specifying the dimension of wrk.
37
 
c           lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
38
 
c   iwrk  : integer array of dimension kwrk. used as workspace.
39
 
c   kwrk  : integer, specifying the dimension of iwrk. kwrk >= mx+my.
40
 
c
41
 
c  output parameters:
42
 
c   z     : real array of dimension (mx*my).
43
 
c           on succesful exit z(my*(i-1)+j) contains the value of the
44
 
c           specified partial derivative of s(x,y) at the point
45
 
c           (x(i),y(j)),i=1,...,mx;j=1,...,my.
46
 
c   ier   : integer error flag
47
 
c    ier=0 : normal return
48
 
c    ier=10: invalid input data (see restrictions)
49
 
c
50
 
c  restrictions:
51
 
c   mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my
52
 
c   lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
53
 
c   tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
54
 
c   ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
55
 
c
56
 
c  other subroutines required:
57
 
c    fpbisp,fpbspl
58
 
c
59
 
c  references :
60
 
c    de boor c : on calculating with b-splines, j. approximation theory
61
 
c                6 (1972) 50-62.
62
 
c   dierckx p. : curve and surface fitting with splines, monographs on
63
 
c                numerical analysis, oxford university press, 1993.
64
 
c
65
 
c  author :
66
 
c    p.dierckx
67
 
c    dept. computer science, k.u.leuven
68
 
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
69
 
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
70
 
c
71
 
c  latest update : march 1989
72
 
c
73
 
c  ..scalar arguments..
74
 
      integer nx,ny,kx,ky,nux,nuy,mx,my,lwrk,kwrk,ier
75
 
c  ..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)
79
 
c  ..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
83
 
c  ..
84
 
c  before starting computations a data check is made. if the input data
85
 
c  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
114
 
c  the partial derivative of order (nux,nuy) of a bivariate spline of
115
 
c  degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
116
 
c  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
172
 
c  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