subroutine fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s, * nxest,nyest,tol,maxit,nc,nx,tx,ny,ty,c,fp,fp0,fpold,reducx, * reducy,fpintx,fpinty,lastdi,nplusx,nplusy,nrx,nry,nrdatx,nrdaty, * wrk,lwrk,ier) c .. c ..scalar arguments.. real xb,xe,yb,ye,s,tol,fp,fp0,fpold,reducx,reducy integer iopt,mx,my,mz,kx,ky,nxest,nyest,maxit,nc,nx,ny,lastdi, * nplusx,nplusy,lwrk,ier c ..array arguments.. real x(mx),y(my),z(mz),tx(nxest),ty(nyest),c(nc),fpintx(nxest), * fpinty(nyest),wrk(lwrk) integer nrdatx(nxest),nrdaty(nyest),nrx(mx),nry(my) c ..local scalars real acc,fpms,f1,f2,f3,p,p1,p2,p3,rn,one,half,con1,con9,con4 integer i,ich1,ich3,ifbx,ifby,ifsx,ifsy,iter,j,kx1,kx2,ky1,ky2, * k3,l,lax,lay,lbx,lby,lq,lri,lsx,lsy,mk1,mm,mpm,mynx,ncof, * nk1x,nk1y,nmaxx,nmaxy,nminx,nminy,nplx,nply,npl1,nrintx, * nrinty,nxe,nxk,nye c ..function references.. real abs,fprati integer max0,min0 c ..subroutine references.. c fpgrre,fpknot c .. c set constants one = 1 half = 0.5e0 con1 = 0.1e0 con9 = 0.9e0 con4 = 0.4e-01 c we partition the working space. kx1 = kx+1 ky1 = ky+1 kx2 = kx1+1 ky2 = ky1+1 lsx = 1 lsy = lsx+mx*kx1 lri = lsy+my*ky1 mm = max0(nxest,my) lq = lri+mm mynx = nxest*my lax = lq+mynx nxk = nxest*kx2 lbx = lax+nxk lay = lbx+nxk lby = lay+nyest*ky2 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c part 1: determination of the number of knots and their position. c c **************************************************************** c c given a set of knots we compute the least-squares spline sinf(x,y), c c and the corresponding sum of squared residuals fp=f(p=inf). c c if iopt=-1 sinf(x,y) is the requested approximation. c c if iopt=0 or iopt=1 we check whether we can accept the knots: c c if fp <=s we will continue with the current set of knots. c c if fp > s we will increase the number of knots and compute the c c corresponding least-squares spline until finally fp<=s. c c the initial choice of knots depends on the value of s and iopt. c c if s=0 we have spline interpolation; in that case the number of c c knots equals nmaxx = mx+kx+1 and nmaxy = my+ky+1. c c if s>0 and c c *iopt=0 we first compute the least-squares polynomial of degree c c kx in x and ky in y; nx=nminx=2*kx+2 and ny=nymin=2*ky+2. c c *iopt=1 we start with the knots found at the last call of the c c routine, except for the case that s > fp0; then we can compute c c the least-squares polynomial directly. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c determine the number of knots for polynomial approximation. nminx = 2*kx1 nminy = 2*ky1 if(iopt.lt.0) go to 120 c acc denotes the absolute tolerance for the root of f(p)=s. acc = tol*s c find nmaxx and nmaxy which denote the number of knots in x- and y- c direction in case of spline interpolation. nmaxx = mx+kx1 nmaxy = my+ky1 c find nxe and nye which denote the maximum number of knots c allowed in each direction nxe = min0(nmaxx,nxest) nye = min0(nmaxy,nyest) if(s.gt.0.) go to 100 c if s = 0, s(x,y) is an interpolating spline. nx = nmaxx ny = nmaxy c test whether the required storage space exceeds the available one. if(ny.gt.nyest .or. nx.gt.nxest) go to 420 c find the position of the interior knots in case of interpolation. c the knots in the x-direction. mk1 = mx-kx1 if(mk1.eq.0) go to 60 k3 = kx/2 i = kx1+1 j = k3+2 if(k3*2.eq.kx) go to 40 do 30 l=1,mk1 tx(i) = x(j) i = i+1 j = j+1 30 continue go to 60 40 do 50 l=1,mk1 tx(i) = (x(j)+x(j-1))*half i = i+1 j = j+1 50 continue c the knots in the y-direction. 60 mk1 = my-ky1 if(mk1.eq.0) go to 120 k3 = ky/2 i = ky1+1 j = k3+2 if(k3*2.eq.ky) go to 80 do 70 l=1,mk1 ty(i) = y(j) i = i+1 j = j+1 70 continue go to 120 80 do 90 l=1,mk1 ty(i) = (y(j)+y(j-1))*half i = i+1 j = j+1 90 continue go to 120 c if s > 0 our initial choice of knots depends on the value of iopt. 100 if(iopt.eq.0) go to 115 if(fp0.le.s) go to 115 c if iopt=1 and fp0 > s we start computing the least- squares spline c according to the set of knots found at the last call of the routine. c we determine the number of grid coordinates x(i) inside each knot c interval (tx(l),tx(l+1)). l = kx2 j = 1 nrdatx(1) = 0 mpm = mx-1 do 105 i=2,mpm nrdatx(j) = nrdatx(j)+1 if(x(i).lt.tx(l)) go to 105 nrdatx(j) = nrdatx(j)-1 l = l+1 j = j+1 nrdatx(j) = 0 105 continue c we determine the number of grid coordinates y(i) inside each knot c interval (ty(l),ty(l+1)). l = ky2 j = 1 nrdaty(1) = 0 mpm = my-1 do 110 i=2,mpm nrdaty(j) = nrdaty(j)+1 if(y(i).lt.ty(l)) go to 110 nrdaty(j) = nrdaty(j)-1 l = l+1 j = j+1 nrdaty(j) = 0 110 continue go to 120 c if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares c polynomial of degree kx in x and ky in y (which is a spline without c interior knots). 115 nx = nminx ny = nminy nrdatx(1) = mx-2 nrdaty(1) = my-2 lastdi = 0 nplusx = 0 nplusy = 0 fp0 = 0. fpold = 0. reducx = 0. reducy = 0. 120 mpm = mx+my ifsx = 0 ifsy = 0 ifbx = 0 ifby = 0 p = -one c main loop for the different sets of knots.mpm=mx+my is a save upper c bound for the number of trials. do 250 iter=1,mpm if(nx.eq.nminx .and. ny.eq.nminy) ier = -2 c find nrintx (nrinty) which is the number of knot intervals in the c x-direction (y-direction). nrintx = nx-nminx+1 nrinty = ny-nminy+1 c find ncof, the number of b-spline coefficients for the current set c of knots. nk1x = nx-kx1 nk1y = ny-ky1 ncof = nk1x*nk1y c find the position of the additional knots which are needed for the c b-spline representation of s(x,y). i = nx do 130 j=1,kx1 tx(j) = xb tx(i) = xe i = i-1 130 continue i = ny do 140 j=1,ky1 ty(j) = yb ty(i) = ye i = i-1 140 continue c find the least-squares spline sinf(x,y) and calculate for each knot c interval tx(j+kx)<=x<=tx(j+kx+1) (ty(j+ky)<=y<=ty(j+ky+1)) the sum c of squared residuals fpintx(j),j=1,2,...,nx-2*kx-1 (fpinty(j),j=1,2, c ...,ny-2*ky-1) for the data points having their absciss (ordinate)- c value belonging to that interval. c fp gives the total sum of squared residuals. call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, * ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), * wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), * nrx,nry) if(ier.eq.(-2)) fp0 = fp c test whether the least-squares spline is an acceptable solution. if(iopt.lt.0) go to 440 fpms = fp-s if(abs(fpms) .lt. acc) go to 440 c if f(p=inf) < s, we accept the choice of knots. if(fpms.lt.0.) go to 300 c if nx=nmaxx and ny=nmaxy, sinf(x,y) is an interpolating spline. if(nx.eq.nmaxx .and. ny.eq.nmaxy) go to 430 c increase the number of knots. c if nx=nxe and ny=nye we cannot further increase the number of knots c because of the storage capacity limitation. if(nx.eq.nxe .and. ny.eq.nye) go to 420 ier = 0 c adjust the parameter reducx or reducy according to the direction c in which the last added knots were located. if(lastdi) 150,170,160 150 reducx = fpold-fp go to 170 160 reducy = fpold-fp c store the sum of squared residuals for the current set of knots. 170 fpold = fp c find nplx, the number of knots we should add in the x-direction. nplx = 1 if(nx.eq.nminx) go to 180 npl1 = nplusx*2 rn = nplusx if(reducx.gt.acc) npl1 = rn*fpms/reducx nplx = min0(nplusx*2,max0(npl1,nplusx/2,1)) c find nply, the number of knots we should add in the y-direction. 180 nply = 1 if(ny.eq.nminy) go to 190 npl1 = nplusy*2 rn = nplusy if(reducy.gt.acc) npl1 = rn*fpms/reducy nply = min0(nplusy*2,max0(npl1,nplusy/2,1)) 190 if(nplx-nply) 210,200,230 200 if(lastdi.lt.0) go to 230 210 if(nx.eq.nxe) go to 230 c addition in the x-direction. lastdi = -1 nplusx = nplx ifsx = 0 do 220 l=1,nplusx c add a new knot in the x-direction call fpknot(x,mx,tx,nx,fpintx,nrdatx,nrintx,nxest,1) c test whether we cannot further increase the number of knots in the c x-direction. if(nx.eq.nxe) go to 250 220 continue go to 250 230 if(ny.eq.nye) go to 210 c addition in the y-direction. lastdi = 1 nplusy = nply ifsy = 0 do 240 l=1,nplusy c add a new knot in the y-direction. call fpknot(y,my,ty,ny,fpinty,nrdaty,nrinty,nyest,1) c test whether we cannot further increase the number of knots in the c y-direction. if(ny.eq.nye) go to 250 240 continue c restart the computations with the new set of knots. 250 continue c test whether the least-squares polynomial is a solution of our c approximation problem. 300 if(ier.eq.(-2)) go to 440 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c part 2: determination of the smoothing spline sp(x,y) c c ***************************************************** c c we have determined the number of knots and their position. we now c c compute the b-spline coefficients of the smoothing spline sp(x,y). c c this smoothing spline varies with the parameter p in such a way thatc c f(p) = sumi=1,mx(sumj=1,my((z(i,j)-sp(x(i),y(j)))**2) c c is a continuous, strictly decreasing function of p. moreover the c c least-squares polynomial corresponds to p=0 and the least-squares c c spline to p=infinity. iteratively we then have to determine the c c positive value of p such that f(p)=s. the process which is proposed c c here makes use of rational interpolation. f(p) is approximated by a c c rational function r(p)=(u*p+v)/(p+w); three values of p (p1,p2,p3) c c with corresponding values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s)c c are used to calculate the new value of p such that r(p)=s. c c convergence is guaranteed by taking f1 > 0 and f3 < 0. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c initial value for p. p1 = 0. f1 = fp0-s p3 = -one f3 = fpms p = one ich1 = 0 ich3 = 0 c iteration process to find the root of f(p)=s. do 350 iter = 1,maxit c find the smoothing spline sp(x,y) and the corresponding sum of c squared residuals fp. call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, * ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), * wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), * nrx,nry) c test whether the approximation sp(x,y) is an acceptable solution. fpms = fp-s if(abs(fpms).lt.acc) go to 440 c test whether the maximum allowable number of iterations has been c reached. if(iter.eq.maxit) go to 400 c carry out one more step of the iteration process. p2 = p f2 = fpms if(ich3.ne.0) go to 320 if((f2-f3).gt.acc) go to 310 c our initial choice of p is too large. p3 = p2 f3 = f2 p = p*con4 if(p.le.p1) p = p1*con9 + p2*con1 go to 350 310 if(f2.lt.0.) ich3 = 1 320 if(ich1.ne.0) go to 340 if((f1-f2).gt.acc) go to 330 c our initial choice of p is too small p1 = p2 f1 = f2 p = p/con4 if(p3.lt.0.) go to 350 if(p.ge.p3) p = p2*con1 + p3*con9 go to 350 c test whether the iteration process proceeds as theoretically c expected. 330 if(f2.gt.0.) ich1 = 1 340 if(f2.ge.f1 .or. f2.le.f3) go to 410 c find the new value of p. p = fprati(p1,f1,p2,f2,p3,f3) 350 continue c error codes and messages. 400 ier = 3 go to 440 410 ier = 2 go to 440 420 ier = 1 go to 440 430 ier = -1 fp = 0. 440 return end