← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpspgr.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/fpspgr.f
1
 
      subroutine fpspgr(iopt,ider,u,mu,v,mv,r,mr,r0,r1,s,nuest,nvest,
2
 
     * tol,maxit,nc,nu,tu,nv,tv,c,fp,fp0,fpold,reducu,reducv,fpintu,
3
 
     * fpintv,dr,step,lastdi,nplusu,nplusv,lastu0,lastu1,nru,nrv,
4
 
     * nrdatu,nrdatv,wrk,lwrk,ier)
5
 
c  ..
6
 
c  ..scalar arguments..
7
 
      integer mu,mv,mr,nuest,nvest,maxit,nc,nu,nv,lastdi,nplusu,nplusv,
8
 
     * lastu0,lastu1,lwrk,ier
9
 
      real*8 r0,r1,s,tol,fp,fp0,fpold,reducu,reducv
10
 
c  ..array arguments..
11
 
      integer iopt(3),ider(4),nrdatu(nuest),nrdatv(nvest),nru(mu),
12
 
     * nrv(mv)
13
 
      real*8 u(mu),v(mv),r(mr),tu(nuest),tv(nvest),c(nc),fpintu(nuest),
14
 
     * fpintv(nvest),dr(6),wrk(lwrk),step(2)
15
 
c  ..local scalars..
16
 
      real*8 acc,fpms,f1,f2,f3,p,per,pi,p1,p2,p3,vb,ve,rmax,rmin,rn,one,
17
 
     *
18
 
     * con1,con4,con9
19
 
      integer i,ich1,ich3,ifbu,ifbv,ifsu,ifsv,istart,iter,i1,i2,j,ju,
20
 
     * ktu,l,l1,l2,l3,l4,mpm,mumin,mu0,mu1,nn,nplu,nplv,npl1,nrintu,
21
 
     * nrintv,nue,numax,nve,nvmax
22
 
c  ..local arrays..
23
 
      integer idd(4)
24
 
      real*8 drr(6)
25
 
c  ..function references..
26
 
      real*8 abs,datan2,fprati
27
 
      integer max0,min0
28
 
c  ..subroutine references..
29
 
c    fpknot,fpopsp
30
 
c  ..
31
 
c   set constants
32
 
      one = 1d0
33
 
      con1 = 0.1e0
34
 
      con9 = 0.9e0
35
 
      con4 = 0.4e-01
36
 
c   initialization
37
 
      ifsu = 0
38
 
      ifsv = 0
39
 
      ifbu = 0
40
 
      ifbv = 0
41
 
      p = -one
42
 
      mumin = 4
43
 
      if(ider(1).ge.0) mumin = mumin-1
44
 
      if(iopt(2).eq.1 .and. ider(2).eq.1) mumin = mumin-1
45
 
      if(ider(3).ge.0) mumin = mumin-1
46
 
      if(iopt(3).eq.1 .and. ider(4).eq.1) mumin = mumin-1
47
 
      if(mumin.eq.0) mumin = 1
48
 
      pi = datan2(0d0,-one)
49
 
      per = pi+pi
50
 
      vb = v(1)
51
 
      ve = vb+per
52
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
53
 
c part 1: determination of the number of knots and their position.     c
54
 
c ****************************************************************     c
55
 
c  given a set of knots we compute the least-squares spline sinf(u,v)  c
56
 
c  and the corresponding sum of squared residuals fp = f(p=inf).       c
57
 
c  if iopt(1)=-1  sinf(u,v) is the requested approximation.            c
58
 
c  if iopt(1)>=0  we check whether we can accept the knots:            c
59
 
c    if fp <= s we will continue with the current set of knots.        c
60
 
c    if fp >  s we will increase the number of knots and compute the   c
61
 
c       corresponding least-squares spline until finally fp <= s.      c
62
 
c    the initial choice of knots depends on the value of s and iopt.   c
63
 
c    if s=0 we have spline interpolation; in that case the number of   c
64
 
c     knots in the u-direction equals nu=numax=mu+6+iopt(2)+iopt(3)    c
65
 
c     and in the v-direction nv=nvmax=mv+7.                            c
66
 
c    if s>0 and                                                        c
67
 
c      iopt(1)=0 we first compute the least-squares polynomial,i.e. a  c
68
 
c       spline without interior knots : nu=8 ; nv=8.                   c
69
 
c      iopt(1)=1 we start with the set of knots found at the last call c
70
 
c       of the routine, except for the case that s > fp0; then we      c
71
 
c       compute the least-squares polynomial directly.                 c
72
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
73
 
      if(iopt(1).lt.0) go to 120
74
 
c  acc denotes the absolute tolerance for the root of f(p)=s.
75
 
      acc = tol*s
76
 
c  numax and nvmax denote the number of knots needed for interpolation.
77
 
      numax = mu+6+iopt(2)+iopt(3)
78
 
      nvmax = mv+7
79
 
      nue = min0(numax,nuest)
80
 
      nve = min0(nvmax,nvest)
81
 
      if(s.gt.0.) go to 100
82
 
c  if s = 0, s(u,v) is an interpolating spline.
83
 
      nu = numax
84
 
      nv = nvmax
85
 
c  test whether the required storage space exceeds the available one.
86
 
      if(nu.gt.nuest .or. nv.gt.nvest) go to 420
87
 
c  find the position of the knots in the v-direction.
88
 
      do 10 l=1,mv
89
 
        tv(l+3) = v(l)
90
 
  10  continue
91
 
      tv(mv+4) = ve
92
 
      l1 = mv-2
93
 
      l2 = mv+5
94
 
      do 20 i=1,3
95
 
         tv(i) = v(l1)-per
96
 
         tv(l2) = v(i+1)+per
97
 
         l1 = l1+1
98
 
         l2 = l2+1
99
 
  20  continue
100
 
c  if not all the derivative values g(i,j) are given, we will first
101
 
c  estimate these values by computing a least-squares spline
102
 
      idd(1) = ider(1)
103
 
      if(idd(1).eq.0) idd(1) = 1
104
 
      if(idd(1).gt.0) dr(1) = r0
105
 
      idd(2) = ider(2)
106
 
      idd(3) = ider(3)
107
 
      if(idd(3).eq.0) idd(3) = 1
108
 
      if(idd(3).gt.0) dr(4) = r1
109
 
      idd(4) = ider(4)
110
 
      if(ider(1).lt.0 .or. ider(3).lt.0) go to 30
111
 
      if(iopt(2).ne.0 .and. ider(2).eq.0) go to 30
112
 
      if(iopt(3).eq.0 .or. ider(4).ne.0) go to 70
113
 
c we set up the knots in the u-direction for computing the least-squares
114
 
c spline.
115
 
  30  i1 = 3
116
 
      i2 = mu-2
117
 
      nu = 4
118
 
      do 40 i=1,mu
119
 
         if(i1.gt.i2) go to 50
120
 
         nu = nu+1
121
 
         tu(nu) = u(i1)
122
 
         i1 = i1+2
123
 
  40  continue
124
 
  50  do 60 i=1,4
125
 
         tu(i) = 0.
126
 
         nu = nu+1
127
 
         tu(nu) = pi
128
 
  60  continue
129
 
c we compute the least-squares spline for estimating the derivatives.
130
 
      call fpopsp(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,r,mr,r0,r1,dr,iopt,idd,
131
 
     *  tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpintu,fpintv,nru,nrv,
132
 
     *  wrk,lwrk)
133
 
      ifsu = 0
134
 
c if all the derivatives at the origin are known, we compute the
135
 
c interpolating spline.
136
 
c we set up the knots in the u-direction, needed for interpolation.
137
 
  70  nn = numax-8
138
 
      if(nn.eq.0) go to 95
139
 
      ju = 2-iopt(2)
140
 
      do 80 l=1,nn
141
 
        tu(l+4) = u(ju)
142
 
        ju = ju+1
143
 
  80  continue
144
 
      nu = numax
145
 
      l = nu
146
 
      do 90 i=1,4
147
 
         tu(i) = 0.
148
 
         tu(l) = pi
149
 
         l = l-1
150
 
  90  continue
151
 
c we compute the interpolating spline.
152
 
  95  call fpopsp(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,r,mr,r0,r1,dr,iopt,idd,
153
 
     *  tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpintu,fpintv,nru,nrv,
154
 
     *  wrk,lwrk)
155
 
      go to 430
156
 
c  if s>0 our initial choice of knots depends on the value of iopt(1).
157
 
 100  ier = 0
158
 
      if(iopt(1).eq.0) go to 115
159
 
      step(1) = -step(1)
160
 
      step(2) = -step(2)
161
 
      if(fp0.le.s) go to 115
162
 
c  if iopt(1)=1 and fp0 > s we start computing the least-squares spline
163
 
c  according to the set of knots found at the last call of the routine.
164
 
c  we determine the number of grid coordinates u(i) inside each knot
165
 
c  interval (tu(l),tu(l+1)).
166
 
      l = 5
167
 
      j = 1
168
 
      nrdatu(1) = 0
169
 
      mu0 = 2-iopt(2)
170
 
      mu1 = mu-1+iopt(3)
171
 
      do 105 i=mu0,mu1
172
 
        nrdatu(j) = nrdatu(j)+1
173
 
        if(u(i).lt.tu(l)) go to 105
174
 
        nrdatu(j) = nrdatu(j)-1
175
 
        l = l+1
176
 
        j = j+1
177
 
        nrdatu(j) = 0
178
 
 105  continue
179
 
c  we determine the number of grid coordinates v(i) inside each knot
180
 
c  interval (tv(l),tv(l+1)).
181
 
      l = 5
182
 
      j = 1
183
 
      nrdatv(1) = 0
184
 
      do 110 i=2,mv
185
 
        nrdatv(j) = nrdatv(j)+1
186
 
        if(v(i).lt.tv(l)) go to 110
187
 
        nrdatv(j) = nrdatv(j)-1
188
 
        l = l+1
189
 
        j = j+1
190
 
        nrdatv(j) = 0
191
 
 110  continue
192
 
      idd(1) = ider(1)
193
 
      idd(2) = ider(2)
194
 
      idd(3) = ider(3)
195
 
      idd(4) = ider(4)
196
 
      go to 120
197
 
c  if iopt(1)=0 or iopt(1)=1 and s >= fp0,we start computing the least-
198
 
c  squares polynomial (which is a spline without interior knots).
199
 
 115  ier = -2
200
 
      idd(1) = ider(1)
201
 
      idd(2) = 1
202
 
      idd(3) = ider(3)
203
 
      idd(4) = 1
204
 
      nu = 8
205
 
      nv = 8
206
 
      nrdatu(1) = mu-2+iopt(2)+iopt(3)
207
 
      nrdatv(1) = mv-1
208
 
      lastdi = 0
209
 
      nplusu = 0
210
 
      nplusv = 0
211
 
      fp0 = 0.
212
 
      fpold = 0.
213
 
      reducu = 0.
214
 
      reducv = 0.
215
 
c  main loop for the different sets of knots.mpm=mu+mv is a save upper
216
 
c  bound for the number of trials.
217
 
 120  mpm = mu+mv
218
 
      do 270 iter=1,mpm
219
 
c  find nrintu (nrintv) which is the number of knot intervals in the
220
 
c  u-direction (v-direction).
221
 
        nrintu = nu-7
222
 
        nrintv = nv-7
223
 
c  find the position of the additional knots which are needed for the
224
 
c  b-spline representation of s(u,v).
225
 
        i = nu
226
 
        do 125 j=1,4
227
 
          tu(j) = 0.
228
 
          tu(i) = pi
229
 
          i = i-1
230
 
 125    continue
231
 
        l1 = 4
232
 
        l2 = l1
233
 
        l3 = nv-3
234
 
        l4 = l3
235
 
        tv(l2) = vb
236
 
        tv(l3) = ve
237
 
        do 130 j=1,3
238
 
          l1 = l1+1
239
 
          l2 = l2-1
240
 
          l3 = l3+1
241
 
          l4 = l4-1
242
 
          tv(l2) = tv(l4)-per
243
 
          tv(l3) = tv(l1)+per
244
 
 130    continue
245
 
c  find an estimate of the range of possible values for the optimal
246
 
c  derivatives at the origin.
247
 
        ktu = nrdatu(1)+2-iopt(2)
248
 
        if(ktu.lt.mumin) ktu = mumin
249
 
        if(ktu.eq.lastu0) go to 140
250
 
         rmin = r0
251
 
         rmax = r0
252
 
         l = mv*ktu
253
 
         do 135 i=1,l
254
 
            if(r(i).lt.rmin) rmin = r(i)
255
 
            if(r(i).gt.rmax) rmax = r(i)
256
 
 135     continue
257
 
         step(1) = rmax-rmin
258
 
         lastu0 = ktu
259
 
 140    ktu = nrdatu(nrintu)+2-iopt(3)
260
 
        if(ktu.lt.mumin) ktu = mumin
261
 
        if(ktu.eq.lastu1) go to 150
262
 
         rmin = r1
263
 
         rmax = r1
264
 
         l = mv*ktu
265
 
         j = mr
266
 
         do 145 i=1,l
267
 
            if(r(j).lt.rmin) rmin = r(j)
268
 
            if(r(j).gt.rmax) rmax = r(j)
269
 
            j = j-1
270
 
 145     continue
271
 
         step(2) = rmax-rmin
272
 
         lastu1 = ktu
273
 
c  find the least-squares spline sinf(u,v).
274
 
 150    call fpopsp(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,r,mr,r0,r1,dr,iopt,
275
 
     *   idd,tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpintu,fpintv,nru,
276
 
     *   nrv,wrk,lwrk)
277
 
        if(step(1).lt.0.) step(1) = -step(1)
278
 
        if(step(2).lt.0.) step(2) = -step(2)
279
 
        if(ier.eq.(-2)) fp0 = fp
280
 
c  test whether the least-squares spline is an acceptable solution.
281
 
        if(iopt(1).lt.0) go to 440
282
 
        fpms = fp-s
283
 
        if(abs(fpms) .lt. acc) go to 440
284
 
c  if f(p=inf) < s, we accept the choice of knots.
285
 
        if(fpms.lt.0.) go to 300
286
 
c  if nu=numax and nv=nvmax, sinf(u,v) is an interpolating spline
287
 
        if(nu.eq.numax .and. nv.eq.nvmax) go to 430
288
 
c  increase the number of knots.
289
 
c  if nu=nue and nv=nve we cannot further increase the number of knots
290
 
c  because of the storage capacity limitation.
291
 
        if(nu.eq.nue .and. nv.eq.nve) go to 420
292
 
        if(ider(1).eq.0) fpintu(1) = fpintu(1)+(r0-dr(1))**2
293
 
        if(ider(3).eq.0) fpintu(nrintu) = fpintu(nrintu)+(r1-dr(4))**2
294
 
        ier = 0
295
 
c  adjust the parameter reducu or reducv according to the direction
296
 
c  in which the last added knots were located.
297
 
        if (lastdi.lt.0) go to 160
298
 
        if (lastdi.eq.0) go to 155
299
 
        go to 170
300
 
 155     nplv = 3
301
 
         idd(2) = ider(2)
302
 
         idd(4) = ider(4)
303
 
         fpold = fp
304
 
         go to 230
305
 
 160    reducu = fpold-fp
306
 
        go to 175
307
 
 170    reducv = fpold-fp
308
 
c  store the sum of squared residuals for the current set of knots.
309
 
 175    fpold = fp
310
 
c  find nplu, the number of knots we should add in the u-direction.
311
 
        nplu = 1
312
 
        if(nu.eq.8) go to 180
313
 
        npl1 = nplusu*2
314
 
        rn = nplusu
315
 
        if(reducu.gt.acc) npl1 = rn*fpms/reducu
316
 
        nplu = min0(nplusu*2,max0(npl1,nplusu/2,1))
317
 
c  find nplv, the number of knots we should add in the v-direction.
318
 
 180    nplv = 3
319
 
        if(nv.eq.8) go to 190
320
 
        npl1 = nplusv*2
321
 
        rn = nplusv
322
 
        if(reducv.gt.acc) npl1 = rn*fpms/reducv
323
 
        nplv = min0(nplusv*2,max0(npl1,nplusv/2,1))
324
 
c  test whether we are going to add knots in the u- or v-direction.
325
 
 190    if (nplu.lt.nplv) go to 210
326
 
        if (nplu.eq.nplv) go to 200
327
 
        go to 230
328
 
 200    if(lastdi.lt.0) go to 230
329
 
 210    if(nu.eq.nue) go to 230
330
 
c  addition in the u-direction.
331
 
        lastdi = -1
332
 
        nplusu = nplu
333
 
        ifsu = 0
334
 
        istart = 0
335
 
        if(iopt(2).eq.0) istart = 1
336
 
        do 220 l=1,nplusu
337
 
c  add a new knot in the u-direction
338
 
          call fpknot(u,mu,tu,nu,fpintu,nrdatu,nrintu,nuest,istart)
339
 
c  test whether we cannot further increase the number of knots in the
340
 
c  u-direction.
341
 
          if(nu.eq.nue) go to 270
342
 
 220    continue
343
 
        go to 270
344
 
 230    if(nv.eq.nve) go to 210
345
 
c  addition in the v-direction.
346
 
        lastdi = 1
347
 
        nplusv = nplv
348
 
        ifsv = 0
349
 
        do 240 l=1,nplusv
350
 
c  add a new knot in the v-direction.
351
 
          call fpknot(v,mv,tv,nv,fpintv,nrdatv,nrintv,nvest,1)
352
 
c  test whether we cannot further increase the number of knots in the
353
 
c  v-direction.
354
 
          if(nv.eq.nve) go to 270
355
 
 240    continue
356
 
c  restart the computations with the new set of knots.
357
 
 270  continue
358
 
c  test whether the least-squares polynomial is a solution of our
359
 
c  approximation problem.
360
 
 300  if(ier.eq.(-2)) go to 440
361
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
362
 
c part 2: determination of the smoothing spline sp(u,v)                c
363
 
c *****************************************************                c
364
 
c  we have determined the number of knots and their position. we now   c
365
 
c  compute the b-spline coefficients of the smoothing spline sp(u,v).  c
366
 
c  this smoothing spline depends on the parameter p in such a way that c
367
 
c    f(p) = sumi=1,mu(sumj=1,mv((z(i,j)-sp(u(i),v(j)))**2)             c
368
 
c  is a continuous, strictly decreasing function of p. moreover the    c
369
 
c  least-squares polynomial corresponds to p=0 and the least-squares   c
370
 
c  spline to p=infinity. then iteratively we have to determine the     c
371
 
c  positive value of p such that f(p)=s. the process which is proposed c
372
 
c  here makes use of rational interpolation. f(p) is approximated by a c
373
 
c  rational function r(p)=(u*p+v)/(p+w); three values of p (p1,p2,p3)  c
374
 
c  with corresponding values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s)c
375
 
c  are used to calculate the new value of p such that r(p)=s.          c
376
 
c  convergence is guaranteed by taking f1 > 0 and f3 < 0.              c
377
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
378
 
c  initial value for p.
379
 
      p1 = 0.
380
 
      f1 = fp0-s
381
 
      p3 = -one
382
 
      f3 = fpms
383
 
      p = one
384
 
      do 305 i=1,6
385
 
        drr(i) = dr(i)
386
 
 305  continue
387
 
      ich1 = 0
388
 
      ich3 = 0
389
 
c  iteration process to find the root of f(p)=s.
390
 
      do 350 iter = 1,maxit
391
 
c  find the smoothing spline sp(u,v) and the corresponding sum f(p).
392
 
        call fpopsp(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,r,mr,r0,r1,drr,iopt,
393
 
     *   idd,tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpintu,fpintv,nru,
394
 
     *   nrv,wrk,lwrk)
395
 
c  test whether the approximation sp(u,v) is an acceptable solution.
396
 
        fpms = fp-s
397
 
        if(abs(fpms).lt.acc) go to 440
398
 
c  test whether the maximum allowable number of iterations has been
399
 
c  reached.
400
 
        if(iter.eq.maxit) go to 400
401
 
c  carry out one more step of the iteration process.
402
 
        p2 = p
403
 
        f2 = fpms
404
 
        if(ich3.ne.0) go to 320
405
 
        if((f2-f3).gt.acc) go to 310
406
 
c  our initial choice of p is too large.
407
 
        p3 = p2
408
 
        f3 = f2
409
 
        p = p*con4
410
 
        if(p.le.p1) p = p1*con9 + p2*con1
411
 
        go to 350
412
 
 310    if(f2.lt.0.) ich3 = 1
413
 
 320    if(ich1.ne.0) go to 340
414
 
        if((f1-f2).gt.acc) go to 330
415
 
c  our initial choice of p is too small
416
 
        p1 = p2
417
 
        f1 = f2
418
 
        p = p/con4
419
 
        if(p3.lt.0.) go to 350
420
 
        if(p.ge.p3) p = p2*con1 + p3*con9
421
 
        go to 350
422
 
c  test whether the iteration process proceeds as theoretically
423
 
c  expected.
424
 
 330    if(f2.gt.0.) ich1 = 1
425
 
 340    if(f2.ge.f1 .or. f2.le.f3) go to 410
426
 
c  find the new value of p.
427
 
        p = fprati(p1,f1,p2,f2,p3,f3)
428
 
 350  continue
429
 
c  error codes and messages.
430
 
 400  ier = 3
431
 
      go to 440
432
 
 410  ier = 2
433
 
      go to 440
434
 
 420  ier = 1
435
 
      go to 440
436
 
 430  ier = -1
437
 
      fp = 0.
438
 
 440  return
439
 
      end