← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpcurf.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/fpcurf.f
1
 
      subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,
2
 
     * n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
3
 
c  ..
4
 
c  ..scalar arguments..
5
 
      real*8 xb,xe,s,tol,fp
6
 
      integer iopt,m,k,nest,maxit,k1,k2,n,ier
7
 
c  ..array arguments..
8
 
      real*8 x(m),y(m),w(m),t(nest),c(nest),fpint(nest),
9
 
     * z(nest),a(nest,k1),b(nest,k2),g(nest,k2),q(m,k1)
10
 
      integer nrdata(nest)
11
 
c  ..local scalars..
12
 
      real*8 acc,con1,con4,con9,cos,half,fpart,fpms,fpold,fp0,f1,f2,f3,
13
 
     * one,p,pinv,piv,p1,p2,p3,rn,sin,store,term,wi,xi,yi
14
 
      integer i,ich1,ich3,it,iter,i1,i2,i3,j,k3,l,l0,
15
 
     * mk1,new,nk1,nmax,nmin,nplus,npl1,nrint,n8
16
 
c  ..local arrays..
17
 
      real*8 h(7)
18
 
c  ..function references
19
 
      real*8 abs,fprati
20
 
      integer max0,min0
21
 
c  ..subroutine references..
22
 
c    fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
23
 
c  ..
24
 
c  set constants
25
 
      one = 0.1d+01
26
 
      con1 = 0.1d0
27
 
      con9 = 0.9d0
28
 
      con4 = 0.4d-01
29
 
      half = 0.5d0
30
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
31
 
c  part 1: determination of the number of knots and their position     c
32
 
c  **************************************************************      c
33
 
c  given a set of knots we compute the least-squares spline sinf(x),   c
34
 
c  and the corresponding sum of squared residuals fp=f(p=inf).         c
35
 
c  if iopt=-1 sinf(x) is the requested approximation.                  c
36
 
c  if iopt=0 or iopt=1 we check whether we can accept the knots:       c
37
 
c    if fp <=s we will continue with the current set of knots.         c
38
 
c    if fp > s we will increase the number of knots and compute the    c
39
 
c       corresponding least-squares spline until finally fp<=s.        c
40
 
c    the initial choice of knots depends on the value of s and iopt.   c
41
 
c    if s=0 we have spline interpolation; in that case the number of   c
42
 
c    knots equals nmax = m+k+1.                                        c
43
 
c    if s > 0 and                                                      c
44
 
c      iopt=0 we first compute the least-squares polynomial of         c
45
 
c      degree k; n = nmin = 2*k+2                                      c
46
 
c      iopt=1 we start with the set of knots found at the last         c
47
 
c      call of the routine, except for the case that s > fp0; then     c
48
 
c      we compute directly the least-squares polynomial of degree k.   c
49
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
50
 
c  determine nmin, the number of knots for polynomial approximation.
51
 
      nmin = 2*k1
52
 
      if(iopt.lt.0) go to 60
53
 
c  calculation of acc, the absolute tolerance for the root of f(p)=s.
54
 
      acc = tol*s
55
 
c  determine nmax, the number of knots for spline interpolation.
56
 
      nmax = m+k1
57
 
      if(s.gt.0.0d0) go to 45
58
 
c  if s=0, s(x) is an interpolating spline.
59
 
c  test whether the required storage space exceeds the available one.
60
 
      n = nmax
61
 
      if(nmax.gt.nest) go to 420
62
 
c  find the position of the interior knots in case of interpolation.
63
 
  10  mk1 = m-k1
64
 
      if(mk1.eq.0) go to 60
65
 
      k3 = k/2
66
 
      i = k2
67
 
      j = k3+2
68
 
      if(k3*2.eq.k) go to 30
69
 
      do 20 l=1,mk1
70
 
        t(i) = x(j)
71
 
        i = i+1
72
 
        j = j+1
73
 
  20  continue
74
 
      go to 60
75
 
  30  do 40 l=1,mk1
76
 
        t(i) = (x(j)+x(j-1))*half
77
 
        i = i+1
78
 
        j = j+1
79
 
  40  continue
80
 
      go to 60
81
 
c  if s>0 our initial choice of knots depends on the value of iopt.
82
 
c  if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
83
 
c  polynomial of degree k which is a spline without interior knots.
84
 
c  if iopt=1 and fp0>s we start computing the least squares spline
85
 
c  according to the set of knots found at the last call of the routine.
86
 
  45  if(iopt.eq.0) go to 50
87
 
      if(n.eq.nmin) go to 50
88
 
      fp0 = fpint(n)
89
 
      fpold = fpint(n-1)
90
 
      nplus = nrdata(n)
91
 
      if(fp0.gt.s) go to 60
92
 
  50  n = nmin
93
 
      fpold = 0.0d0
94
 
      nplus = 0
95
 
      nrdata(1) = m-2
96
 
c  main loop for the different sets of knots. m is a save upper bound
97
 
c  for the number of trials.
98
 
  60  do 200 iter = 1,m
99
 
        if(n.eq.nmin) ier = -2
100
 
c  find nrint, tne number of knot intervals.
101
 
        nrint = n-nmin+1
102
 
c  find the position of the additional knots which are needed for
103
 
c  the b-spline representation of s(x).
104
 
        nk1 = n-k1
105
 
        i = n
106
 
        do 70 j=1,k1
107
 
          t(j) = xb
108
 
          t(i) = xe
109
 
          i = i-1
110
 
  70    continue
111
 
c  compute the b-spline coefficients of the least-squares spline
112
 
c  sinf(x). the observation matrix a is built up row by row and
113
 
c  reduced to upper triangular form by givens transformations.
114
 
c  at the same time fp=f(p=inf) is computed.
115
 
        fp = 0.0d0
116
 
c  initialize the observation matrix a.
117
 
        do 80 i=1,nk1
118
 
          z(i) = 0.0d0
119
 
          do 80 j=1,k1
120
 
            a(i,j) = 0.0d0
121
 
  80    continue
122
 
        l = k1
123
 
        do 130 it=1,m
124
 
c  fetch the current data point x(it),y(it).
125
 
          xi = x(it)
126
 
          wi = w(it)
127
 
          yi = y(it)*wi
128
 
c  search for knot interval t(l) <= xi < t(l+1).
129
 
  85      if(xi.lt.t(l+1) .or. l.eq.nk1) go to 90
130
 
          l = l+1
131
 
          go to 85
132
 
c  evaluate the (k+1) non-zero b-splines at xi and store them in q.
133
 
  90      call fpbspl(t,n,k,xi,l,h)
134
 
          do 95 i=1,k1
135
 
            q(it,i) = h(i)
136
 
            h(i) = h(i)*wi
137
 
  95      continue
138
 
c  rotate the new row of the observation matrix into triangle.
139
 
          j = l-k1
140
 
          do 110 i=1,k1
141
 
            j = j+1
142
 
            piv = h(i)
143
 
            if(piv.eq.0.0d0) go to 110
144
 
c  calculate the parameters of the givens transformation.
145
 
            call fpgivs(piv,a(j,1),cos,sin)
146
 
c  transformations to right hand side.
147
 
            call fprota(cos,sin,yi,z(j))
148
 
            if(i.eq.k1) go to 120
149
 
            i2 = 1
150
 
            i3 = i+1
151
 
            do 100 i1 = i3,k1
152
 
              i2 = i2+1
153
 
c  transformations to left hand side.
154
 
              call fprota(cos,sin,h(i1),a(j,i2))
155
 
 100        continue
156
 
 110      continue
157
 
c  add contribution of this row to the sum of squares of residual
158
 
c  right hand sides.
159
 
 120      fp = fp+yi*yi
160
 
 130    continue
161
 
        if(ier.eq.(-2)) fp0 = fp
162
 
        fpint(n) = fp0
163
 
        fpint(n-1) = fpold
164
 
        nrdata(n) = nplus
165
 
c  backward substitution to obtain the b-spline coefficients.
166
 
        call fpback(a,z,nk1,k1,c,nest)
167
 
c  test whether the approximation sinf(x) is an acceptable solution.
168
 
        if(iopt.lt.0) go to 440
169
 
        fpms = fp-s
170
 
        if(abs(fpms).lt.acc) go to 440
171
 
c  if f(p=inf) < s accept the choice of knots.
172
 
        if(fpms.lt.0.0d0) go to 250
173
 
c  if n = nmax, sinf(x) is an interpolating spline.
174
 
        if(n.eq.nmax) go to 430
175
 
c  increase the number of knots.
176
 
c  if n=nest we cannot increase the number of knots because of
177
 
c  the storage capacity limitation.
178
 
        if(n.eq.nest) go to 420
179
 
c  determine the number of knots nplus we are going to add.
180
 
        if(ier.eq.0) go to 140
181
 
        nplus = 1
182
 
        ier = 0
183
 
        go to 150
184
 
 140    npl1 = nplus*2
185
 
        rn = nplus
186
 
        if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
187
 
        nplus = min0(nplus*2,max0(npl1,nplus/2,1))
188
 
 150    fpold = fp
189
 
c  compute the sum((w(i)*(y(i)-s(x(i))))**2) for each knot interval
190
 
c  t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
191
 
        fpart = 0.0d0
192
 
        i = 1
193
 
        l = k2
194
 
        new = 0
195
 
        do 180 it=1,m
196
 
          if(x(it).lt.t(l) .or. l.gt.nk1) go to 160
197
 
          new = 1
198
 
          l = l+1
199
 
 160      term = 0.0d0
200
 
          l0 = l-k2
201
 
          do 170 j=1,k1
202
 
            l0 = l0+1
203
 
            term = term+c(l0)*q(it,j)
204
 
 170      continue
205
 
          term = (w(it)*(term-y(it)))**2
206
 
          fpart = fpart+term
207
 
          if(new.eq.0) go to 180
208
 
          store = term*half
209
 
          fpint(i) = fpart-store
210
 
          i = i+1
211
 
          fpart = store
212
 
          new = 0
213
 
 180    continue
214
 
        fpint(nrint) = fpart
215
 
        do 190 l=1,nplus
216
 
c  add a new knot.
217
 
          call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
218
 
c  if n=nmax we locate the knots as for interpolation.
219
 
          if(n.eq.nmax) go to 10
220
 
c  test whether we cannot further increase the number of knots.
221
 
          if(n.eq.nest) go to 200
222
 
 190    continue
223
 
c  restart the computations with the new set of knots.
224
 
 200  continue
225
 
c  test whether the least-squares kth degree polynomial is a solution
226
 
c  of our approximation problem.
227
 
 250  if(ier.eq.(-2)) go to 440
228
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
229
 
c  part 2: determination of the smoothing spline sp(x).                c
230
 
c  ***************************************************                 c
231
 
c  we have determined the number of knots and their position.          c
232
 
c  we now compute the b-spline coefficients of the smoothing spline    c
233
 
c  sp(x). the observation matrix a is extended by the rows of matrix   c
234
 
c  b expressing that the kth derivative discontinuities of sp(x) at    c
235
 
c  the interior knots t(k+2),...t(n-k-1) must be zero. the corres-     c
236
 
c  ponding weights of these additional rows are set to 1/p.            c
237
 
c  iteratively we then have to determine the value of p such that      c
238
 
c  f(p)=sum((w(i)*(y(i)-sp(x(i))))**2) be = s. we already know that    c
239
 
c  the least-squares kth degree polynomial corresponds to p=0, and     c
240
 
c  that the least-squares spline corresponds to p=infinity. the        c
241
 
c  iteration process which is proposed here, makes use of rational     c
242
 
c  interpolation. since f(p) is a convex and strictly decreasing       c
243
 
c  function of p, it can be approximated by a rational function        c
244
 
c  r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond-  c
245
 
c  ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used      c
246
 
c  to calculate the new value of p such that r(p)=s. convergence is    c
247
 
c  guaranteed by taking f1>0 and f3<0.                                 c
248
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
249
 
c  evaluate the discontinuity jump of the kth derivative of the
250
 
c  b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
251
 
      call fpdisc(t,n,k2,b,nest)
252
 
c  initial value for p.
253
 
      p1 = 0.0d0
254
 
      f1 = fp0-s
255
 
      p3 = -one
256
 
      f3 = fpms
257
 
      p = 0.
258
 
      do 255 i=1,nk1
259
 
         p = p+a(i,1)
260
 
 255  continue
261
 
      rn = nk1
262
 
      p = rn/p
263
 
      ich1 = 0
264
 
      ich3 = 0
265
 
      n8 = n-nmin
266
 
c  iteration process to find the root of f(p) = s.
267
 
      do 360 iter=1,maxit
268
 
c  the rows of matrix b with weight 1/p are rotated into the
269
 
c  triangularised observation matrix a which is stored in g.
270
 
        pinv = one/p
271
 
        do 260 i=1,nk1
272
 
          c(i) = z(i)
273
 
          g(i,k2) = 0.0d0
274
 
          do 260 j=1,k1
275
 
            g(i,j) = a(i,j)
276
 
 260    continue
277
 
        do 300 it=1,n8
278
 
c  the row of matrix b is rotated into triangle by givens transformation
279
 
          do 270 i=1,k2
280
 
            h(i) = b(it,i)*pinv
281
 
 270      continue
282
 
          yi = 0.0d0
283
 
          do 290 j=it,nk1
284
 
            piv = h(1)
285
 
c  calculate the parameters of the givens transformation.
286
 
            call fpgivs(piv,g(j,1),cos,sin)
287
 
c  transformations to right hand side.
288
 
            call fprota(cos,sin,yi,c(j))
289
 
            if(j.eq.nk1) go to 300
290
 
            i2 = k1
291
 
            if(j.gt.n8) i2 = nk1-j
292
 
            do 280 i=1,i2
293
 
c  transformations to left hand side.
294
 
              i1 = i+1
295
 
              call fprota(cos,sin,h(i1),g(j,i1))
296
 
              h(i) = h(i1)
297
 
 280        continue
298
 
            h(i2+1) = 0.0d0
299
 
 290      continue
300
 
 300    continue
301
 
c  backward substitution to obtain the b-spline coefficients.
302
 
        call fpback(g,c,nk1,k2,c,nest)
303
 
c  computation of f(p).
304
 
        fp = 0.0d0
305
 
        l = k2
306
 
        do 330 it=1,m
307
 
          if(x(it).lt.t(l) .or. l.gt.nk1) go to 310
308
 
          l = l+1
309
 
 310      l0 = l-k2
310
 
          term = 0.0d0
311
 
          do 320 j=1,k1
312
 
            l0 = l0+1
313
 
            term = term+c(l0)*q(it,j)
314
 
 320      continue
315
 
          fp = fp+(w(it)*(term-y(it)))**2
316
 
 330    continue
317
 
c  test whether the approximation sp(x) is an acceptable solution.
318
 
        fpms = fp-s
319
 
        if(abs(fpms).lt.acc) go to 440
320
 
c  test whether the maximal number of iterations is reached.
321
 
        if(iter.eq.maxit) go to 400
322
 
c  carry out one more step of the iteration process.
323
 
        p2 = p
324
 
        f2 = fpms
325
 
        if(ich3.ne.0) go to 340
326
 
        if((f2-f3).gt.acc) go to 335
327
 
c  our initial choice of p is too large.
328
 
        p3 = p2
329
 
        f3 = f2
330
 
        p = p*con4
331
 
        if(p.le.p1) p=p1*con9 + p2*con1
332
 
        go to 360
333
 
 335    if(f2.lt.0.0d0) ich3=1
334
 
 340    if(ich1.ne.0) go to 350
335
 
        if((f1-f2).gt.acc) go to 345
336
 
c  our initial choice of p is too small
337
 
        p1 = p2
338
 
        f1 = f2
339
 
        p = p/con4
340
 
        if(p3.lt.0.) go to 360
341
 
        if(p.ge.p3) p = p2*con1 + p3*con9
342
 
        go to 360
343
 
 345    if(f2.gt.0.0d0) ich1=1
344
 
c  test whether the iteration process proceeds as theoretically
345
 
c  expected.
346
 
 350    if(f2.ge.f1 .or. f2.le.f3) go to 410
347
 
c  find the new value for p.
348
 
        p = fprati(p1,f1,p2,f2,p3,f3)
349
 
 360  continue
350
 
c  error codes and messages.
351
 
 400  ier = 3
352
 
      go to 440
353
 
 410  ier = 2
354
 
      go to 440
355
 
 420  ier = 1
356
 
      go to 440
357
 
 430  ier = -1
358
 
 440  return
359
 
      end