subroutine fpgrsp(ifsu,ifsv,ifbu,ifbv,iback,u,mu,v,mv,r,mr,dr, * iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,spu,spv, * right,q,au,av1,av2,bu,bv,a0,a1,b0,b1,c0,c1,cosi,nru,nrv) c .. c ..scalar arguments.. real*8 p,sq,fp integer ifsu,ifsv,ifbu,ifbv,iback,mu,mv,mr,iop0,iop1,nu,nv,nc, * mm,mvnu c ..array arguments.. real*8 u(mu),v(mv),r(mr),dr(6),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv) *, * spu(mu,4),spv(mv,4),right(mm),q(mvnu),au(nu,5),av1(nv,6),c0(nv), * av2(nv,4),a0(2,mv),b0(2,nv),cosi(2,nv),bu(nu,5),bv(nv,5),c1(nv), * a1(2,mv),b1(2,nv) integer nru(mu),nrv(mv) c ..local scalars.. real*8 arg,co,dr01,dr02,dr03,dr11,dr12,dr13,fac,fac0,fac1,pinv,piv *, * si,term,one,three,half integer i,ic,ii,ij,ik,iq,irot,it,ir,i0,i1,i2,i3,j,jj,jk,jper, * j0,j1,k,k1,k2,l,l0,l1,l2,mvv,ncof,nrold,nroldu,nroldv,number, * numu,numu1,numv,numv1,nuu,nu4,nu7,nu8,nu9,nv11,nv4,nv7,nv8,n1 c ..local arrays.. real*8 h(5),h1(5),h2(4) c ..function references.. integer min0 real*8 cos,sin c ..subroutine references.. c fpback,fpbspl,fpgivs,fpcyt1,fpcyt2,fpdisc,fpbacp,fprota c .. c let c | (spu) | | (spv) | c (au) = | -------------- | (av) = | -------------- | c | sqrt(1/p) (bu) | | sqrt(1/p) (bv) | c c | r ' 0 | c q = | ------ | c | 0 ' 0 | c c with c : the (nu-4) x (nv-4) matrix which contains the b-spline c coefficients. c r : the mu x mv matrix which contains the function values. c spu,spv: the mu x (nu-4), resp. mv x (nv-4) observation matrices c according to the least-squares problems in the u-,resp. c v-direction. c bu,bv : the (nu-7) x (nu-4),resp. (nv-7) x (nv-4) matrices c containing the discontinuity jumps of the derivatives c of the b-splines in the u-,resp.v-variable at the knots c the b-spline coefficients of the smoothing spline are then calculated c as the least-squares solution of the following over-determined linear c system of equations c c (1) (av) c (au)' = q c c subject to the constraints c c (2) c(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4 c c (3) if iop0 = 0 c(1,j) = dr(1) c iop0 = 1 c(1,j) = dr(1) c c(2,j) = dr(1)+(dr(2)*cosi(1,j)+dr(3)*cosi(2,j))* c tu(5)/3. = c0(j) , j=1,2,...nv-4 c c (4) if iop1 = 0 c(nu-4,j) = dr(4) c iop1 = 1 c(nu-4,j) = dr(4) c c(nu-5,j) = dr(4)+(dr(5)*cosi(1,j)+dr(6)*cosi(2,j)) c *(tu(nu-4)-tu(nu-3))/3. = c1(j) c c set constants one = 1 three = 3 half = 0.5 c initialization nu4 = nu-4 nu7 = nu-7 nu8 = nu-8 nu9 = nu-9 nv4 = nv-4 nv7 = nv-7 nv8 = nv-8 nv11 = nv-11 nuu = nu4-iop0-iop1-2 if(p.gt.0.) pinv = one/p c it depends on the value of the flags ifsu,ifsv,ifbu,ifbv,iop0,iop1 c and on the value of p whether the matrices (spu), (spv), (bu), (bv), c (cosi) still must be determined. if(ifsu.ne.0) go to 30 c calculate the non-zero elements of the matrix (spu) which is the ob- c servation matrix according to the least-squares spline approximation c problem in the u-direction. l = 4 l1 = 5 number = 0 do 25 it=1,mu arg = u(it) 10 if(arg.lt.tu(l1) .or. l.eq.nu4) go to 15 l = l1 l1 = l+1 number = number+1 go to 10 15 call fpbspl(tu,nu,3,arg,l,h) do 20 i=1,4 spu(it,i) = h(i) 20 continue nru(it) = number 25 continue ifsu = 1 c calculate the non-zero elements of the matrix (spv) which is the ob- c servation matrix according to the least-squares spline approximation c problem in the v-direction. 30 if(ifsv.ne.0) go to 85 l = 4 l1 = 5 number = 0 do 50 it=1,mv arg = v(it) 35 if(arg.lt.tv(l1) .or. l.eq.nv4) go to 40 l = l1 l1 = l+1 number = number+1 go to 35 40 call fpbspl(tv,nv,3,arg,l,h) do 45 i=1,4 spv(it,i) = h(i) 45 continue nrv(it) = number 50 continue ifsv = 1 if(iop0.eq.0 .and. iop1.eq.0) go to 85 c calculate the coefficients of the interpolating splines for cos(v) c and sin(v). do 55 i=1,nv4 cosi(1,i) = 0. cosi(2,i) = 0. 55 continue if(nv7.lt.4) go to 85 do 65 i=1,nv7 l = i+3 arg = tv(l) call fpbspl(tv,nv,3,arg,l,h) do 60 j=1,3 av1(i,j) = h(j) 60 continue cosi(1,i) = cos(arg) cosi(2,i) = sin(arg) 65 continue call fpcyt1(av1,nv7,nv) do 80 j=1,2 do 70 i=1,nv7 right(i) = cosi(j,i) 70 continue call fpcyt2(av1,nv7,right,right,nv) do 75 i=1,nv7 cosi(j,i+1) = right(i) 75 continue cosi(j,1) = cosi(j,nv7+1) cosi(j,nv7+2) = cosi(j,2) cosi(j,nv4) = cosi(j,3) 80 continue 85 if(p.le.0.) go to 150 c calculate the non-zero elements of the matrix (bu). if(ifbu.ne.0 .or. nu8.eq.0) go to 90 call fpdisc(tu,nu,5,bu,nu) ifbu = 1 c calculate the non-zero elements of the matrix (bv). 90 if(ifbv.ne.0 .or. nv8.eq.0) go to 150 call fpdisc(tv,nv,5,bv,nv) ifbv = 1 c substituting (2),(3) and (4) into (1), we obtain the overdetermined c system c (5) (avv) (cc) (auu)' = (qq) c from which the nuu*nv7 remaining coefficients c c(i,j) , i=2+iop0,3+iop0,...,nu-5-iop1,j=1,2,...,nv-7. c the elements of (cc), are then determined in the least-squares sense. c simultaneously, we compute the resulting sum of squared residuals sq. 150 dr01 = dr(1) dr11 = dr(4) do 155 i=1,mv a0(1,i) = dr01 a1(1,i) = dr11 155 continue if(nv8.eq.0 .or. p.le.0.) go to 165 do 160 i=1,nv8 b0(1,i) = 0. b1(1,i) = 0. 160 continue 165 mvv = mv if(iop0.eq.0) go to 195 fac = (tu(5)-tu(4))/three dr02 = dr(2)*fac dr03 = dr(3)*fac do 170 i=1,nv4 c0(i) = dr01+dr02*cosi(1,i)+dr03*cosi(2,i) 170 continue do 180 i=1,mv number = nrv(i) fac = 0. do 175 j=1,4 number = number+1 fac = fac+c0(number)*spv(i,j) 175 continue a0(2,i) = fac 180 continue if(nv8.eq.0 .or. p.le.0.) go to 195 do 190 i=1,nv8 number = i fac = 0. do 185 j=1,5 fac = fac+c0(number)*bv(i,j) number = number+1 185 continue b0(2,i) = fac*pinv 190 continue mvv = mv+nv8 195 if(iop1.eq.0) go to 225 fac = (tu(nu4)-tu(nu4+1))/three dr12 = dr(5)*fac dr13 = dr(6)*fac do 200 i=1,nv4 c1(i) = dr11+dr12*cosi(1,i)+dr13*cosi(2,i) 200 continue do 210 i=1,mv number = nrv(i) fac = 0. do 205 j=1,4 number = number+1 fac = fac+c1(number)*spv(i,j) 205 continue a1(2,i) = fac 210 continue if(nv8.eq.0 .or. p.le.0.) go to 225 do 220 i=1,nv8 number = i fac = 0. do 215 j=1,5 fac = fac+c1(number)*bv(i,j) number = number+1 215 continue b1(2,i) = fac*pinv 220 continue mvv = mv+nv8 c we first determine the matrices (auu) and (qq). then we reduce the c matrix (auu) to an unit upper triangular form (ru) using givens c rotations without square roots. we apply the same transformations to c the rows of matrix qq to obtain the mv x nuu matrix g. c we store matrix (ru) into au and g into q. 225 l = mvv*nuu c initialization. sq = 0. if(l.eq.0) go to 245 do 230 i=1,l q(i) = 0. 230 continue do 240 i=1,nuu do 240 j=1,5 au(i,j) = 0. 240 continue l = 0 245 nrold = 0 n1 = nrold+1 do 420 it=1,mu number = nru(it) c find the appropriate column of q. 250 do 260 j=1,mvv right(j) = 0. 260 continue if(nrold.eq.number) go to 280 if(p.le.0.) go to 410 c fetch a new row of matrix (bu). do 270 j=1,5 h(j) = bu(n1,j)*pinv 270 continue i0 = 1 i1 = 5 go to 310 c fetch a new row of matrix (spu). 280 do 290 j=1,4 h(j) = spu(it,j) 290 continue c find the appropriate column of q. do 300 j=1,mv l = l+1 right(j) = r(l) 300 continue i0 = 1 i1 = 4 310 j0 = n1 j1 = nu7-number c take into account that we eliminate the constraints (3) 315 if(j0-1.gt.iop0) go to 335 fac0 = h(i0) do 320 j=1,mv right(j) = right(j)-fac0*a0(j0,j) 320 continue if(mv.eq.mvv) go to 330 j = mv do 325 jj=1,nv8 j = j+1 right(j) = right(j)-fac0*b0(j0,jj) 325 continue 330 j0 = j0+1 i0 = i0+1 go to 315 c take into account that we eliminate the constraints (4) 335 if(j1-1.gt.iop1) go to 360 fac1 = h(i1) do 340 j=1,mv right(j) = right(j)-fac1*a1(j1,j) 340 continue if(mv.eq.mvv) go to 350 j = mv do 345 jj=1,nv8 j = j+1 right(j) = right(j)-fac1*b1(j1,jj) 345 continue 350 j1 = j1+1 i1 = i1-1 go to 335 360 irot = nrold-iop0-1 if(irot.lt.0) irot = 0 c rotate the new row of matrix (auu) into triangle. if(i0.gt.i1) go to 390 do 385 i=i0,i1 irot = irot+1 piv = h(i) if(piv.eq.0.) go to 385 c calculate the parameters of the givens transformation. call fpgivs(piv,au(irot,1),co,si) c apply that transformation to the rows of matrix (qq). iq = (irot-1)*mvv do 370 j=1,mvv iq = iq+1 call fprota(co,si,right(j),q(iq)) 370 continue c apply that transformation to the columns of (auu). if(i.eq.i1) go to 385 i2 = 1 i3 = i+1 do 380 j=i3,i1 i2 = i2+1 call fprota(co,si,h(j),au(irot,i2)) 380 continue 385 continue c we update the sum of squared residuals. 390 do 395 j=1,mvv sq = sq+right(j)**2 395 continue if(nrold.eq.number) go to 420 410 nrold = n1 n1 = n1+1 go to 250 420 continue if(nuu.eq.0) go to 800 c we determine the matrix (avv) and then we reduce her to an unit c upper triangular form (rv) using givens rotations without square c roots. we apply the same transformations to the columns of matrix c g to obtain the (nv-7) x (nu-6-iop0-iop1) matrix h. c we store matrix (rv) into av1 and av2, h into c. c the nv7 x nv7 triangular unit upper matrix (rv) has the form c | av1 ' | c (rv) = | ' av2 | c | 0 ' | c with (av2) a nv7 x 4 matrix and (av1) a nv11 x nv11 unit upper c triangular matrix of bandwidth 5. ncof = nuu*nv7 c initialization. do 430 i=1,ncof c(i) = 0. 430 continue do 440 i=1,nv4 av1(i,5) = 0. do 440 j=1,4 av1(i,j) = 0. av2(i,j) = 0. 440 continue jper = 0 nrold = 0 do 770 it=1,mv number = nrv(it) 450 if(nrold.eq.number) go to 480 if(p.le.0.) go to 760 c fetch a new row of matrix (bv). n1 = nrold+1 do 460 j=1,5 h(j) = bv(n1,j)*pinv 460 continue c find the appropiate row of g. do 465 j=1,nuu right(j) = 0. 465 continue if(mv.eq.mvv) go to 510 l = mv+n1 do 470 j=1,nuu right(j) = q(l) l = l+mvv 470 continue go to 510 c fetch a new row of matrix (spv) 480 h(5) = 0. do 490 j=1,4 h(j) = spv(it,j) 490 continue c find the appropiate row of g. l = it do 500 j=1,nuu right(j) = q(l) l = l+mvv 500 continue c test whether there are non-zero values in the new row of (avv) c corresponding to the b-splines n(j;v),j=nv7+1,...,nv4. 510 if(nrold.lt.nv11) go to 710 if(jper.ne.0) go to 550 c initialize the matrix (av2). jk = nv11+1 do 540 i=1,4 ik = jk do 520 j=1,5 if(ik.le.0) go to 530 av2(ik,i) = av1(ik,j) ik = ik-1 520 continue 530 jk = jk+1 540 continue jper = 1 c if one of the non-zero elements of the new row corresponds to one of c the b-splines n(j;v),j=nv7+1,...,nv4, we take account of condition c (2) for setting up this row of (avv). the row is stored in h1( the c part with respect to av1) and h2 (the part with respect to av2). 550 do 560 i=1,4 h1(i) = 0. h2(i) = 0. 560 continue h1(5) = 0. j = nrold-nv11 do 600 i=1,5 j = j+1 l0 = j 570 l1 = l0-4 if(l1.le.0) go to 590 if(l1.le.nv11) go to 580 l0 = l1-nv11 go to 570 580 h1(l1) = h(i) go to 600 590 h2(l0) = h2(l0) + h(i) 600 continue c rotate the new row of (avv) into triangle. if(nv11.le.0) go to 670 c rotations with the rows 1,2,...,nv11 of (avv). do 660 j=1,nv11 piv = h1(1) i2 = min0(nv11-j,4) if(piv.eq.0.) go to 640 c calculate the parameters of the givens transformation. call fpgivs(piv,av1(j,1),co,si) c apply that transformation to the columns of matrix g. ic = j do 610 i=1,nuu call fprota(co,si,right(i),c(ic)) ic = ic+nv7 610 continue c apply that transformation to the rows of (avv) with respect to av2. do 620 i=1,4 call fprota(co,si,h2(i),av2(j,i)) 620 continue c apply that transformation to the rows of (avv) with respect to av1. if(i2.eq.0) go to 670 do 630 i=1,i2 i1 = i+1 call fprota(co,si,h1(i1),av1(j,i1)) 630 continue 640 do 650 i=1,i2 h1(i) = h1(i+1) 650 continue h1(i2+1) = 0. 660 continue c rotations with the rows nv11+1,...,nv7 of avv. 670 do 700 j=1,4 ij = nv11+j if(ij.le.0) go to 700 piv = h2(j) if(piv.eq.0.) go to 700 c calculate the parameters of the givens transformation. call fpgivs(piv,av2(ij,j),co,si) c apply that transformation to the columns of matrix g. ic = ij do 680 i=1,nuu call fprota(co,si,right(i),c(ic)) ic = ic+nv7 680 continue if(j.eq.4) go to 700 c apply that transformation to the rows of (avv) with respect to av2. j1 = j+1 do 690 i=j1,4 call fprota(co,si,h2(i),av2(ij,i)) 690 continue 700 continue c we update the sum of squared residuals. do 705 i=1,nuu sq = sq+right(i)**2 705 continue go to 750 c rotation into triangle of the new row of (avv), in case the elements c corresponding to the b-splines n(j;v),j=nv7+1,...,nv4 are all zero. 710 irot =nrold do 740 i=1,5 irot = irot+1 piv = h(i) if(piv.eq.0.) go to 740 c calculate the parameters of the givens transformation. call fpgivs(piv,av1(irot,1),co,si) c apply that transformation to the columns of matrix g. ic = irot do 720 j=1,nuu call fprota(co,si,right(j),c(ic)) ic = ic+nv7 720 continue c apply that transformation to the rows of (avv). if(i.eq.5) go to 740 i2 = 1 i3 = i+1 do 730 j=i3,5 i2 = i2+1 call fprota(co,si,h(j),av1(irot,i2)) 730 continue 740 continue c we update the sum of squared residuals. do 745 i=1,nuu sq = sq+right(i)**2 745 continue 750 if(nrold.eq.number) go to 770 760 nrold = nrold+1 go to 450 770 continue c test whether the b-spline coefficients must be determined. if(iback.ne.0) return c backward substitution to obtain the b-spline coefficients as the c solution of the linear system (rv) (cr) (ru)' = h. c first step: solve the system (rv) (c1) = h. k = 1 do 780 i=1,nuu call fpbacp(av1,av2,c(k),nv7,4,c(k),5,nv) k = k+nv7 780 continue c second step: solve the system (cr) (ru)' = (c1). k = 0 do 795 j=1,nv7 k = k+1 l = k do 785 i=1,nuu right(i) = c(l) l = l+nv7 785 continue call fpback(au,right,nuu,5,right,nu) l = k do 790 i=1,nuu c(l) = right(i) l = l+nv7 790 continue 795 continue c calculate from the conditions (2)-(3)-(4), the remaining b-spline c coefficients. 800 ncof = nu4*nv4 j = ncof do 805 l=1,nv4 q(l) = dr01 q(j) = dr11 j = j-1 805 continue i = nv4 j = 0 if(iop0.eq.0) go to 815 do 810 l=1,nv4 i = i+1 q(i) = c0(l) 810 continue 815 if(nuu.eq.0) go to 835 do 830 l=1,nuu ii = i do 820 k=1,nv7 i = i+1 j = j+1 q(i) = c(j) 820 continue do 825 k=1,3 ii = ii+1 i = i+1 q(i) = q(ii) 825 continue 830 continue 835 if(iop1.eq.0) go to 845 do 840 l=1,nv4 i = i+1 q(i) = c1(l) 840 continue 845 do 850 i=1,ncof c(i) = q(i) 850 continue c calculate the quantities c res(i,j) = (r(i,j) - s(u(i),v(j)))**2 , i=1,2,..,mu;j=1,2,..,mv c fp = sumi=1,mu(sumj=1,mv(res(i,j))) c fpu(r) = sum''i(sumj=1,mv(res(i,j))) , r=1,2,...,nu-7 c tu(r+3) <= u(i) <= tu(r+4) c fpv(r) = sumi=1,mu(sum''j(res(i,j))) , r=1,2,...,nv-7 c tv(r+3) <= v(j) <= tv(r+4) fp = 0. do 890 i=1,nu fpu(i) = 0. 890 continue do 900 i=1,nv fpv(i) = 0. 900 continue ir = 0 nroldu = 0 c main loop for the different grid points. do 950 i1=1,mu numu = nru(i1) numu1 = numu+1 nroldv = 0 do 940 i2=1,mv numv = nrv(i2) numv1 = numv+1 ir = ir+1 c evaluate s(u,v) at the current grid point by making the sum of the c cross products of the non-zero b-splines at (u,v), multiplied with c the appropiate b-spline coefficients. term = 0. k1 = numu*nv4+numv do 920 l1=1,4 k2 = k1 fac = spu(i1,l1) do 910 l2=1,4 k2 = k2+1 term = term+fac*spv(i2,l2)*c(k2) 910 continue k1 = k1+nv4 920 continue c calculate the squared residual at the current grid point. term = (r(ir)-term)**2 c adjust the different parameters. fp = fp+term fpu(numu1) = fpu(numu1)+term fpv(numv1) = fpv(numv1)+term fac = term*half if(numv.eq.nroldv) go to 930 fpv(numv1) = fpv(numv1)-fac fpv(numv) = fpv(numv)+fac 930 nroldv = numv if(numu.eq.nroldu) go to 940 fpu(numu1) = fpu(numu1)-fac fpu(numu) = fpu(numu)+fac 940 continue nroldu = numu 950 continue return end