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