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)
7
integer ifsu,ifsv,ifbu,ifbv,iback,mu,mv,mr,iop0,iop1,nu,nv,nc,
10
real*8 u(mu),v(mv),r(mr),dr(6),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv)
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),
15
integer nru(mu),nrv(mv)
17
real*8 arg,co,dr01,dr02,dr03,dr11,dr12,dr13,fac,fac0,fac1,pinv,piv
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
24
real*8 h(5),h1(5),h2(4)
25
c ..function references..
28
c ..subroutine references..
29
c fpback,fpbspl,fpgivs,fpcyt1,fpcyt2,fpdisc,fpbacp,fprota
33
c (au) = | -------------- | (av) = | -------------- |
34
c | sqrt(1/p) (bu) | | sqrt(1/p) (bv) |
40
c with c : the (nu-4) x (nv-4) matrix which contains the b-spline
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.
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
53
c (1) (av) c (au)' = q
55
c subject to the constraints
57
c (2) c(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4
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
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)
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.
96
10 if(arg.lt.tu(l1) .or. l.eq.nu4) go to 15
101
15 call fpbspl(tu,nu,3,arg,l,h)
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
117
35 if(arg.lt.tv(l1) .or. l.eq.nv4) go to 40
122
40 call fpbspl(tv,nv,3,arg,l,h)
129
if(iop0.eq.0 .and. iop1.eq.0) go to 85
130
c calculate the coefficients of the interpolating splines for cos(v)
136
if(nv7.lt.4) go to 85
140
call fpbspl(tv,nv,3,arg,l,h)
147
call fpcyt1(av1,nv7,nv)
152
call fpcyt2(av1,nv7,right,right,nv)
154
cosi(j,i+1) = right(i)
156
cosi(j,1) = cosi(j,nv7+1)
157
cosi(j,nv7+2) = cosi(j,2)
158
cosi(j,nv4) = cosi(j,3)
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)
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)
169
c substituting (2),(3) and (4) into (1), we obtain the overdetermined
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.
182
if(nv8.eq.0 .or. p.le.0.) go to 165
188
if(iop0.eq.0) go to 195
189
fac = (tu(5)-tu(4))/three
193
c0(i) = dr01+dr02*cosi(1,i)+dr03*cosi(2,i)
200
fac = fac+c0(number)*spv(i,j)
204
if(nv8.eq.0 .or. p.le.0.) go to 195
209
fac = fac+c0(number)*bv(i,j)
215
195 if(iop1.eq.0) go to 225
216
fac = (tu(nu4)-tu(nu4+1))/three
220
c1(i) = dr11+dr12*cosi(1,i)+dr13*cosi(2,i)
227
fac = fac+c1(number)*spv(i,j)
231
if(nv8.eq.0 .or. p.le.0.) go to 225
236
fac = fac+c1(number)*bv(i,j)
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.
263
c find the appropriate column of q.
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).
276
c fetch a new row of matrix (spu).
280
c find the appropriate column of q.
289
c take into account that we eliminate the constraints (3)
290
315 if(j0-1.gt.iop0) go to 335
293
right(j) = right(j)-fac0*a0(j0,j)
295
if(mv.eq.mvv) go to 330
299
right(j) = right(j)-fac0*b0(j0,jj)
304
c take into account that we eliminate the constraints (4)
305
335 if(j1-1.gt.iop1) go to 360
308
right(j) = right(j)-fac1*a1(j1,j)
310
if(mv.eq.mvv) go to 350
314
right(j) = right(j)-fac1*b1(j1,jj)
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
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).
333
call fprota(co,si,right(j),q(iq))
335
c apply that transformation to the columns of (auu).
336
if(i.eq.i1) go to 385
341
call fprota(co,si,h(j),au(irot,i2))
344
c we update the sum of squared residuals.
348
400 if(nrold.eq.number) go to 420
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
363
c with (av2) a nv7 x 4 matrix and (av1) a nv11 x nv11 unit upper
364
c triangular matrix of bandwidth 5.
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).
387
c find the appropiate row of g.
391
if(mv.eq.mvv) go to 510
398
c fetch a new row of matrix (spv)
403
c find the appropiate row of g.
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).
418
if(ik.le.0) go to 530
419
av2(ik,i) = av1(ik,j)
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).
439
if(l1.le.0) go to 590
440
if(l1.le.nv11) go to 580
445
590 h2(l0) = h2(l0) + h(i)
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).
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.
459
call fprota(co,si,right(i),c(ic))
462
c apply that transformation to the rows of (avv) with respect to av2.
464
call fprota(co,si,h2(i),av2(j,i))
466
c apply that transformation to the rows of (avv) with respect to av1.
467
if(i2.eq.0) go to 670
470
call fprota(co,si,h1(i1),av1(j,i1))
477
c rotations with the rows nv11+1,...,nv7 of avv.
480
if(ij.le.0) go to 700
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.
488
call fprota(co,si,right(i),c(ic))
492
c apply that transformation to the rows of (avv) with respect to av2.
495
call fprota(co,si,h2(i),av2(ij,i))
498
c we update the sum of squared residuals.
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.
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.
515
call fprota(co,si,right(j),c(ic))
518
c apply that transformation to the rows of (avv).
524
call fprota(co,si,h(j),av1(irot,i2))
527
c we update the sum of squared residuals.
531
750 if(nrold.eq.number) go to 770
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.
542
call fpbacp(av1,av2,c(k),nv7,4,c(k),5,nv)
545
c second step: solve the system (cr) (ru)' = (c1).
554
call fpback(au,right,nuu,5,right,nu)
561
c calculate from the conditions (2)-(3)-(4), the remaining b-spline
572
if(iop0.eq.0) go to 815
577
815 if(nuu.eq.0) go to 835
591
835 if(iop1.eq.0) go to 845
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)
615
c main loop for the different grid points.
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.
634
term = term+fac*spv(i2,l2)*c(k2)
638
c calculate the squared residual at the current grid point.
639
term = (r(ir)-term)**2
640
c adjust the different parameters.
642
fpu(numu1) = fpu(numu1)+term
643
fpv(numv1) = fpv(numv1)+term
645
if(numv.eq.nroldv) go to 930
646
fpv(numv1) = fpv(numv1)-fac
647
fpv(numv) = fpv(numv)+fac
649
if(numu.eq.nroldu) go to 940
650
fpu(numu1) = fpu(numu1)-fac
651
fpu(numu) = fpu(numu)+fac