← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/spline/fitpack/fpgrsp.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/sandbox/spline/fitpack/fpgrsp.f
1
 
      subroutine fpgrsp(ifsu,ifsv,ifbu,ifbv,iback,u,mu,v,mv,r,mr,dr,
2
 
     * iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,spu,spv,
3
 
     * right,q,au,av1,av2,bu,bv,a0,a1,b0,b1,c0,c1,cosi,nru,nrv)
4
 
c  ..
5
 
c  ..scalar arguments..
6
 
      real*8 p,sq,fp
7
 
      integer ifsu,ifsv,ifbu,ifbv,iback,mu,mv,mr,iop0,iop1,nu,nv,nc,
8
 
     * mm,mvnu
9
 
c  ..array arguments..
10
 
      real*8 u(mu),v(mv),r(mr),dr(6),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv)
11
 
     *,
12
 
     * spu(mu,4),spv(mv,4),right(mm),q(mvnu),au(nu,5),av1(nv,6),c0(nv),
13
 
     * av2(nv,4),a0(2,mv),b0(2,nv),cosi(2,nv),bu(nu,5),bv(nv,5),c1(nv),
14
 
     * a1(2,mv),b1(2,nv)
15
 
      integer nru(mu),nrv(mv)
16
 
c  ..local scalars..
17
 
      real*8 arg,co,dr01,dr02,dr03,dr11,dr12,dr13,fac,fac0,fac1,pinv,piv
18
 
     *,
19
 
     * si,term,one,three,half
20
 
      integer i,ic,ii,ij,ik,iq,irot,it,ir,i0,i1,i2,i3,j,jj,jk,jper,
21
 
     * j0,j1,k,k1,k2,l,l0,l1,l2,mvv,ncof,nrold,nroldu,nroldv,number,
22
 
     * numu,numu1,numv,numv1,nuu,nu4,nu7,nu8,nu9,nv11,nv4,nv7,nv8,n1
23
 
c  ..local arrays..
24
 
      real*8 h(5),h1(5),h2(4)
25
 
c  ..function references..
26
 
      integer min0
27
 
      real*8 cos,sin
28
 
c  ..subroutine references..
29
 
c    fpback,fpbspl,fpgivs,fpcyt1,fpcyt2,fpdisc,fpbacp,fprota
30
 
c  ..
31
 
c  let
32
 
c               |     (spu)      |            |     (spv)      |
33
 
c        (au) = | -------------- |     (av) = | -------------- |
34
 
c               | sqrt(1/p) (bu) |            | sqrt(1/p) (bv) |
35
 
c
36
 
c                                | r  ' 0 |
37
 
c                            q = | ------ |
38
 
c                                | 0  ' 0 |
39
 
c
40
 
c  with c      : the (nu-4) x (nv-4) matrix which contains the b-spline
41
 
c                coefficients.
42
 
c       r      : the mu x mv matrix which contains the function values.
43
 
c       spu,spv: the mu x (nu-4), resp. mv x (nv-4) observation matrices
44
 
c                according to the least-squares problems in the u-,resp.
45
 
c                v-direction.
46
 
c       bu,bv  : the (nu-7) x (nu-4),resp. (nv-7) x (nv-4) matrices
47
 
c                containing the discontinuity jumps of the derivatives
48
 
c                of the b-splines in the u-,resp.v-variable at the knots
49
 
c  the b-spline coefficients of the smoothing spline are then calculated
50
 
c  as the least-squares solution of the following over-determined linear
51
 
c  system of equations
52
 
c
53
 
c  (1)  (av) c (au)' = q
54
 
c
55
 
c  subject to the constraints
56
 
c
57
 
c  (2)  c(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4
58
 
c
59
 
c  (3)  if iop0 = 0  c(1,j) = dr(1)
60
 
c          iop0 = 1  c(1,j) = dr(1)
61
 
c                    c(2,j) = dr(1)+(dr(2)*cosi(1,j)+dr(3)*cosi(2,j))*
62
 
c                            tu(5)/3. = c0(j) , j=1,2,...nv-4
63
 
c
64
 
c  (4)  if iop1 = 0  c(nu-4,j) = dr(4)
65
 
c          iop1 = 1  c(nu-4,j) = dr(4)
66
 
c                    c(nu-5,j) = dr(4)+(dr(5)*cosi(1,j)+dr(6)*cosi(2,j))
67
 
c                                *(tu(nu-4)-tu(nu-3))/3. = c1(j)
68
 
c
69
 
c  set constants
70
 
      one = 1
71
 
      three = 3
72
 
      half = 0.5
73
 
c  initialization
74
 
      nu4 = nu-4
75
 
      nu7 = nu-7
76
 
      nu8 = nu-8
77
 
      nu9 = nu-9
78
 
      nv4 = nv-4
79
 
      nv7 = nv-7
80
 
      nv8 = nv-8
81
 
      nv11 = nv-11
82
 
      nuu = nu4-iop0-iop1-2
83
 
      if(p.gt.0.) pinv = one/p
84
 
c  it depends on the value of the flags ifsu,ifsv,ifbu,ifbv,iop0,iop1
85
 
c  and on the value of p whether the matrices (spu), (spv), (bu), (bv),
86
 
c  (cosi) still must be determined.
87
 
      if(ifsu.ne.0) go to 30
88
 
c  calculate the non-zero elements of the matrix (spu) which is the ob-
89
 
c  servation matrix according to the least-squares spline approximation
90
 
c  problem in the u-direction.
91
 
      l = 4
92
 
      l1 = 5
93
 
      number = 0
94
 
      do 25 it=1,mu
95
 
        arg = u(it)
96
 
  10    if(arg.lt.tu(l1) .or. l.eq.nu4) go to 15
97
 
        l = l1
98
 
        l1 = l+1
99
 
        number = number+1
100
 
        go to 10
101
 
  15    call fpbspl(tu,nu,3,arg,l,h)
102
 
        do 20 i=1,4
103
 
          spu(it,i) = h(i)
104
 
  20    continue
105
 
        nru(it) = number
106
 
  25  continue
107
 
      ifsu = 1
108
 
c  calculate the non-zero elements of the matrix (spv) which is the ob-
109
 
c  servation matrix according to the least-squares spline approximation
110
 
c  problem in the v-direction.
111
 
  30  if(ifsv.ne.0) go to 85
112
 
      l = 4
113
 
      l1 = 5
114
 
      number = 0
115
 
      do 50 it=1,mv
116
 
        arg = v(it)
117
 
  35    if(arg.lt.tv(l1) .or. l.eq.nv4) go to 40
118
 
        l = l1
119
 
        l1 = l+1
120
 
        number = number+1
121
 
        go to 35
122
 
  40    call fpbspl(tv,nv,3,arg,l,h)
123
 
        do 45 i=1,4
124
 
          spv(it,i) = h(i)
125
 
  45    continue
126
 
        nrv(it) = number
127
 
  50  continue
128
 
      ifsv = 1
129
 
      if(iop0.eq.0 .and. iop1.eq.0) go to 85
130
 
c  calculate the coefficients of the interpolating splines for cos(v)
131
 
c  and sin(v).
132
 
      do 55 i=1,nv4
133
 
         cosi(1,i) = 0.
134
 
         cosi(2,i) = 0.
135
 
  55  continue
136
 
      if(nv7.lt.4) go to 85
137
 
      do 65 i=1,nv7
138
 
         l = i+3
139
 
         arg = tv(l)
140
 
         call fpbspl(tv,nv,3,arg,l,h)
141
 
         do 60 j=1,3
142
 
            av1(i,j) = h(j)
143
 
  60     continue
144
 
         cosi(1,i) = cos(arg)
145
 
         cosi(2,i) = sin(arg)
146
 
  65  continue
147
 
      call fpcyt1(av1,nv7,nv)
148
 
      do 80 j=1,2
149
 
         do 70 i=1,nv7
150
 
            right(i) = cosi(j,i)
151
 
  70     continue
152
 
         call fpcyt2(av1,nv7,right,right,nv)
153
 
         do 75 i=1,nv7
154
 
            cosi(j,i+1) = right(i)
155
 
  75     continue
156
 
         cosi(j,1) = cosi(j,nv7+1)
157
 
         cosi(j,nv7+2) = cosi(j,2)
158
 
         cosi(j,nv4) = cosi(j,3)
159
 
  80  continue
160
 
  85  if(p.le.0.) go to  150
161
 
c  calculate the non-zero elements of the matrix (bu).
162
 
      if(ifbu.ne.0 .or. nu8.eq.0) go to 90
163
 
      call fpdisc(tu,nu,5,bu,nu)
164
 
      ifbu = 1
165
 
c  calculate the non-zero elements of the matrix (bv).
166
 
  90  if(ifbv.ne.0 .or. nv8.eq.0) go to 150
167
 
      call fpdisc(tv,nv,5,bv,nv)
168
 
      ifbv = 1
169
 
c  substituting (2),(3) and (4) into (1), we obtain the overdetermined
170
 
c  system
171
 
c         (5)  (avv) (cc) (auu)' = (qq)
172
 
c  from which the nuu*nv7 remaining coefficients
173
 
c         c(i,j) , i=2+iop0,3+iop0,...,nu-5-iop1,j=1,2,...,nv-7.
174
 
c  the elements of (cc), are then determined in the least-squares sense.
175
 
c  simultaneously, we compute the resulting sum of squared residuals sq.
176
 
 150  dr01 = dr(1)
177
 
      dr11 = dr(4)
178
 
      do 155 i=1,mv
179
 
         a0(1,i) = dr01
180
 
         a1(1,i) = dr11
181
 
 155  continue
182
 
      if(nv8.eq.0 .or. p.le.0.) go to 165
183
 
      do 160 i=1,nv8
184
 
         b0(1,i) = 0.
185
 
         b1(1,i) = 0.
186
 
 160  continue
187
 
 165  mvv = mv
188
 
      if(iop0.eq.0) go to 195
189
 
      fac = (tu(5)-tu(4))/three
190
 
      dr02 = dr(2)*fac
191
 
      dr03 = dr(3)*fac
192
 
      do 170 i=1,nv4
193
 
         c0(i) = dr01+dr02*cosi(1,i)+dr03*cosi(2,i)
194
 
 170  continue
195
 
      do 180 i=1,mv
196
 
         number = nrv(i)
197
 
         fac = 0.
198
 
         do 175 j=1,4
199
 
            number = number+1
200
 
            fac = fac+c0(number)*spv(i,j)
201
 
 175     continue
202
 
         a0(2,i) = fac
203
 
 180  continue
204
 
      if(nv8.eq.0 .or. p.le.0.) go to 195
205
 
      do 190 i=1,nv8
206
 
         number = i
207
 
         fac = 0.
208
 
         do 185 j=1,5
209
 
            fac = fac+c0(number)*bv(i,j)
210
 
            number = number+1
211
 
 185     continue
212
 
         b0(2,i) = fac*pinv
213
 
 190  continue
214
 
      mvv = mv+nv8
215
 
 195  if(iop1.eq.0) go to 225
216
 
      fac = (tu(nu4)-tu(nu4+1))/three
217
 
      dr12 = dr(5)*fac
218
 
      dr13 = dr(6)*fac
219
 
      do 200 i=1,nv4
220
 
         c1(i) = dr11+dr12*cosi(1,i)+dr13*cosi(2,i)
221
 
 200  continue
222
 
      do 210 i=1,mv
223
 
         number = nrv(i)
224
 
         fac = 0.
225
 
         do 205 j=1,4
226
 
            number = number+1
227
 
            fac = fac+c1(number)*spv(i,j)
228
 
 205     continue
229
 
         a1(2,i) = fac
230
 
 210  continue
231
 
      if(nv8.eq.0 .or. p.le.0.) go to 225
232
 
      do 220 i=1,nv8
233
 
         number = i
234
 
         fac = 0.
235
 
         do 215 j=1,5
236
 
            fac = fac+c1(number)*bv(i,j)
237
 
            number = number+1
238
 
 215     continue
239
 
         b1(2,i) = fac*pinv
240
 
 220  continue
241
 
      mvv = mv+nv8
242
 
c  we first determine the matrices (auu) and (qq). then we reduce the
243
 
c  matrix (auu) to an unit upper triangular form (ru) using givens
244
 
c  rotations without square roots. we apply the same transformations to
245
 
c  the rows of matrix qq to obtain the mv x nuu matrix g.
246
 
c  we store matrix (ru) into au and g into q.
247
 
 225  l = mvv*nuu
248
 
c  initialization.
249
 
      sq = 0.
250
 
      if(l.eq.0) go to 245
251
 
      do 230 i=1,l
252
 
        q(i) = 0.
253
 
 230  continue
254
 
      do 240 i=1,nuu
255
 
        do 240 j=1,5
256
 
          au(i,j) = 0.
257
 
 240  continue
258
 
      l = 0
259
 
 245  nrold = 0
260
 
      n1 = nrold+1
261
 
      do 420 it=1,mu
262
 
        number = nru(it)
263
 
c  find the appropriate column of q.
264
 
 250    do 260 j=1,mvv
265
 
           right(j) = 0.
266
 
 260    continue
267
 
        if(nrold.eq.number) go to 280
268
 
        if(p.le.0.) go to 410
269
 
c  fetch a new row of matrix (bu).
270
 
        do 270 j=1,5
271
 
          h(j) = bu(n1,j)*pinv
272
 
 270    continue
273
 
        i0 = 1
274
 
        i1 = 5
275
 
        go to 310
276
 
c  fetch a new row of matrix (spu).
277
 
 280    do 290 j=1,4
278
 
          h(j) = spu(it,j)
279
 
 290    continue
280
 
c  find the appropriate column of q.
281
 
        do 300 j=1,mv
282
 
          l = l+1
283
 
          right(j) = r(l)
284
 
 300    continue
285
 
        i0 = 1
286
 
        i1 = 4
287
 
 310    j0 = n1
288
 
        j1 = nu7-number
289
 
c  take into account that we eliminate the constraints (3)
290
 
 315     if(j0-1.gt.iop0) go to 335
291
 
         fac0 = h(i0)
292
 
         do 320 j=1,mv
293
 
            right(j) = right(j)-fac0*a0(j0,j)
294
 
 320     continue
295
 
         if(mv.eq.mvv) go to 330
296
 
         j = mv
297
 
         do 325 jj=1,nv8
298
 
            j = j+1
299
 
            right(j) = right(j)-fac0*b0(j0,jj)
300
 
 325     continue
301
 
 330     j0 = j0+1
302
 
         i0 = i0+1
303
 
         go to 315
304
 
c  take into account that we eliminate the constraints (4)
305
 
 335     if(j1-1.gt.iop1) go to 360
306
 
         fac1 = h(i1)
307
 
         do 340 j=1,mv
308
 
            right(j) = right(j)-fac1*a1(j1,j)
309
 
 340     continue
310
 
         if(mv.eq.mvv) go to 350
311
 
         j = mv
312
 
         do 345 jj=1,nv8
313
 
            j = j+1
314
 
            right(j) = right(j)-fac1*b1(j1,jj)
315
 
 345     continue
316
 
 350     j1 = j1+1
317
 
         i1 = i1-1
318
 
         go to 335
319
 
 360     irot = nrold-iop0-1
320
 
         if(irot.lt.0) irot = 0
321
 
c  rotate the new row of matrix (auu) into triangle.
322
 
        if(i0.gt.i1) go to 390
323
 
        do 385 i=i0,i1
324
 
          irot = irot+1
325
 
          piv = h(i)
326
 
          if(piv.eq.0.) go to 385
327
 
c  calculate the parameters of the givens transformation.
328
 
          call fpgivs(piv,au(irot,1),co,si)
329
 
c  apply that transformation to the rows of matrix (qq).
330
 
          iq = (irot-1)*mvv
331
 
          do 370 j=1,mvv
332
 
            iq = iq+1
333
 
            call fprota(co,si,right(j),q(iq))
334
 
 370      continue
335
 
c  apply that transformation to the columns of (auu).
336
 
          if(i.eq.i1) go to 385
337
 
          i2 = 1
338
 
          i3 = i+1
339
 
          do 380 j=i3,i1
340
 
            i2 = i2+1
341
 
            call fprota(co,si,h(j),au(irot,i2))
342
 
 380      continue
343
 
 385    continue
344
 
c  we update the sum of squared residuals.
345
 
 390    do 395 j=1,mvv
346
 
          sq = sq+right(j)**2
347
 
 395    continue
348
 
 400    if(nrold.eq.number) go to 420
349
 
 410    nrold = n1
350
 
        n1 = n1+1
351
 
        go to 250
352
 
 420  continue
353
 
      if(nuu.eq.0) go to 800
354
 
c  we determine the matrix (avv) and then we reduce her to an unit
355
 
c  upper triangular form (rv) using givens rotations without square
356
 
c  roots. we apply the same transformations to the columns of matrix
357
 
c  g to obtain the (nv-7) x (nu-6-iop0-iop1) matrix h.
358
 
c  we store matrix (rv) into av1 and av2, h into c.
359
 
c  the nv7 x nv7 triangular unit upper matrix (rv) has the form
360
 
c              | av1 '     |
361
 
c       (rv) = |     ' av2 |
362
 
c              |  0  '     |
363
 
c  with (av2) a nv7 x 4 matrix and (av1) a nv11 x nv11 unit upper
364
 
c  triangular matrix of bandwidth 5.
365
 
      ncof = nuu*nv7
366
 
c  initialization.
367
 
      do 430 i=1,ncof
368
 
        c(i) = 0.
369
 
 430  continue
370
 
      do 440 i=1,nv4
371
 
        av1(i,5) = 0.
372
 
        do 440 j=1,4
373
 
          av1(i,j) = 0.
374
 
          av2(i,j) = 0.
375
 
 440  continue
376
 
      jper = 0
377
 
      nrold = 0
378
 
      do 770 it=1,mv
379
 
        number = nrv(it)
380
 
 450    if(nrold.eq.number) go to 480
381
 
        if(p.le.0.) go to 760
382
 
c  fetch a new row of matrix (bv).
383
 
        n1 = nrold+1
384
 
        do 460 j=1,5
385
 
          h(j) = bv(n1,j)*pinv
386
 
 460    continue
387
 
c  find the appropiate row of g.
388
 
        do 465 j=1,nuu
389
 
          right(j) = 0.
390
 
 465    continue
391
 
        if(mv.eq.mvv) go to 510
392
 
        l = mv+n1
393
 
        do 470 j=1,nuu
394
 
          right(j) = q(l)
395
 
          l = l+mvv
396
 
 470    continue
397
 
        go to 510
398
 
c  fetch a new row of matrix (spv)
399
 
 480    h(5) = 0.
400
 
        do 490 j=1,4
401
 
          h(j) = spv(it,j)
402
 
 490    continue
403
 
c  find the appropiate row of g.
404
 
        l = it
405
 
        do 500 j=1,nuu
406
 
          right(j) = q(l)
407
 
          l = l+mvv
408
 
 500    continue
409
 
c  test whether there are non-zero values in the new row of (avv)
410
 
c  corresponding to the b-splines n(j;v),j=nv7+1,...,nv4.
411
 
 510     if(nrold.lt.nv11) go to 710
412
 
         if(jper.ne.0) go to 550
413
 
c  initialize the matrix (av2).
414
 
         jk = nv11+1
415
 
         do 540 i=1,4
416
 
            ik = jk
417
 
            do 520 j=1,5
418
 
               if(ik.le.0) go to 530
419
 
               av2(ik,i) = av1(ik,j)
420
 
               ik = ik-1
421
 
 520        continue
422
 
 530        jk = jk+1
423
 
 540     continue
424
 
         jper = 1
425
 
c  if one of the non-zero elements of the new row corresponds to one of
426
 
c  the b-splines n(j;v),j=nv7+1,...,nv4, we take account of condition
427
 
c  (2) for setting up this row of (avv). the row is stored in h1( the
428
 
c  part with respect to av1) and h2 (the part with respect to av2).
429
 
 550     do 560 i=1,4
430
 
            h1(i) = 0.
431
 
            h2(i) = 0.
432
 
 560     continue
433
 
         h1(5) = 0.
434
 
         j = nrold-nv11
435
 
         do 600 i=1,5
436
 
            j = j+1
437
 
            l0 = j
438
 
 570        l1 = l0-4
439
 
            if(l1.le.0) go to 590
440
 
            if(l1.le.nv11) go to 580
441
 
            l0 = l1-nv11
442
 
            go to 570
443
 
 580        h1(l1) = h(i)
444
 
            go to 600
445
 
 590        h2(l0) = h2(l0) + h(i)
446
 
 600     continue
447
 
c  rotate the new row of (avv) into triangle.
448
 
         if(nv11.le.0) go to 670
449
 
c  rotations with the rows 1,2,...,nv11 of (avv).
450
 
         do 660 j=1,nv11
451
 
            piv = h1(1)
452
 
            i2 = min0(nv11-j,4)
453
 
            if(piv.eq.0.) go to 640
454
 
c  calculate the parameters of the givens transformation.
455
 
            call fpgivs(piv,av1(j,1),co,si)
456
 
c  apply that transformation to the columns of matrix g.
457
 
            ic = j
458
 
            do 610 i=1,nuu
459
 
               call fprota(co,si,right(i),c(ic))
460
 
               ic = ic+nv7
461
 
 610        continue
462
 
c  apply that transformation to the rows of (avv) with respect to av2.
463
 
            do 620 i=1,4
464
 
               call fprota(co,si,h2(i),av2(j,i))
465
 
 620        continue
466
 
c  apply that transformation to the rows of (avv) with respect to av1.
467
 
            if(i2.eq.0) go to 670
468
 
            do 630 i=1,i2
469
 
               i1 = i+1
470
 
               call fprota(co,si,h1(i1),av1(j,i1))
471
 
 630        continue
472
 
 640        do 650 i=1,i2
473
 
               h1(i) = h1(i+1)
474
 
 650        continue
475
 
            h1(i2+1) = 0.
476
 
 660     continue
477
 
c  rotations with the rows nv11+1,...,nv7 of avv.
478
 
 670     do 700 j=1,4
479
 
            ij = nv11+j
480
 
            if(ij.le.0) go to 700
481
 
            piv = h2(j)
482
 
            if(piv.eq.0.) go to 700
483
 
c  calculate the parameters of the givens transformation.
484
 
            call fpgivs(piv,av2(ij,j),co,si)
485
 
c  apply that transformation to the columns of matrix g.
486
 
            ic = ij
487
 
            do 680 i=1,nuu
488
 
               call fprota(co,si,right(i),c(ic))
489
 
               ic = ic+nv7
490
 
 680        continue
491
 
            if(j.eq.4) go to 700
492
 
c  apply that transformation to the rows of (avv) with respect to av2.
493
 
            j1 = j+1
494
 
            do 690 i=j1,4
495
 
               call fprota(co,si,h2(i),av2(ij,i))
496
 
 690        continue
497
 
 700     continue
498
 
c  we update the sum of squared residuals.
499
 
         do 705 i=1,nuu
500
 
           sq = sq+right(i)**2
501
 
 705     continue
502
 
         go to 750
503
 
c  rotation into triangle of the new row of (avv), in case the elements
504
 
c  corresponding to the b-splines n(j;v),j=nv7+1,...,nv4 are all zero.
505
 
 710     irot =nrold
506
 
         do 740 i=1,5
507
 
            irot = irot+1
508
 
            piv = h(i)
509
 
            if(piv.eq.0.) go to 740
510
 
c  calculate the parameters of the givens transformation.
511
 
            call fpgivs(piv,av1(irot,1),co,si)
512
 
c  apply that transformation to the columns of matrix g.
513
 
            ic = irot
514
 
            do 720 j=1,nuu
515
 
               call fprota(co,si,right(j),c(ic))
516
 
               ic = ic+nv7
517
 
 720        continue
518
 
c  apply that transformation to the rows of (avv).
519
 
            if(i.eq.5) go to 740
520
 
            i2 = 1
521
 
            i3 = i+1
522
 
            do 730 j=i3,5
523
 
               i2 = i2+1
524
 
               call fprota(co,si,h(j),av1(irot,i2))
525
 
 730        continue
526
 
 740     continue
527
 
c  we update the sum of squared residuals.
528
 
         do 745 i=1,nuu
529
 
           sq = sq+right(i)**2
530
 
 745     continue
531
 
 750     if(nrold.eq.number) go to 770
532
 
 760     nrold = nrold+1
533
 
         go to 450
534
 
 770  continue
535
 
c  test whether the b-spline coefficients must be determined.
536
 
      if(iback.ne.0) return
537
 
c  backward substitution to obtain the b-spline coefficients as the
538
 
c  solution of the linear system    (rv) (cr) (ru)' = h.
539
 
c  first step: solve the system  (rv) (c1) = h.
540
 
      k = 1
541
 
      do 780 i=1,nuu
542
 
         call fpbacp(av1,av2,c(k),nv7,4,c(k),5,nv)
543
 
         k = k+nv7
544
 
 780  continue
545
 
c  second step: solve the system  (cr) (ru)' = (c1).
546
 
      k = 0
547
 
      do 795 j=1,nv7
548
 
        k = k+1
549
 
        l = k
550
 
        do 785 i=1,nuu
551
 
          right(i) = c(l)
552
 
          l = l+nv7
553
 
 785    continue
554
 
        call fpback(au,right,nuu,5,right,nu)
555
 
        l = k
556
 
        do 790 i=1,nuu
557
 
           c(l) = right(i)
558
 
           l = l+nv7
559
 
 790    continue
560
 
 795  continue
561
 
c  calculate from the conditions (2)-(3)-(4), the remaining b-spline
562
 
c  coefficients.
563
 
 800  ncof = nu4*nv4
564
 
      j = ncof
565
 
      do 805 l=1,nv4
566
 
         q(l) = dr01
567
 
         q(j) = dr11
568
 
         j = j-1
569
 
 805  continue
570
 
      i = nv4
571
 
      j = 0
572
 
      if(iop0.eq.0) go to 815
573
 
      do 810 l=1,nv4
574
 
         i = i+1
575
 
         q(i) = c0(l)
576
 
 810  continue
577
 
 815  if(nuu.eq.0) go to 835
578
 
      do 830 l=1,nuu
579
 
         ii = i
580
 
         do 820 k=1,nv7
581
 
            i = i+1
582
 
            j = j+1
583
 
            q(i) = c(j)
584
 
 820     continue
585
 
         do 825 k=1,3
586
 
            ii = ii+1
587
 
            i = i+1
588
 
            q(i) = q(ii)
589
 
 825     continue
590
 
 830  continue
591
 
 835  if(iop1.eq.0) go to 845
592
 
      do 840 l=1,nv4
593
 
         i = i+1
594
 
         q(i) = c1(l)
595
 
 840  continue
596
 
 845  do 850 i=1,ncof
597
 
         c(i) = q(i)
598
 
 850  continue
599
 
c  calculate the quantities
600
 
c    res(i,j) = (r(i,j) - s(u(i),v(j)))**2 , i=1,2,..,mu;j=1,2,..,mv
601
 
c    fp = sumi=1,mu(sumj=1,mv(res(i,j)))
602
 
c    fpu(r) = sum''i(sumj=1,mv(res(i,j))) , r=1,2,...,nu-7
603
 
c                  tu(r+3) <= u(i) <= tu(r+4)
604
 
c    fpv(r) = sumi=1,mu(sum''j(res(i,j))) , r=1,2,...,nv-7
605
 
c                  tv(r+3) <= v(j) <= tv(r+4)
606
 
      fp = 0.
607
 
      do 890 i=1,nu
608
 
        fpu(i) = 0.
609
 
 890  continue
610
 
      do 900 i=1,nv
611
 
        fpv(i) = 0.
612
 
 900  continue
613
 
      ir = 0
614
 
      nroldu = 0
615
 
c  main loop for the different grid points.
616
 
      do 950 i1=1,mu
617
 
        numu = nru(i1)
618
 
        numu1 = numu+1
619
 
        nroldv = 0
620
 
        do 940 i2=1,mv
621
 
          numv = nrv(i2)
622
 
          numv1 = numv+1
623
 
          ir = ir+1
624
 
c  evaluate s(u,v) at the current grid point by making the sum of the
625
 
c  cross products of the non-zero b-splines at (u,v), multiplied with
626
 
c  the appropiate b-spline coefficients.
627
 
          term = 0.
628
 
          k1 = numu*nv4+numv
629
 
          do 920 l1=1,4
630
 
            k2 = k1
631
 
            fac = spu(i1,l1)
632
 
            do 910 l2=1,4
633
 
              k2 = k2+1
634
 
              term = term+fac*spv(i2,l2)*c(k2)
635
 
 910        continue
636
 
            k1 = k1+nv4
637
 
 920      continue
638
 
c  calculate the squared residual at the current grid point.
639
 
          term = (r(ir)-term)**2
640
 
c  adjust the different parameters.
641
 
          fp = fp+term
642
 
          fpu(numu1) = fpu(numu1)+term
643
 
          fpv(numv1) = fpv(numv1)+term
644
 
          fac = term*half
645
 
          if(numv.eq.nroldv) go to 930
646
 
          fpv(numv1) = fpv(numv1)-fac
647
 
          fpv(numv) = fpv(numv)+fac
648
 
 930      nroldv = numv
649
 
          if(numu.eq.nroldu) go to 940
650
 
          fpu(numu1) = fpu(numu1)-fac
651
 
          fpu(numu) = fpu(numu)+fac
652
 
 940    continue
653
 
        nroldu = numu
654
 
 950  continue
655
 
      return
656
 
      end