← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpclos.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/fpclos.f
1
 
      subroutine fpclos(iopt,idim,m,u,mx,x,w,k,s,nest,tol,maxit,k1,k2,
2
 
     * n,t,nc,c,fp,fpint,z,a1,a2,b,g1,g2,q,nrdata,ier)
3
 
c  ..
4
 
c  ..scalar arguments..
5
 
      real*8 s,tol,fp
6
 
      integer iopt,idim,m,mx,k,nest,maxit,k1,k2,n,nc,ier
7
 
c  ..array arguments..
8
 
      real*8 u(m),x(mx),w(m),t(nest),c(nc),fpint(nest),z(nc),a1(nest,k1)
9
 
     *,
10
 
     * a2(nest,k),b(nest,k2),g1(nest,k2),g2(nest,k1),q(m,k1)
11
 
      integer nrdata(nest)
12
 
c  ..local scalars..
13
 
      real*8 acc,cos,d1,fac,fpart,fpms,fpold,fp0,f1,f2,f3,p,per,pinv,piv
14
 
     *,
15
 
     * p1,p2,p3,sin,store,term,ui,wi,rn,one,con1,con4,con9,half
16
 
      integer i,ich1,ich3,ij,ik,it,iter,i1,i2,i3,j,jj,jk,jper,j1,j2,kk,
17
 
     * kk1,k3,l,l0,l1,l5,mm,m1,new,nk1,nk2,nmax,nmin,nplus,npl1,
18
 
     * nrint,n10,n11,n7,n8
19
 
c  ..local arrays..
20
 
      real*8 h(6),h1(7),h2(6),xi(10)
21
 
c  ..function references..
22
 
      real*8 abs,fprati
23
 
      integer max0,min0
24
 
c  ..subroutine references..
25
 
c    fpbacp,fpbspl,fpgivs,fpdisc,fpknot,fprota
26
 
c  ..
27
 
c  set constants
28
 
      one = 0.1e+01
29
 
      con1 = 0.1e0
30
 
      con9 = 0.9e0
31
 
      con4 = 0.4e-01
32
 
      half = 0.5e0
33
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
34
 
c  part 1: determination of the number of knots and their position     c
35
 
c  **************************************************************      c
36
 
c  given a set of knots we compute the least-squares closed curve      c
37
 
c  sinf(u). if the sum f(p=inf) <= s we accept the choice of knots.    c
38
 
c  if iopt=-1 sinf(u) is the requested curve                           c
39
 
c  if iopt=0 or iopt=1 we check whether we can accept the knots:       c
40
 
c    if fp <=s we will continue with the current set of knots.         c
41
 
c    if fp > s we will increase the number of knots and compute the    c
42
 
c       corresponding least-squares curve until finally fp<=s.         c
43
 
c  the initial choice of knots depends on the value of s and iopt.     c
44
 
c    if s=0 we have spline interpolation; in that case the number of   c
45
 
c    knots equals nmax = m+2*k.                                        c
46
 
c    if s > 0 and                                                      c
47
 
c      iopt=0 we first compute the least-squares polynomial curve of   c
48
 
c      degree k; n = nmin = 2*k+2. since s(u) must be periodic we      c
49
 
c      find that s(u) reduces to a fixed point.                        c
50
 
c      iopt=1 we start with the set of knots found at the last         c
51
 
c      call of the routine, except for the case that s > fp0; then     c
52
 
c      we compute directly the least-squares polynomial curve.         c
53
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
54
 
      m1 = m-1
55
 
      kk = k
56
 
      kk1 = k1
57
 
      k3 = 3*k+1
58
 
      nmin = 2*k1
59
 
c  determine the length of the period of the splines.
60
 
      per = u(m)-u(1)
61
 
      if(iopt.lt.0) go to 50
62
 
c  calculation of acc, the absolute tolerance for the root of f(p)=s.
63
 
      acc = tol*s
64
 
c  determine nmax, the number of knots for periodic spline interpolation
65
 
      nmax = m+2*k
66
 
      if(s.gt.0. .or. nmax.eq.nmin) go to 30
67
 
c  if s=0, s(u) is an interpolating curve.
68
 
      n = nmax
69
 
c  test whether the required storage space exceeds the available one.
70
 
      if(n.gt.nest) go to 620
71
 
c  find the position of the interior knots in case of interpolation.
72
 
   5  if((k/2)*2 .eq.k) go to 20
73
 
      do 10 i=2,m1
74
 
        j = i+k
75
 
        t(j) = u(i)
76
 
  10  continue
77
 
      if(s.gt.0.) go to 50
78
 
      kk = k-1
79
 
      kk1 = k
80
 
      if(kk.gt.0) go to 50
81
 
      t(1) = t(m)-per
82
 
      t(2) = u(1)
83
 
      t(m+1) = u(m)
84
 
      t(m+2) = t(3)+per
85
 
      jj = 0
86
 
      do 15 i=1,m1
87
 
        j = i
88
 
        do 12 j1=1,idim
89
 
          jj = jj+1
90
 
          c(j) = x(jj)
91
 
          j = j+n
92
 
  12    continue
93
 
  15  continue
94
 
      jj = 1
95
 
      j = m
96
 
      do 17 j1=1,idim
97
 
        c(j) = c(jj)
98
 
        j = j+n
99
 
        jj = jj+n
100
 
  17  continue
101
 
      fp = 0.
102
 
      fpint(n) = fp0
103
 
      fpint(n-1) = 0.
104
 
      nrdata(n) = 0
105
 
      go to 630
106
 
  20  do 25 i=2,m1
107
 
         j = i+k
108
 
         t(j) = (u(i)+u(i-1))*half
109
 
  25  continue
110
 
      go to 50
111
 
c  if s > 0 our initial choice depends on the value of iopt.
112
 
c  if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
113
 
c  polynomial curve. (i.e. a constant point).
114
 
c  if iopt=1 and fp0>s we start computing the least-squares closed
115
 
c  curve according the set of knots found at the last call of the
116
 
c  routine.
117
 
  30  if(iopt.eq.0) go to 35
118
 
      if(n.eq.nmin) go to 35
119
 
      fp0 = fpint(n)
120
 
      fpold = fpint(n-1)
121
 
      nplus = nrdata(n)
122
 
      if(fp0.gt.s) go to 50
123
 
c  the case that s(u) is a fixed point is treated separetely.
124
 
c  fp0 denotes the corresponding sum of squared residuals.
125
 
  35  fp0 = 0.
126
 
      d1 = 0.
127
 
      do 37 j=1,idim
128
 
        z(j) = 0.
129
 
  37  continue
130
 
      jj = 0
131
 
      do 45 it=1,m1
132
 
        wi = w(it)
133
 
        call fpgivs(wi,d1,cos,sin)
134
 
        do 40 j=1,idim
135
 
          jj = jj+1
136
 
          fac = wi*x(jj)
137
 
          call fprota(cos,sin,fac,z(j))
138
 
          fp0 = fp0+fac**2
139
 
  40    continue
140
 
  45  continue
141
 
      do 47 j=1,idim
142
 
        z(j) = z(j)/d1
143
 
  47  continue
144
 
c  test whether that fixed point is a solution of our problem.
145
 
      fpms = fp0-s
146
 
      if(fpms.lt.acc .or. nmax.eq.nmin) go to 640
147
 
      fpold = fp0
148
 
c  test whether the required storage space exceeds the available one.
149
 
      if(n.ge.nest) go to 620
150
 
c  start computing the least-squares closed curve with one
151
 
c  interior knot.
152
 
      nplus = 1
153
 
      n = nmin+1
154
 
      mm = (m+1)/2
155
 
      t(k2) = u(mm)
156
 
      nrdata(1) = mm-2
157
 
      nrdata(2) = m1-mm
158
 
c  main loop for the different sets of knots. m is a save upper
159
 
c  bound for the number of trials.
160
 
  50  do 340 iter=1,m
161
 
c  find nrint, the number of knot intervals.
162
 
        nrint = n-nmin+1
163
 
c  find the position of the additional knots which are needed for
164
 
c  the b-spline representation of s(u). if we take
165
 
c      t(k+1) = u(1), t(n-k) = u(m)
166
 
c      t(k+1-j) = t(n-k-j) - per, j=1,2,...k
167
 
c      t(n-k+j) = t(k+1+j) + per, j=1,2,...k
168
 
c  then s(u) will be a smooth closed curve if the b-spline
169
 
c  coefficients satisfy the following conditions
170
 
c      c((i-1)*n+n7+j) = c((i-1)*n+j), j=1,...k,i=1,2,...,idim (**)
171
 
c  with n7=n-2*k-1.
172
 
        t(k1) = u(1)
173
 
        nk1 = n-k1
174
 
        nk2 = nk1+1
175
 
        t(nk2) = u(m)
176
 
        do 60 j=1,k
177
 
          i1 = nk2+j
178
 
          i2 = nk2-j
179
 
          j1 = k1+j
180
 
          j2 = k1-j
181
 
          t(i1) = t(j1)+per
182
 
          t(j2) = t(i2)-per
183
 
  60    continue
184
 
c  compute the b-spline coefficients of the least-squares closed curve
185
 
c  sinf(u). the observation matrix a is built up row by row while
186
 
c  taking into account condition (**) and is reduced to triangular
187
 
c  form by givens transformations .
188
 
c  at the same time fp=f(p=inf) is computed.
189
 
c  the n7 x n7 triangularised upper matrix a has the form
190
 
c            ! a1 '    !
191
 
c        a = !    ' a2 !
192
 
c            ! 0  '    !
193
 
c  with a2 a n7 x k matrix and a1 a n10 x n10 upper triangular
194
 
c  matrix of bandwith k+1 ( n10 = n7-k).
195
 
c  initialization.
196
 
        do 65 i=1,nc
197
 
          z(i) = 0.
198
 
  65    continue
199
 
        do 70 i=1,nk1
200
 
          do 70 j=1,kk1
201
 
            a1(i,j) = 0.
202
 
  70    continue
203
 
        n7 = nk1-k
204
 
        n10 = n7-kk
205
 
        jper = 0
206
 
        fp = 0.
207
 
        l = k1
208
 
        jj = 0
209
 
        do 290 it=1,m1
210
 
c  fetch the current data point u(it),x(it)
211
 
          ui = u(it)
212
 
          wi = w(it)
213
 
          do 75 j=1,idim
214
 
            jj = jj+1
215
 
            xi(j) = x(jj)*wi
216
 
  75      continue
217
 
c  search for knot interval t(l) <= ui < t(l+1).
218
 
  80      if(ui.lt.t(l+1)) go to 85
219
 
          l = l+1
220
 
          go to 80
221
 
c  evaluate the (k+1) non-zero b-splines at ui and store them in q.
222
 
  85      call fpbspl(t,n,k,ui,l,h)
223
 
          do 90 i=1,k1
224
 
            q(it,i) = h(i)
225
 
            h(i) = h(i)*wi
226
 
  90      continue
227
 
          l5 = l-k1
228
 
c  test whether the b-splines nj,k+1(u),j=1+n7,...nk1 are all zero at ui
229
 
          if(l5.lt.n10) go to 285
230
 
          if(jper.ne.0) go to 160
231
 
c  initialize the matrix a2.
232
 
          do 95 i=1,n7
233
 
          do 95 j=1,kk
234
 
              a2(i,j) = 0.
235
 
  95      continue
236
 
          jk = n10+1
237
 
          do 110 i=1,kk
238
 
            ik = jk
239
 
            do 100 j=1,kk1
240
 
              if(ik.le.0) go to 105
241
 
              a2(ik,i) = a1(ik,j)
242
 
              ik = ik-1
243
 
 100        continue
244
 
 105        jk = jk+1
245
 
 110      continue
246
 
          jper = 1
247
 
c  if one of the b-splines nj,k+1(u),j=n7+1,...nk1 is not zero at ui
248
 
c  we take account of condition (**) for setting up the new row
249
 
c  of the observation matrix a. this row is stored in the arrays h1
250
 
c  (the part with respect to a1) and h2 (the part with
251
 
c  respect to a2).
252
 
 160      do 170 i=1,kk
253
 
            h1(i) = 0.
254
 
            h2(i) = 0.
255
 
 170      continue
256
 
          h1(kk1) = 0.
257
 
          j = l5-n10
258
 
          do 210 i=1,kk1
259
 
            j = j+1
260
 
            l0 = j
261
 
 180        l1 = l0-kk
262
 
            if(l1.le.0) go to 200
263
 
            if(l1.le.n10) go to 190
264
 
            l0 = l1-n10
265
 
            go to 180
266
 
 190        h1(l1) = h(i)
267
 
            go to 210
268
 
 200        h2(l0) = h2(l0)+h(i)
269
 
 210      continue
270
 
c  rotate the new row of the observation matrix into triangle
271
 
c  by givens transformations.
272
 
          if(n10.le.0) go to 250
273
 
c  rotation with the rows 1,2,...n10 of matrix a.
274
 
          do 240 j=1,n10
275
 
            piv = h1(1)
276
 
            if(piv.ne.0.) go to 214
277
 
            do 212 i=1,kk
278
 
              h1(i) = h1(i+1)
279
 
 212        continue
280
 
            h1(kk1) = 0.
281
 
            go to 240
282
 
c  calculate the parameters of the givens transformation.
283
 
 214        call fpgivs(piv,a1(j,1),cos,sin)
284
 
c  transformation to the right hand side.
285
 
            j1 = j
286
 
            do 217 j2=1,idim
287
 
              call fprota(cos,sin,xi(j2),z(j1))
288
 
              j1 = j1+n
289
 
 217        continue
290
 
c  transformations to the left hand side with respect to a2.
291
 
            do 220 i=1,kk
292
 
              call fprota(cos,sin,h2(i),a2(j,i))
293
 
 220        continue
294
 
            if(j.eq.n10) go to 250
295
 
            i2 = min0(n10-j,kk)
296
 
c  transformations to the left hand side with respect to a1.
297
 
            do 230 i=1,i2
298
 
              i1 = i+1
299
 
              call fprota(cos,sin,h1(i1),a1(j,i1))
300
 
              h1(i) = h1(i1)
301
 
 230        continue
302
 
            h1(i1) = 0.
303
 
 240      continue
304
 
c  rotation with the rows n10+1,...n7 of matrix a.
305
 
 250      do 270 j=1,kk
306
 
            ij = n10+j
307
 
            if(ij.le.0) go to 270
308
 
            piv = h2(j)
309
 
            if(piv.eq.0.) go to 270
310
 
c  calculate the parameters of the givens transformation.
311
 
            call fpgivs(piv,a2(ij,j),cos,sin)
312
 
c  transformations to right hand side.
313
 
            j1 = ij
314
 
            do 255 j2=1,idim
315
 
              call fprota(cos,sin,xi(j2),z(j1))
316
 
              j1 = j1+n
317
 
 255        continue
318
 
            if(j.eq.kk) go to 280
319
 
            j1 = j+1
320
 
c  transformations to left hand side.
321
 
            do 260 i=j1,kk
322
 
              call fprota(cos,sin,h2(i),a2(ij,i))
323
 
 260        continue
324
 
 270      continue
325
 
c  add contribution of this row to the sum of squares of residual
326
 
c  right hand sides.
327
 
 280      do 282 j2=1,idim
328
 
            fp = fp+xi(j2)**2
329
 
 282      continue
330
 
          go to 290
331
 
c  rotation of the new row of the observation matrix into
332
 
c  triangle in case the b-splines nj,k+1(u),j=n7+1,...n-k-1 are all zero
333
 
c  at ui.
334
 
 285      j = l5
335
 
          do 140 i=1,kk1
336
 
            j = j+1
337
 
            piv = h(i)
338
 
            if(piv.eq.0.) go to 140
339
 
c  calculate the parameters of the givens transformation.
340
 
            call fpgivs(piv,a1(j,1),cos,sin)
341
 
c  transformations to right hand side.
342
 
            j1 = j
343
 
            do 125 j2=1,idim
344
 
              call fprota(cos,sin,xi(j2),z(j1))
345
 
              j1 = j1+n
346
 
 125        continue
347
 
            if(i.eq.kk1) go to 150
348
 
            i2 = 1
349
 
            i3 = i+1
350
 
c  transformations to left hand side.
351
 
            do 130 i1=i3,kk1
352
 
              i2 = i2+1
353
 
              call fprota(cos,sin,h(i1),a1(j,i2))
354
 
 130        continue
355
 
 140      continue
356
 
c  add contribution of this row to the sum of squares of residual
357
 
c  right hand sides.
358
 
 150      do 155 j2=1,idim
359
 
            fp = fp+xi(j2)**2
360
 
 155      continue
361
 
 290    continue
362
 
        fpint(n) = fp0
363
 
        fpint(n-1) = fpold
364
 
        nrdata(n) = nplus
365
 
c  backward substitution to obtain the b-spline coefficients .
366
 
        j1 = 1
367
 
        do 292 j2=1,idim
368
 
           call fpbacp(a1,a2,z(j1),n7,kk,c(j1),kk1,nest)
369
 
           j1 = j1+n
370
 
 292    continue
371
 
c  calculate from condition (**) the remaining coefficients.
372
 
        do 297 i=1,k
373
 
          j1 = i
374
 
          do 295 j=1,idim
375
 
            j2 = j1+n7
376
 
            c(j2) = c(j1)
377
 
            j1 = j1+n
378
 
 295      continue
379
 
 297    continue
380
 
        if(iopt.lt.0) go to 660
381
 
c  test whether the approximation sinf(u) is an acceptable solution.
382
 
        fpms = fp-s
383
 
        if(abs(fpms).lt.acc) go to 660
384
 
c  if f(p=inf) < s accept the choice of knots.
385
 
        if(fpms.lt.0.) go to 350
386
 
c  if n=nmax, sinf(u) is an interpolating curve.
387
 
        if(n.eq.nmax) go to 630
388
 
c  increase the number of knots.
389
 
c  if n=nest we cannot increase the number of knots because of the
390
 
c  storage capacity limitation.
391
 
        if(n.eq.nest) go to 620
392
 
c  determine the number of knots nplus we are going to add.
393
 
        npl1 = nplus*2
394
 
        rn = nplus
395
 
        if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
396
 
        nplus = min0(nplus*2,max0(npl1,nplus/2,1))
397
 
        fpold = fp
398
 
c  compute the sum of squared residuals for each knot interval
399
 
c  t(j+k) <= ui <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
400
 
        fpart = 0.
401
 
        i = 1
402
 
        l = k1
403
 
        jj = 0
404
 
        do 320 it=1,m1
405
 
          if(u(it).lt.t(l)) go to 300
406
 
          new = 1
407
 
          l = l+1
408
 
 300      term = 0.
409
 
          l0 = l-k2
410
 
          do 310 j2=1,idim
411
 
            fac = 0.
412
 
            j1 = l0
413
 
            do 305 j=1,k1
414
 
              j1 = j1+1
415
 
              fac = fac+c(j1)*q(it,j)
416
 
 305        continue
417
 
            jj = jj+1
418
 
            term = term+(w(it)*(fac-x(jj)))**2
419
 
            l0 = l0+n
420
 
 310      continue
421
 
          fpart = fpart+term
422
 
          if(new.eq.0) go to 320
423
 
          if(l.gt.k2) go to 315
424
 
          fpint(nrint) = term
425
 
          new = 0
426
 
          go to 320
427
 
 315      store = term*half
428
 
          fpint(i) = fpart-store
429
 
          i = i+1
430
 
          fpart = store
431
 
          new = 0
432
 
 320    continue
433
 
        fpint(nrint) = fpint(nrint)+fpart
434
 
        do 330 l=1,nplus
435
 
c  add a new knot
436
 
          call fpknot(u,m,t,n,fpint,nrdata,nrint,nest,1)
437
 
c  if n=nmax we locate the knots as for interpolation
438
 
          if(n.eq.nmax) go to 5
439
 
c  test whether we cannot further increase the number of knots.
440
 
          if(n.eq.nest) go to 340
441
 
 330    continue
442
 
c  restart the computations with the new set of knots.
443
 
 340  continue
444
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
445
 
c  part 2: determination of the smoothing closed curve sp(u).          c
446
 
c  **********************************************************          c
447
 
c  we have determined the number of knots and their position.          c
448
 
c  we now compute the b-spline coefficients of the smoothing curve     c
449
 
c  sp(u). the observation matrix a is extended by the rows of matrix   c
450
 
c  b expressing that the kth derivative discontinuities of sp(u) at    c
451
 
c  the interior knots t(k+2),...t(n-k-1) must be zero. the corres-     c
452
 
c  ponding weights of these additional rows are set to 1/p.            c
453
 
c  iteratively we then have to determine the value of p such that f(p),c
454
 
c  the sum of squared residuals be = s. we already know that the least-c
455
 
c  squares polynomial curve corresponds to p=0, and that the least-    c
456
 
c  squares periodic spline curve corresponds to p=infinity. the        c
457
 
c  iteration process which is proposed here, makes use of rational     c
458
 
c  interpolation. since f(p) is a convex and strictly decreasing       c
459
 
c  function of p, it can be approximated by a rational function        c
460
 
c  r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond-  c
461
 
c  ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used      c
462
 
c  to calculate the new value of p such that r(p)=s. convergence is    c
463
 
c  guaranteed by taking f1>0 and f3<0.                                 c
464
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
465
 
c  evaluate the discontinuity jump of the kth derivative of the
466
 
c  b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
467
 
 350  call fpdisc(t,n,k2,b,nest)
468
 
c  initial value for p.
469
 
      p1 = 0.
470
 
      f1 = fp0-s
471
 
      p3 = -one
472
 
      f3 = fpms
473
 
      n11 = n10-1
474
 
      n8 = n7-1
475
 
      p = 0.
476
 
      l = n7
477
 
      do 352 i=1,k
478
 
         j = k+1-i
479
 
         p = p+a2(l,j)
480
 
         l = l-1
481
 
         if(l.eq.0) go to 356
482
 
 352  continue
483
 
      do 354 i=1,n10
484
 
         p = p+a1(i,1)
485
 
 354  continue
486
 
 356  rn = n7
487
 
      p = rn/p
488
 
      ich1 = 0
489
 
      ich3 = 0
490
 
c  iteration process to find the root of f(p) = s.
491
 
      do 595 iter=1,maxit
492
 
c  form the matrix g  as the matrix a extended by the rows of matrix b.
493
 
c  the rows of matrix b with weight 1/p are rotated into
494
 
c  the triangularised observation matrix a.
495
 
c  after triangularisation our n7 x n7 matrix g takes the form
496
 
c            ! g1 '    !
497
 
c        g = !    ' g2 !
498
 
c            ! 0  '    !
499
 
c  with g2 a n7 x (k+1) matrix and g1 a n11 x n11 upper triangular
500
 
c  matrix of bandwidth k+2. ( n11 = n7-k-1)
501
 
        pinv = one/p
502
 
c  store matrix a into g
503
 
        do 358 i=1,nc
504
 
          c(i) = z(i)
505
 
 358    continue
506
 
        do 360 i=1,n7
507
 
          g1(i,k1) = a1(i,k1)
508
 
          g1(i,k2) = 0.
509
 
          g2(i,1) = 0.
510
 
          do 360 j=1,k
511
 
            g1(i,j) = a1(i,j)
512
 
            g2(i,j+1) = a2(i,j)
513
 
 360    continue
514
 
        l = n10
515
 
        do 370 j=1,k1
516
 
          if(l.le.0) go to 375
517
 
          g2(l,1) = a1(l,j)
518
 
          l = l-1
519
 
 370    continue
520
 
 375    do 540 it=1,n8
521
 
c  fetch a new row of matrix b and store it in the arrays h1 (the part
522
 
c  with respect to g1) and h2 (the part with respect to g2).
523
 
          do 380 j=1,idim
524
 
            xi(j) = 0.
525
 
 380      continue
526
 
          do 385 i=1,k1
527
 
            h1(i) = 0.
528
 
            h2(i) = 0.
529
 
 385      continue
530
 
          h1(k2) = 0.
531
 
          if(it.gt.n11) go to 420
532
 
          l = it
533
 
          l0 = it
534
 
          do 390 j=1,k2
535
 
            if(l0.eq.n10) go to 400
536
 
            h1(j) = b(it,j)*pinv
537
 
            l0 = l0+1
538
 
 390      continue
539
 
          go to 470
540
 
 400      l0 = 1
541
 
          do 410 l1=j,k2
542
 
            h2(l0) = b(it,l1)*pinv
543
 
            l0 = l0+1
544
 
 410      continue
545
 
          go to 470
546
 
 420      l = 1
547
 
          i = it-n10
548
 
          do 460 j=1,k2
549
 
            i = i+1
550
 
            l0 = i
551
 
 430        l1 = l0-k1
552
 
            if(l1.le.0) go to 450
553
 
            if(l1.le.n11) go to 440
554
 
            l0 = l1-n11
555
 
            go to 430
556
 
 440        h1(l1) = b(it,j)*pinv
557
 
            go to 460
558
 
 450        h2(l0) = h2(l0)+b(it,j)*pinv
559
 
 460      continue
560
 
          if(n11.le.0) go to 510
561
 
c  rotate this row into triangle by givens transformations
562
 
c  rotation with the rows l,l+1,...n11.
563
 
 470      do 500 j=l,n11
564
 
            piv = h1(1)
565
 
c  calculate the parameters of the givens transformation.
566
 
            call fpgivs(piv,g1(j,1),cos,sin)
567
 
c  transformation to right hand side.
568
 
            j1 = j
569
 
            do 475 j2=1,idim
570
 
              call fprota(cos,sin,xi(j2),c(j1))
571
 
              j1 = j1+n
572
 
 475        continue
573
 
c  transformation to the left hand side with respect to g2.
574
 
            do 480 i=1,k1
575
 
              call fprota(cos,sin,h2(i),g2(j,i))
576
 
 480        continue
577
 
            if(j.eq.n11) go to 510
578
 
            i2 = min0(n11-j,k1)
579
 
c  transformation to the left hand side with respect to g1.
580
 
            do 490 i=1,i2
581
 
              i1 = i+1
582
 
              call fprota(cos,sin,h1(i1),g1(j,i1))
583
 
              h1(i) = h1(i1)
584
 
 490        continue
585
 
            h1(i1) = 0.
586
 
 500      continue
587
 
c  rotation with the rows n11+1,...n7
588
 
 510      do 530 j=1,k1
589
 
            ij = n11+j
590
 
            if(ij.le.0) go to 530
591
 
            piv = h2(j)
592
 
c  calculate the parameters of the givens transformation
593
 
            call fpgivs(piv,g2(ij,j),cos,sin)
594
 
c  transformation to the right hand side.
595
 
            j1 = ij
596
 
            do 515 j2=1,idim
597
 
              call fprota(cos,sin,xi(j2),c(j1))
598
 
              j1 = j1+n
599
 
 515        continue
600
 
            if(j.eq.k1) go to 540
601
 
            j1 = j+1
602
 
c  transformation to the left hand side.
603
 
            do 520 i=j1,k1
604
 
              call fprota(cos,sin,h2(i),g2(ij,i))
605
 
 520        continue
606
 
 530      continue
607
 
 540    continue
608
 
c  backward substitution to obtain the b-spline coefficients
609
 
        j1 = 1
610
 
        do 542 j2=1,idim
611
 
          call fpbacp(g1,g2,c(j1),n7,k1,c(j1),k2,nest)
612
 
          j1 = j1+n
613
 
 542    continue
614
 
c  calculate from condition (**) the remaining b-spline coefficients.
615
 
        do 547 i=1,k
616
 
          j1 = i
617
 
          do 545 j=1,idim
618
 
            j2 = j1+n7
619
 
            c(j2) = c(j1)
620
 
            j1 = j1+n
621
 
 545      continue
622
 
 547    continue
623
 
c  computation of f(p).
624
 
        fp = 0.
625
 
        l = k1
626
 
        jj = 0
627
 
        do 570 it=1,m1
628
 
          if(u(it).lt.t(l)) go to 550
629
 
          l = l+1
630
 
 550      l0 = l-k2
631
 
          term = 0.
632
 
          do 565 j2=1,idim
633
 
            fac = 0.
634
 
            j1 = l0
635
 
            do 560 j=1,k1
636
 
              j1 = j1+1
637
 
              fac = fac+c(j1)*q(it,j)
638
 
 560        continue
639
 
            jj = jj+1
640
 
            term = term+(fac-x(jj))**2
641
 
            l0 = l0+n
642
 
 565      continue
643
 
          fp = fp+term*w(it)**2
644
 
 570    continue
645
 
c  test whether the approximation sp(u) is an acceptable solution.
646
 
        fpms = fp-s
647
 
        if(abs(fpms).lt.acc) go to 660
648
 
c  test whether the maximal number of iterations is reached.
649
 
        if(iter.eq.maxit) go to 600
650
 
c  carry out one more step of the iteration process.
651
 
        p2 = p
652
 
        f2 = fpms
653
 
        if(ich3.ne.0) go to 580
654
 
        if((f2-f3) .gt. acc) go to 575
655
 
c  our initial choice of p is too large.
656
 
        p3 = p2
657
 
        f3 = f2
658
 
        p = p*con4
659
 
        if(p.le.p1) p = p1*con9 +p2*con1
660
 
        go to 595
661
 
 575    if(f2.lt.0.) ich3 = 1
662
 
 580    if(ich1.ne.0) go to 590
663
 
        if((f1-f2) .gt. acc) go to 585
664
 
c  our initial choice of p is too small
665
 
        p1 = p2
666
 
        f1 = f2
667
 
        p = p/con4
668
 
        if(p3.lt.0.) go to 595
669
 
        if(p.ge.p3) p = p2*con1 +p3*con9
670
 
        go to 595
671
 
 585    if(f2.gt.0.) ich1 = 1
672
 
c  test whether the iteration process proceeds as theoretically
673
 
c  expected.
674
 
 590    if(f2.ge.f1 .or. f2.le.f3) go to 610
675
 
c  find the new value for p.
676
 
        p = fprati(p1,f1,p2,f2,p3,f3)
677
 
 595  continue
678
 
c  error codes and messages.
679
 
 600  ier = 3
680
 
      go to 660
681
 
 610  ier = 2
682
 
      go to 660
683
 
 620  ier = 1
684
 
      go to 660
685
 
 630  ier = -1
686
 
      go to 660
687
 
 640  ier = -2
688
 
c  the point (z(1),z(2),...,z(idim)) is a solution of our problem.
689
 
c  a constant function is a spline of degree k with all b-spline
690
 
c  coefficients equal to that constant.
691
 
      do 650 i=1,k1
692
 
        rn = k1-i
693
 
        t(i) = u(1)-rn*per
694
 
        j = i+k1
695
 
        rn = i-1
696
 
        t(j) = u(m)+rn*per
697
 
 650  continue
698
 
      n = nmin
699
 
      j1 = 0
700
 
      do 658 j=1,idim
701
 
        fac = z(j)
702
 
        j2 = j1
703
 
        do 654 i=1,k1
704
 
          j2 = j2+1
705
 
          c(j2) = fac
706
 
 654    continue
707
 
        j1 = j1+n
708
 
 658  continue
709
 
      fp = fp0
710
 
      fpint(n) = fp0
711
 
      fpint(n-1) = 0.
712
 
      nrdata(n) = 0
713
 
 660  return
714
 
      end