← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/sandbox/spline/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:
scipy/sandbox/spline/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