← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fppola.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/fppola.f
1
 
      subroutine fppola(iopt1,iopt2,iopt3,m,u,v,z,w,rad,s,nuest,nvest,
2
 
     * eta,tol,maxit,ib1,ib3,nc,ncc,intest,nrest,nu,tu,nv,tv,c,fp,sup,
3
 
     * fpint,coord,f,ff,row,cs,cosi,a,q,bu,bv,spu,spv,h,index,nummer,
4
 
     * wrk,lwrk,ier)
5
 
c  ..scalar arguments..
6
 
      integer iopt1,iopt2,iopt3,m,nuest,nvest,maxit,ib1,ib3,nc,ncc,
7
 
     * intest,nrest,nu,nv,lwrk,ier
8
 
      real*8 s,eta,tol,fp,sup
9
 
c  ..array arguments..
10
 
      integer index(nrest),nummer(m)
11
 
      real*8 u(m),v(m),z(m),w(m),tu(nuest),tv(nvest),c(nc),fpint(intest)
12
 
     *,
13
 
     * coord(intest),f(ncc),ff(nc),row(nvest),cs(nvest),cosi(5,nvest),
14
 
     * a(ncc,ib1),q(ncc,ib3),bu(nuest,5),bv(nvest,5),spu(m,4),spv(m,4),
15
 
     * h(ib3),wrk(lwrk)
16
 
c  ..user supplied function..
17
 
      real*8 rad
18
 
c  ..local scalars..
19
 
      real*8 acc,arg,co,c1,c2,c3,c4,dmax,eps,fac,fac1,fac2,fpmax,fpms,
20
 
     * f1,f2,f3,hui,huj,p,pi,pinv,piv,pi2,p1,p2,p3,r,ratio,si,sigma,
21
 
     * sq,store,uu,u2,u3,wi,zi,rn,one,two,three,con1,con4,con9,half,ten
22
 
      integer i,iband,iband3,iband4,ich1,ich3,ii,il,in,ipar,ipar1,irot,
23
 
     * iter,i1,i2,i3,j,jl,jrot,j1,j2,k,l,la,lf,lh,ll,lu,lv,lwest,l1,l2,
24
 
     * l3,l4,ncof,ncoff,nvv,nv4,nreg,nrint,nrr,nr1,nuu,nu4,num,num1,
25
 
     * numin,nvmin,rank,iband1
26
 
c  ..local arrays..
27
 
      real*8 hu(4),hv(4)
28
 
c  ..function references..
29
 
      real*8 abs,atan,cos,fprati,sin,sqrt
30
 
      integer min0
31
 
c  ..subroutine references..
32
 
c    fporde,fpbspl,fpback,fpgivs,fprota,fprank,fpdisc,fprppo
33
 
c  ..
34
 
c  set constants
35
 
      one = 1
36
 
      two = 2
37
 
      three = 3
38
 
      ten = 10
39
 
      half = 0.5e0
40
 
      con1 = 0.1e0
41
 
      con9 = 0.9e0
42
 
      con4 = 0.4e-01
43
 
      pi = atan(one)*4
44
 
      pi2 = pi+pi
45
 
      ipar = iopt2*(iopt2+3)/2
46
 
      ipar1 = ipar+1
47
 
      eps = sqrt(eta)
48
 
      if(iopt1.lt.0) go to 90
49
 
      numin = 9
50
 
      nvmin = 9+iopt2*(iopt2+1)
51
 
c  calculation of acc, the absolute tolerance for the root of f(p)=s.
52
 
      acc = tol*s
53
 
      if(iopt1.eq.0) go to 10
54
 
      if(s.lt.sup) then
55
 
        if (nv.lt.nvmin) go to 70
56
 
        go to 90
57
 
      endif
58
 
c  if iopt1 = 0 we begin by computing the weighted least-squares
59
 
c  polymomial of the form
60
 
c     s(u,v) = f(1)*(1-u**3)+f(2)*u**3+f(3)*(u**2-u**3)+f(4)*(u-u**3)
61
 
c  where f(4) = 0 if iopt2> 0 , f(3) = 0 if iopt2 > 1 and
62
 
c        f(2) = 0 if iopt3> 0.
63
 
c  the corresponding weighted sum of squared residuals gives the upper
64
 
c  bound sup for the smoothing factor s.
65
 
  10  sup = 0.
66
 
      do 20 i=1,4
67
 
         f(i) = 0.
68
 
         do 20 j=1,4
69
 
            a(i,j) = 0.
70
 
 20   continue
71
 
      do 50 i=1,m
72
 
         wi = w(i)
73
 
         zi = z(i)*wi
74
 
         uu = u(i)
75
 
         u2 = uu*uu
76
 
         u3 = uu*u2
77
 
         h(1) = (one-u3)*wi
78
 
         h(2) = u3*wi
79
 
         h(3) = u2*(one-uu)*wi
80
 
         h(4) = uu*(one-u2)*wi
81
 
         if(iopt3.ne.0) h(2) = 0.
82
 
         if(iopt2.gt.1) h(3) = 0.
83
 
         if(iopt2.gt.0) h(4) = 0.
84
 
         do 40 j=1,4
85
 
            piv = h(j)
86
 
            if(piv.eq.0.) go to 40
87
 
            call fpgivs(piv,a(j,1),co,si)
88
 
            call fprota(co,si,zi,f(j))
89
 
            if(j.eq.4) go to 40
90
 
            j1 = j+1
91
 
            j2 = 1
92
 
            do 30 l=j1,4
93
 
               j2 = j2+1
94
 
               call fprota(co,si,h(l),a(j,j2))
95
 
  30        continue
96
 
  40     continue
97
 
         sup = sup+zi*zi
98
 
  50  continue
99
 
      if(a(4,1).ne.0.) f(4) = f(4)/a(4,1)
100
 
      if(a(3,1).ne.0.) f(3) = (f(3)-a(3,2)*f(4))/a(3,1)
101
 
      if(a(2,1).ne.0.) f(2) = (f(2)-a(2,2)*f(3)-a(2,3)*f(4))/a(2,1)
102
 
      if(a(1,1).ne.0.)
103
 
     * f(1) = (f(1)-a(1,2)*f(2)-a(1,3)*f(3)-a(1,4)*f(4))/a(1,1)
104
 
c  find the b-spline representation of this least-squares polynomial
105
 
      c1 = f(1)
106
 
      c4 = f(2)
107
 
      c2 = f(4)/three+c1
108
 
      c3 = (f(3)+two*f(4))/three+c1
109
 
      nu = 8
110
 
      nv = 8
111
 
      do 60 i=1,4
112
 
         c(i) = c1
113
 
         c(i+4) = c2
114
 
         c(i+8) = c3
115
 
         c(i+12) = c4
116
 
         tu(i) = 0.
117
 
         tu(i+4) = one
118
 
         rn = 2*i-9
119
 
         tv(i) = rn*pi
120
 
         rn = 2*i-1
121
 
         tv(i+4) = rn*pi
122
 
  60  continue
123
 
      fp = sup
124
 
c  test whether the least-squares polynomial is an acceptable solution
125
 
      fpms = sup-s
126
 
      if(fpms.lt.acc) go to 960
127
 
c  test whether we cannot further increase the number of knots.
128
 
  70  if(nuest.lt.numin .or. nvest.lt.nvmin) go to 950
129
 
c  find the initial set of interior knots of the spline in case iopt1=0.
130
 
      nu = numin
131
 
      nv = nvmin
132
 
      tu(5) = half
133
 
      nvv = nv-8
134
 
      rn = nvv+1
135
 
      fac = pi2/rn
136
 
      do 80 i=1,nvv
137
 
         rn = i
138
 
         tv(i+4) = rn*fac-pi
139
 
  80  continue
140
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
141
 
c  part 1 : computation of least-squares bicubic splines.              c
142
 
c  ******************************************************              c
143
 
c  if iopt1<0 we compute the least-squares bicubic spline according    c
144
 
c  to the given set of knots.                                          c
145
 
c  if iopt1>=0 we compute least-squares bicubic splines with in-       c
146
 
c  creasing numbers of knots until the corresponding sum f(p=inf)<=s.  c
147
 
c  the initial set of knots then depends on the value of iopt1         c
148
 
c    if iopt1=0 we start with one interior knot in the u-direction     c
149
 
c              (0.5) and 1+iopt2*(iopt2+1) in the v-direction.         c
150
 
c    if iopt1>0 we start with the set of knots found at the last       c
151
 
c              call of the routine.                                    c
152
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
153
 
c  main loop for the different sets of knots. m is a save upper bound
154
 
c  for the number of trials.
155
 
  90  do 570 iter=1,m
156
 
c  find the position of the additional knots which are needed for the
157
 
c  b-spline representation of s(u,v).
158
 
         l1 = 4
159
 
         l2 = l1
160
 
         l3 = nv-3
161
 
         l4 = l3
162
 
         tv(l2) = -pi
163
 
         tv(l3) = pi
164
 
         do 120 i=1,3
165
 
            l1 = l1+1
166
 
            l2 = l2-1
167
 
            l3 = l3+1
168
 
            l4 = l4-1
169
 
            tv(l2) = tv(l4)-pi2
170
 
            tv(l3) = tv(l1)+pi2
171
 
 120     continue
172
 
        l = nu
173
 
        do 130 i=1,4
174
 
          tu(i) = 0.
175
 
          tu(l) = one
176
 
          l = l-1
177
 
 130    continue
178
 
c  find nrint, the total number of knot intervals and nreg, the number
179
 
c  of panels in which the approximation domain is subdivided by the
180
 
c  intersection of knots.
181
 
        nuu = nu-7
182
 
        nvv = nv-7
183
 
        nrr = nvv/2
184
 
        nr1 = nrr+1
185
 
        nrint = nuu+nvv
186
 
        nreg = nuu*nvv
187
 
c  arrange the data points according to the panel they belong to.
188
 
        call fporde(u,v,m,3,3,tu,nu,tv,nv,nummer,index,nreg)
189
 
        if(iopt2.eq.0) go to 195
190
 
c  find the b-spline coefficients cosi of the cubic spline
191
 
c  approximations for cr(v)=rad(v)*cos(v) and sr(v) = rad(v)*sin(v)
192
 
c  if iopt2=1, and additionally also for cr(v)**2,sr(v)**2 and
193
 
c  2*cr(v)*sr(v) if iopt2=2
194
 
        do 140 i=1,nvv
195
 
           do 135 j=1,ipar
196
 
              cosi(j,i) = 0.
197
 
 135       continue
198
 
           do 140 j=1,nvv
199
 
              a(i,j) = 0.
200
 
 140    continue
201
 
c  the coefficients cosi are obtained from interpolation conditions
202
 
c  at the knots tv(i),i=4,5,...nv-4.
203
 
        do 175 i=1,nvv
204
 
           l2 = i+3
205
 
           arg = tv(l2)
206
 
           call fpbspl(tv,nv,3,arg,l2,hv)
207
 
           do 145 j=1,nvv
208
 
              row(j) = 0.
209
 
 145       continue
210
 
           ll = i
211
 
           do 150 j=1,3
212
 
              if(ll.gt.nvv) ll= 1
213
 
              row(ll) = row(ll)+hv(j)
214
 
              ll = ll+1
215
 
 150       continue
216
 
           co = cos(arg)
217
 
           si = sin(arg)
218
 
           r = rad(arg)
219
 
           cs(1) = co*r
220
 
           cs(2) = si*r
221
 
           if(iopt2.eq.1) go to 155
222
 
           cs(3) = cs(1)*cs(1)
223
 
           cs(4) = cs(2)*cs(2)
224
 
           cs(5) = cs(1)*cs(2)
225
 
 155       do 170 j=1,nvv
226
 
              piv = row(j)
227
 
              if(piv.eq.0.) go to 170
228
 
              call fpgivs(piv,a(j,1),co,si)
229
 
              do 160 l=1,ipar
230
 
                 call fprota(co,si,cs(l),cosi(l,j))
231
 
 160          continue
232
 
              if(j.eq.nvv) go to 175
233
 
              j1 = j+1
234
 
              j2 = 1
235
 
              do 165 l=j1,nvv
236
 
                 j2 = j2+1
237
 
                 call fprota(co,si,row(l),a(j,j2))
238
 
 165          continue
239
 
 170       continue
240
 
 175    continue
241
 
         do 190 l=1,ipar
242
 
            do 180 j=1,nvv
243
 
               cs(j) = cosi(l,j)
244
 
 180        continue
245
 
            call fpback(a,cs,nvv,nvv,cs,ncc)
246
 
            do 185 j=1,nvv
247
 
               cosi(l,j) = cs(j)
248
 
 185        continue
249
 
 190     continue
250
 
c  find ncof, the dimension of the spline and ncoff, the number
251
 
c  of coefficients in the standard b-spline representation.
252
 
 195    nu4 = nu-4
253
 
        nv4 = nv-4
254
 
        ncoff = nu4*nv4
255
 
        ncof = ipar1+nvv*(nu4-1-iopt2-iopt3)
256
 
c  find the bandwidth of the observation matrix a.
257
 
        iband = 4*nvv
258
 
        if(nuu-iopt2-iopt3.le.1) iband = ncof
259
 
        iband1 = iband-1
260
 
c  initialize the observation matrix a.
261
 
        do 200 i=1,ncof
262
 
          f(i) = 0.
263
 
          do 200 j=1,iband
264
 
            a(i,j) = 0.
265
 
 200    continue
266
 
c  initialize the sum of squared residuals.
267
 
        fp = 0.
268
 
        ratio = one+tu(6)/tu(5)
269
 
c  fetch the data points in the new order. main loop for the
270
 
c  different panels.
271
 
        do 380 num=1,nreg
272
 
c  fix certain constants for the current panel; jrot records the column
273
 
c  number of the first non-zero element in a row of the observation
274
 
c  matrix according to a data point of the panel.
275
 
          num1 = num-1
276
 
          lu = num1/nvv
277
 
          l1 = lu+4
278
 
          lv = num1-lu*nvv+1
279
 
          l2 = lv+3
280
 
          jrot = 0
281
 
          if(lu.gt.iopt2) jrot = ipar1+(lu-iopt2-1)*nvv
282
 
          lu = lu+1
283
 
c  test whether there are still data points in the current panel.
284
 
          in = index(num)
285
 
 210      if(in.eq.0) go to 380
286
 
c  fetch a new data point.
287
 
          wi = w(in)
288
 
          zi = z(in)*wi
289
 
c  evaluate for the u-direction, the 4 non-zero b-splines at u(in)
290
 
          call fpbspl(tu,nu,3,u(in),l1,hu)
291
 
c  evaluate for the v-direction, the 4 non-zero b-splines at v(in)
292
 
          call fpbspl(tv,nv,3,v(in),l2,hv)
293
 
c  store the value of these b-splines in spu and spv resp.
294
 
          do 220 i=1,4
295
 
            spu(in,i) = hu(i)
296
 
            spv(in,i) = hv(i)
297
 
 220      continue
298
 
c  initialize the new row of observation matrix.
299
 
          do 240 i=1,iband
300
 
            h(i) = 0.
301
 
 240      continue
302
 
c  calculate the non-zero elements of the new row by making the cross
303
 
c  products of the non-zero b-splines in u- and v-direction and
304
 
c  by taking into account the conditions of the splines.
305
 
          do 250 i=1,nvv
306
 
             row(i) = 0.
307
 
 250      continue
308
 
c  take into account the periodicity condition of the bicubic splines.
309
 
          ll = lv
310
 
          do 260 i=1,4
311
 
             if(ll.gt.nvv) ll=1
312
 
             row(ll) = row(ll)+hv(i)
313
 
             ll = ll+1
314
 
 260      continue
315
 
c  take into account the other conditions of the splines.
316
 
          if(iopt2.eq.0 .or. lu.gt.iopt2+1) go to 280
317
 
          do 270 l=1,ipar
318
 
             cs(l) = 0.
319
 
             do 270 i=1,nvv
320
 
                cs(l) = cs(l)+row(i)*cosi(l,i)
321
 
 270     continue
322
 
c  fill in the non-zero elements of the new row.
323
 
 280     j1 = 0
324
 
         do 330 j =1,4
325
 
            jlu = j+lu
326
 
            huj = hu(j)
327
 
            if(jlu.gt.iopt2+2) go to 320
328
 
            go to (290,290,300,310),jlu
329
 
 290        h(1) = huj
330
 
            j1 = 1
331
 
            go to 330
332
 
 300        h(1) = h(1)+huj
333
 
            h(2) = huj*cs(1)
334
 
            h(3) = huj*cs(2)
335
 
            j1 = 3
336
 
            go to 330
337
 
 310        h(1) = h(1)+huj
338
 
            h(2) = h(2)+huj*ratio*cs(1)
339
 
            h(3) = h(3)+huj*ratio*cs(2)
340
 
            h(4) = huj*cs(3)
341
 
            h(5) = huj*cs(4)
342
 
            h(6) = huj*cs(5)
343
 
            j1 = 6
344
 
            go to 330
345
 
 320        if(jlu.gt.nu4 .and. iopt3.ne.0) go to 330
346
 
            do 325 i=1,nvv
347
 
               j1 = j1+1
348
 
               h(j1) = row(i)*huj
349
 
 325        continue
350
 
 330      continue
351
 
          do 335 i=1,iband
352
 
            h(i) = h(i)*wi
353
 
 335      continue
354
 
c  rotate the row into triangle by givens transformations.
355
 
          irot = jrot
356
 
          do 350 i=1,iband
357
 
            irot = irot+1
358
 
            piv = h(i)
359
 
            if(piv.eq.0.) go to 350
360
 
c  calculate the parameters of the givens transformation.
361
 
            call fpgivs(piv,a(irot,1),co,si)
362
 
c  apply that transformation to the right hand side.
363
 
            call fprota(co,si,zi,f(irot))
364
 
            if(i.eq.iband) go to 360
365
 
c  apply that transformation to the left hand side.
366
 
            i2 = 1
367
 
            i3 = i+1
368
 
            do 340 j=i3,iband
369
 
              i2 = i2+1
370
 
              call fprota(co,si,h(j),a(irot,i2))
371
 
 340        continue
372
 
 350      continue
373
 
c  add the contribution of the row to the sum of squares of residual
374
 
c  right hand sides.
375
 
 360      fp = fp+zi**2
376
 
c  find the number of the next data point in the panel.
377
 
 370      in = nummer(in)
378
 
          go to 210
379
 
 380    continue
380
 
c  find dmax, the maximum value for the diagonal elements in the reduced
381
 
c  triangle.
382
 
        dmax = 0.
383
 
        do 390 i=1,ncof
384
 
          if(a(i,1).le.dmax) go to 390
385
 
          dmax = a(i,1)
386
 
 390    continue
387
 
c  check whether the observation matrix is rank deficient.
388
 
        sigma = eps*dmax
389
 
        do 400 i=1,ncof
390
 
          if(a(i,1).le.sigma) go to 410
391
 
 400    continue
392
 
c  backward substitution in case of full rank.
393
 
        call fpback(a,f,ncof,iband,c,ncc)
394
 
        rank = ncof
395
 
        do 405 i=1,ncof
396
 
          q(i,1) = a(i,1)/dmax
397
 
 405    continue
398
 
        go to 430
399
 
c  in case of rank deficiency, find the minimum norm solution.
400
 
 410    lwest = ncof*iband+ncof+iband
401
 
        if(lwrk.lt.lwest) go to 925
402
 
        lf = 1
403
 
        lh = lf+ncof
404
 
        la = lh+iband
405
 
        do 420 i=1,ncof
406
 
          ff(i) = f(i)
407
 
          do 420 j=1,iband
408
 
            q(i,j) = a(i,j)
409
 
 420    continue
410
 
        call fprank(q,ff,ncof,iband,ncc,sigma,c,sq,rank,wrk(la),
411
 
     *   wrk(lf),wrk(lh))
412
 
        do 425 i=1,ncof
413
 
          q(i,1) = q(i,1)/dmax
414
 
 425    continue
415
 
c  add to the sum of squared residuals, the contribution of reducing
416
 
c  the rank.
417
 
        fp = fp+sq
418
 
c  find the coefficients in the standard b-spline representation of
419
 
c  the spline.
420
 
 430    call fprppo(nu,nv,iopt2,iopt3,cosi,ratio,c,ff,ncoff)
421
 
c  test whether the least-squares spline is an acceptable solution.
422
 
        if(iopt1.lt.0) then
423
 
          if (fp.le.0) go to 970
424
 
          go to 980
425
 
        endif
426
 
        fpms = fp-s
427
 
        if(abs(fpms).le.acc) then
428
 
            if (fp.le.0) go to 970
429
 
            go to 980
430
 
        endif
431
 
c  if f(p=inf) < s, accept the choice of knots.
432
 
        if(fpms.lt.0.) go to 580
433
 
c  test whether we cannot further increase the number of knots
434
 
        if(m.lt.ncof) go to 935
435
 
c  search where to add a new knot.
436
 
c  find for each interval the sum of squared residuals fpint for the
437
 
c  data points having the coordinate belonging to that knot interval.
438
 
c  calculate also coord which is the same sum, weighted by the position
439
 
c  of the data points considered.
440
 
 440    do 450 i=1,nrint
441
 
          fpint(i) = 0.
442
 
          coord(i) = 0.
443
 
 450    continue
444
 
        do 490 num=1,nreg
445
 
          num1 = num-1
446
 
          lu = num1/nvv
447
 
          l1 = lu+1
448
 
          lv = num1-lu*nvv
449
 
          l2 = lv+1+nuu
450
 
          jrot = lu*nv4+lv
451
 
          in = index(num)
452
 
 460      if(in.eq.0) go to 490
453
 
          store = 0.
454
 
          i1 = jrot
455
 
          do 480 i=1,4
456
 
            hui = spu(in,i)
457
 
            j1 = i1
458
 
            do 470 j=1,4
459
 
              j1 = j1+1
460
 
              store = store+hui*spv(in,j)*c(j1)
461
 
 470        continue
462
 
            i1 = i1+nv4
463
 
 480      continue
464
 
          store = (w(in)*(z(in)-store))**2
465
 
          fpint(l1) = fpint(l1)+store
466
 
          coord(l1) = coord(l1)+store*u(in)
467
 
          fpint(l2) = fpint(l2)+store
468
 
          coord(l2) = coord(l2)+store*v(in)
469
 
          in = nummer(in)
470
 
          go to 460
471
 
 490    continue
472
 
c bring together the information concerning knot panels which are
473
 
c symmetric with respect to the origin.
474
 
        do 495 i=1,nrr
475
 
          l1 = nuu+i
476
 
          l2 = l1+nrr
477
 
          fpint(l1) = fpint(l1)+fpint(l2)
478
 
          coord(l1) = coord(l1)+coord(l2)-pi*fpint(l2)
479
 
 495    continue
480
 
c  find the interval for which fpint is maximal on the condition that
481
 
c  there still can be added a knot.
482
 
        l1 = 1
483
 
        l2 = nuu+nrr
484
 
        if(nuest.lt.nu+1) l1=nuu+1
485
 
        if(nvest.lt.nv+2) l2=nuu
486
 
c  test whether we cannot further increase the number of knots.
487
 
        if(l1.gt.l2) go to 950
488
 
 500    fpmax = 0.
489
 
        l = 0
490
 
        do 510 i=l1,l2
491
 
          if(fpmax.ge.fpint(i)) go to 510
492
 
          l = i
493
 
          fpmax = fpint(i)
494
 
 510    continue
495
 
        if(l.eq.0) go to 930
496
 
c  calculate the position of the new knot.
497
 
        arg = coord(l)/fpint(l)
498
 
c  test in what direction the new knot is going to be added.
499
 
        if(l.gt.nuu) go to 530
500
 
c  addition in the u-direction
501
 
        l4 = l+4
502
 
        fpint(l) = 0.
503
 
        fac1 = tu(l4)-arg
504
 
        fac2 = arg-tu(l4-1)
505
 
        if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 500
506
 
        j = nu
507
 
        do 520 i=l4,nu
508
 
          tu(j+1) = tu(j)
509
 
          j = j-1
510
 
 520    continue
511
 
        tu(l4) = arg
512
 
        nu = nu+1
513
 
        go to 570
514
 
c  addition in the v-direction
515
 
 530    l4 = l+4-nuu
516
 
        fpint(l) = 0.
517
 
        fac1 = tv(l4)-arg
518
 
        fac2 = arg-tv(l4-1)
519
 
        if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 500
520
 
        ll = nrr+4
521
 
        j = ll
522
 
        do 550 i=l4,ll
523
 
          tv(j+1) = tv(j)
524
 
          j = j-1
525
 
 550    continue
526
 
        tv(l4) = arg
527
 
        nv = nv+2
528
 
        nrr = nrr+1
529
 
        do 560 i=5,ll
530
 
          j = i+nrr
531
 
          tv(j) = tv(i)+pi
532
 
 560    continue
533
 
c  restart the computations with the new set of knots.
534
 
 570  continue
535
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
536
 
c part 2: determination of the smoothing bicubic spline.               c
537
 
c ******************************************************               c
538
 
c we have determined the number of knots and their position. we now    c
539
 
c compute the coefficients of the smoothing spline sp(u,v).            c
540
 
c the observation matrix a is extended by the rows of a matrix, expres-c
541
 
c sing that sp(u,v) must be a constant function in the variable        c
542
 
c v and a cubic polynomial in the variable u. the corresponding        c
543
 
c weights of these additional rows are set to 1/(p). iteratively       c
544
 
c we than have to determine the value of p such that f(p) = sum((w(i)* c
545
 
c (z(i)-sp(u(i),v(i))))**2)  be = s.                                   c
546
 
c we already know that the least-squares polynomial corresponds to p=0,c
547
 
c and that the least-squares bicubic spline corresponds to p=infin.    c
548
 
c the iteration process makes use of rational interpolation. since f(p)c
549
 
c is a convex and strictly decreasing function of p, it can be approx- c
550
 
c imated by a rational function of the form r(p) = (u*p+v)/(p+w).      c
551
 
c three values of p (p1,p2,p3) with corresponding values of f(p) (f1=  c
552
 
c f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the new value   c
553
 
c of p such that r(p)=s. convergence is guaranteed by taking f1>0,f3<0.c
554
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
555
 
c  evaluate the discontinuity jumps of the 3-th order derivative of
556
 
c  the b-splines at the knots tu(l),l=5,...,nu-4.
557
 
 580  call fpdisc(tu,nu,5,bu,nuest)
558
 
c  evaluate the discontinuity jumps of the 3-th order derivative of
559
 
c  the b-splines at the knots tv(l),l=5,...,nv-4.
560
 
      call fpdisc(tv,nv,5,bv,nvest)
561
 
c  initial value for p.
562
 
      p1 = 0.
563
 
      f1 = sup-s
564
 
      p3 = -one
565
 
      f3 = fpms
566
 
      p = 0.
567
 
      do 590 i=1,ncof
568
 
        p = p+a(i,1)
569
 
 590  continue
570
 
      rn = ncof
571
 
      p = rn/p
572
 
c  find the bandwidth of the extended observation matrix.
573
 
      iband4 = iband+ipar1
574
 
      if(iband4.gt.ncof) iband4 = ncof
575
 
      iband3 = iband4 -1
576
 
      ich1 = 0
577
 
      ich3 = 0
578
 
      nuu = nu4-iopt3-1
579
 
c  iteration process to find the root of f(p)=s.
580
 
      do 920 iter=1,maxit
581
 
        pinv = one/p
582
 
c  store the triangularized observation matrix into q.
583
 
        do 630 i=1,ncof
584
 
          ff(i) = f(i)
585
 
          do 620 j=1,iband4
586
 
            q(i,j) = 0.
587
 
 620      continue
588
 
          do 630 j=1,iband
589
 
            q(i,j) = a(i,j)
590
 
 630    continue
591
 
c  extend the observation matrix with the rows of a matrix, expressing
592
 
c  that for u=constant sp(u,v) must be a constant function.
593
 
        do 720 i=5,nv4
594
 
          ii = i-4
595
 
          do 635 l=1,nvv
596
 
             row(l) = 0.
597
 
 635      continue
598
 
          ll = ii
599
 
          do 640  l=1,5
600
 
             if(ll.gt.nvv) ll=1
601
 
             row(ll) = row(ll)+bv(ii,l)
602
 
             ll = ll+1
603
 
 640      continue
604
 
          do 720 j=1,nuu
605
 
c  initialize the new row.
606
 
            do 645 l=1,iband
607
 
              h(l) = 0.
608
 
 645        continue
609
 
c  fill in the non-zero elements of the row. jrot records the column
610
 
c  number of the first non-zero element in the row.
611
 
            if(j.gt.iopt2) go to 665
612
 
            if(j.eq.2) go to 655
613
 
            do 650 k=1,2
614
 
               cs(k) = 0.
615
 
               do 650 l=1,nvv
616
 
                  cs(k) = cs(k)+cosi(k,l)*row(l)
617
 
 650        continue
618
 
            h(1) = cs(1)
619
 
            h(2) = cs(2)
620
 
            jrot = 2
621
 
            go to 675
622
 
 655        do 660 k=3,5
623
 
               cs(k) = 0.
624
 
               do 660 l=1,nvv
625
 
                  cs(k) = cs(k)+cosi(k,l)*row(l)
626
 
 660        continue
627
 
            h(1) = cs(1)*ratio
628
 
            h(2) = cs(2)*ratio
629
 
            h(3) = cs(3)
630
 
            h(4) = cs(4)
631
 
            h(5) = cs(5)
632
 
            jrot = 2
633
 
            go to 675
634
 
 665        do 670 l=1,nvv
635
 
               h(l) = row(l)
636
 
 670        continue
637
 
            jrot = ipar1+1+(j-iopt2-1)*nvv
638
 
 675        do 677 l=1,iband
639
 
              h(l) = h(l)*pinv
640
 
 677        continue
641
 
            zi = 0.
642
 
c  rotate the new row into triangle by givens transformations.
643
 
            do 710 irot=jrot,ncof
644
 
              piv = h(1)
645
 
              i2 = min0(iband1,ncof-irot)
646
 
              if(piv.eq.0.) then
647
 
                 if (i2.le.0) go to 720
648
 
                 go to 690
649
 
              endif
650
 
c  calculate the parameters of the givens transformation.
651
 
              call fpgivs(piv,q(irot,1),co,si)
652
 
c  apply that givens transformation to the right hand side.
653
 
              call fprota(co,si,zi,ff(irot))
654
 
              if(i2.eq.0) go to 720
655
 
c  apply that givens transformation to the left hand side.
656
 
              do 680 l=1,i2
657
 
                l1 = l+1
658
 
                call fprota(co,si,h(l1),q(irot,l1))
659
 
 680          continue
660
 
 690          do 700 l=1,i2
661
 
                h(l) = h(l+1)
662
 
 700          continue
663
 
              h(i2+1) = 0.
664
 
 710        continue
665
 
 720    continue
666
 
c  extend the observation matrix with the rows of a matrix expressing
667
 
c  that for v=constant. sp(u,v) must be a cubic polynomial.
668
 
        do 810 i=5,nu4
669
 
          ii = i-4
670
 
          do 810 j=1,nvv
671
 
c  initialize the new row
672
 
            do 730 l=1,iband4
673
 
              h(l) = 0.
674
 
 730        continue
675
 
c  fill in the non-zero elements of the row. jrot records the column
676
 
c  number of the first non-zero element in the row.
677
 
            j1 = 1
678
 
            do 760 l=1,5
679
 
               il = ii+l-1
680
 
               if(il.eq.nu4 .and. iopt3.ne.0) go to 760
681
 
               if(il.gt.iopt2+1) go to 750
682
 
               go to (735,740,745),il
683
 
 735           h(1) = bu(ii,l)
684
 
               j1 = j+1
685
 
               go to 760
686
 
 740           h(1) = h(1)+bu(ii,l)
687
 
               h(2) = bu(ii,l)*cosi(1,j)
688
 
               h(3) = bu(ii,l)*cosi(2,j)
689
 
               j1 = j+3
690
 
               go to 760
691
 
 745           h(1) = h(1)+bu(ii,l)
692
 
               h(2) = bu(ii,l)*cosi(1,j)*ratio
693
 
               h(3) = bu(ii,l)*cosi(2,j)*ratio
694
 
               h(4) = bu(ii,l)*cosi(3,j)
695
 
               h(5) = bu(ii,l)*cosi(4,j)
696
 
               h(6) = bu(ii,l)*cosi(5,j)
697
 
               j1 = j+6
698
 
               go to 760
699
 
 750           h(j1) = bu(ii,l)
700
 
               j1 = j1+nvv
701
 
 760        continue
702
 
            do 765 l=1,iband4
703
 
              h(l) = h(l)*pinv
704
 
 765        continue
705
 
            zi = 0.
706
 
            jrot = 1
707
 
            if(ii.gt.iopt2+1) jrot = ipar1+(ii-iopt2-2)*nvv+j
708
 
c  rotate the new row into triangle by givens transformations.
709
 
            do 800 irot=jrot,ncof
710
 
              piv = h(1)
711
 
              i2 = min0(iband3,ncof-irot)
712
 
              if(piv.eq.0.) then
713
 
                if (i2.le.0) go to 810
714
 
                go to 780
715
 
              endif
716
 
c  calculate the parameters of the givens transformation.
717
 
              call fpgivs(piv,q(irot,1),co,si)
718
 
c  apply that givens transformation to the right hand side.
719
 
              call fprota(co,si,zi,ff(irot))
720
 
              if(i2.eq.0) go to 810
721
 
c  apply that givens transformation to the left hand side.
722
 
              do 770 l=1,i2
723
 
                l1 = l+1
724
 
                call fprota(co,si,h(l1),q(irot,l1))
725
 
 770          continue
726
 
 780          do 790 l=1,i2
727
 
                h(l) = h(l+1)
728
 
 790          continue
729
 
              h(i2+1) = 0.
730
 
 800        continue
731
 
 810    continue
732
 
c  find dmax, the maximum value for the diagonal elements in the
733
 
c  reduced triangle.
734
 
        dmax = 0.
735
 
        do 820 i=1,ncof
736
 
          if(q(i,1).le.dmax) go to 820
737
 
          dmax = q(i,1)
738
 
 820    continue
739
 
c  check whether the matrix is rank deficient.
740
 
        sigma = eps*dmax
741
 
        do 830 i=1,ncof
742
 
          if(q(i,1).le.sigma) go to 840
743
 
 830    continue
744
 
c  backward substitution in case of full rank.
745
 
        call fpback(q,ff,ncof,iband4,c,ncc)
746
 
        rank = ncof
747
 
        go to 845
748
 
c  in case of rank deficiency, find the minimum norm solution.
749
 
 840    lwest = ncof*iband4+ncof+iband4
750
 
        if(lwrk.lt.lwest) go to 925
751
 
        lf = 1
752
 
        lh = lf+ncof
753
 
        la = lh+iband4
754
 
        call fprank(q,ff,ncof,iband4,ncc,sigma,c,sq,rank,wrk(la),
755
 
     *   wrk(lf),wrk(lh))
756
 
 845    do 850 i=1,ncof
757
 
           q(i,1) = q(i,1)/dmax
758
 
 850    continue
759
 
c  find the coefficients in the standard b-spline representation of
760
 
c  the polar spline.
761
 
        call fprppo(nu,nv,iopt2,iopt3,cosi,ratio,c,ff,ncoff)
762
 
c  compute f(p).
763
 
        fp = 0.
764
 
        do 890 num = 1,nreg
765
 
          num1 = num-1
766
 
          lu = num1/nvv
767
 
          lv = num1-lu*nvv
768
 
          jrot = lu*nv4+lv
769
 
          in = index(num)
770
 
 860      if(in.eq.0) go to 890
771
 
          store = 0.
772
 
          i1 = jrot
773
 
          do 880 i=1,4
774
 
            hui = spu(in,i)
775
 
            j1 = i1
776
 
            do 870 j=1,4
777
 
              j1 = j1+1
778
 
              store = store+hui*spv(in,j)*c(j1)
779
 
 870        continue
780
 
            i1 = i1+nv4
781
 
 880      continue
782
 
          fp = fp+(w(in)*(z(in)-store))**2
783
 
          in = nummer(in)
784
 
          go to 860
785
 
 890    continue
786
 
c  test whether the approximation sp(u,v) is an acceptable solution
787
 
        fpms = fp-s
788
 
        if(abs(fpms).le.acc) go to 980
789
 
c  test whether the maximum allowable number of iterations has been
790
 
c  reached.
791
 
        if(iter.eq.maxit) go to 940
792
 
c  carry out one more step of the iteration process.
793
 
        p2 = p
794
 
        f2 = fpms
795
 
        if(ich3.ne.0) go to 900
796
 
        if((f2-f3).gt.acc) go to 895
797
 
c  our initial choice of p is too large.
798
 
        p3 = p2
799
 
        f3 = f2
800
 
        p = p*con4
801
 
        if(p.le.p1) p = p1*con9 + p2*con1
802
 
        go to 920
803
 
 895    if(f2.lt.0.) ich3 = 1
804
 
 900    if(ich1.ne.0) go to 910
805
 
        if((f1-f2).gt.acc) go to 905
806
 
c  our initial choice of p is too small
807
 
        p1 = p2
808
 
        f1 = f2
809
 
        p = p/con4
810
 
        if(p3.lt.0.) go to 920
811
 
        if(p.ge.p3) p = p2*con1 +p3*con9
812
 
        go to 920
813
 
 905    if(f2.gt.0.) ich1 = 1
814
 
c  test whether the iteration process proceeds as theoretically
815
 
c  expected.
816
 
 910    if(f2.ge.f1 .or. f2.le.f3) go to 945
817
 
c  find the new value of p.
818
 
        p = fprati(p1,f1,p2,f2,p3,f3)
819
 
 920  continue
820
 
c  error codes and messages.
821
 
 925  ier = lwest
822
 
      go to 990
823
 
 930  ier = 5
824
 
      go to 990
825
 
 935  ier = 4
826
 
      go to 990
827
 
 940  ier = 3
828
 
      go to 990
829
 
 945  ier = 2
830
 
      go to 990
831
 
 950  ier = 1
832
 
      go to 990
833
 
 960  ier = -2
834
 
      go to 990
835
 
 970  ier = -1
836
 
      fp = 0.
837
 
 980  if(ncof.ne.rank) ier = -rank
838
 
 990  return
839
 
      end
840