subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest,
* nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest,
* nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h,
* index,nummer,wrk,lwrk,ier)
c ..
c ..scalar arguments..
real xb,xe,yb,ye,s,eta,tol,fp,fp0
integer iopt,m,kxx,kyy,nxest,nyest,maxit,nmax,km1,km2,ib1,ib3,
* nc,intest,nrest,nx0,ny0,lwrk,ier
c ..array arguments..
real x(m),y(m),z(m),w(m),tx(nmax),ty(nmax),c(nc),fpint(intest),
* coord(intest),f(nc),ff(nc),a(nc,ib1),q(nc,ib3),bx(nmax,km2),
* by(nmax,km2),spx(m,km1),spy(m,km1),h(ib3),wrk(lwrk)
integer index(nrest),nummer(m)
c ..local scalars..
real acc,arg,cos,dmax,fac1,fac2,fpmax,fpms,f1,f2,f3,hxi,p,pinv,
* piv,p1,p2,p3,sigma,sin,sq,store,wi,x0,x1,y0,y1,zi,eps,
* rn,one,con1,con9,con4,half,ten
integer i,iband,iband1,iband3,iband4,ibb,ichang,ich1,ich3,ii,
* in,irot,iter,i1,i2,i3,j,jrot,jxy,j1,kx,kx1,kx2,ky,ky1,ky2,l,
* la,lf,lh,lwest,lx,ly,l1,l2,n,ncof,nk1x,nk1y,nminx,nminy,nreg,
* nrint,num,num1,nx,nxe,nxx,ny,nye,nyy,n1,rank
c ..local arrays..
real hx(6),hy(6)
c ..function references..
real abs,fprati,sqrt
integer min0
c ..subroutine references..
c fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota
c ..
c set constants
one = 0.1e+01
con1 = 0.1e0
con9 = 0.9e0
con4 = 0.4e-01
half = 0.5e0
ten = 0.1e+02
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 weighted 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 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=nminy=2*ky+2. c
c fp0=f(0) denotes the corresponding weighted sum of squared c
c residuals c
c if 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
c eventually the independent variables x and y (and the corresponding c
c parameters) will be switched if this can reduce the bandwidth of the c
c system to be solved. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ichang denotes whether(1) or not(-1) the directions have been inter-
c changed.
ichang = -1
x0 = xb
x1 = xe
y0 = yb
y1 = ye
kx = kxx
ky = kyy
kx1 = kx+1
ky1 = ky+1
nxe = nxest
nye = nyest
eps = sqrt(eta)
if(iopt.lt.0) go to 20
c calculation of acc, the absolute tolerance for the root of f(p)=s.
acc = tol*s
if(iopt.eq.0) go to 10
if(fp0.gt.s) go to 20
c initialization for the least-squares polynomial.
10 nminx = 2*kx1
nminy = 2*ky1
nx = nminx
ny = nminy
ier = -2
go to 30
20 nx = nx0
ny = ny0
c main loop for the different sets of knots. m is a save upper bound
c for the number of trials.
30 do 420 iter=1,m
c find the position of the additional knots which are needed for the
c b-spline representation of s(x,y).
l = nx
do 40 i=1,kx1
tx(i) = x0
tx(l) = x1
l = l-1
40 continue
l = ny
do 50 i=1,ky1
ty(i) = y0
ty(l) = y1
l = l-1
50 continue
c find nrint, the total number of knot intervals and nreg, the number
c of panels in which the approximation domain is subdivided by the
c intersection of knots.
nxx = nx-2*kx1+1
nyy = ny-2*ky1+1
nrint = nxx+nyy
nreg = nxx*nyy
c find the bandwidth of the observation matrix a.
c if necessary, interchange the variables x and y, in order to obtain
c a minimal bandwidth.
iband1 = kx*(ny-ky1)+ky
l = ky*(nx-kx1)+kx
if(iband1.le.l) go to 130
iband1 = l
ichang = -ichang
do 60 i=1,m
store = x(i)
x(i) = y(i)
y(i) = store
60 continue
store = x0
x0 = y0
y0 = store
store = x1
x1 = y1
y1 = store
n = min0(nx,ny)
do 70 i=1,n
store = tx(i)
tx(i) = ty(i)
ty(i) = store
70 continue
n1 = n+1
if(nx-ny) 80,120,100
80 do 90 i=n1,ny
tx(i) = ty(i)
90 continue
go to 120
100 do 110 i=n1,nx
ty(i) = tx(i)
110 continue
120 l = nx
nx = ny
ny = l
l = nxe
nxe = nye
nye = l
l = nxx
nxx = nyy
nyy = l
l = kx
kx = ky
ky = l
kx1 = kx+1
ky1 = ky+1
130 iband = iband1+1
c arrange the data points according to the panel they belong to.
call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg)
c find ncof, the number of b-spline coefficients.
nk1x = nx-kx1
nk1y = ny-ky1
ncof = nk1x*nk1y
c initialize the observation matrix a.
do 140 i=1,ncof
f(i) = 0.
do 140 j=1,iband
a(i,j) = 0.
140 continue
c initialize the sum of squared residuals.
fp = 0.
c fetch the data points in the new order. main loop for the
c different panels.
do 250 num=1,nreg
c fix certain constants for the current panel; jrot records the column
c number of the first non-zero element in a row of the observation
c matrix according to a data point of the panel.
num1 = num-1
lx = num1/nyy
l1 = lx+kx1
ly = num1-lx*nyy
l2 = ly+ky1
jrot = lx*nk1y+ly
c test whether there are still data points in the panel.
in = index(num)
150 if(in.eq.0) go to 250
c fetch a new data point.
wi = w(in)
zi = z(in)*wi
c evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in).
call fpbspl(tx,nx,kx,x(in),l1,hx)
c evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in).
call fpbspl(ty,ny,ky,y(in),l2,hy)
c store the value of these b-splines in spx and spy respectively.
do 160 i=1,kx1
spx(in,i) = hx(i)
160 continue
do 170 i=1,ky1
spy(in,i) = hy(i)
170 continue
c initialize the new row of observation matrix.
do 180 i=1,iband
h(i) = 0.
180 continue
c calculate the non-zero elements of the new row by making the cross
c products of the non-zero b-splines in x- and y-direction.
i1 = 0
do 200 i=1,kx1
hxi = hx(i)
j1 = i1
do 190 j=1,ky1
j1 = j1+1
h(j1) = hxi*hy(j)*wi
190 continue
i1 = i1+nk1y
200 continue
c rotate the row into triangle by givens transformations .
irot = jrot
do 220 i=1,iband
irot = irot+1
piv = h(i)
if(piv.eq.0.) go to 220
c calculate the parameters of the givens transformation.
call fpgivs(piv,a(irot,1),cos,sin)
c apply that transformation to the right hand side.
call fprota(cos,sin,zi,f(irot))
if(i.eq.iband) go to 230
c apply that transformation to the left hand side.
i2 = 1
i3 = i+1
do 210 j=i3,iband
i2 = i2+1
call fprota(cos,sin,h(j),a(irot,i2))
210 continue
220 continue
c add the contribution of the row to the sum of squares of residual
c right hand sides.
230 fp = fp+zi**2
c find the number of the next data point in the panel.
240 in = nummer(in)
go to 150
250 continue
c find dmax, the maximum value for the diagonal elements in the reduced
c triangle.
dmax = 0.
do 260 i=1,ncof
if(a(i,1).le.dmax) go to 260
dmax = a(i,1)
260 continue
c check whether the observation matrix is rank deficient.
sigma = eps*dmax
do 270 i=1,ncof
if(a(i,1).le.sigma) go to 280
270 continue
c backward substitution in case of full rank.
call fpback(a,f,ncof,iband,c,nc)
rank = ncof
do 275 i=1,ncof
q(i,1) = a(i,1)/dmax
275 continue
go to 300
c in case of rank deficiency, find the minimum norm solution.
c check whether there is sufficient working space
280 lwest = ncof*iband+ncof+iband
if(lwrk.lt.lwest) go to 780
do 290 i=1,ncof
ff(i) = f(i)
do 290 j=1,iband
q(i,j) = a(i,j)
290 continue
lf =1
lh = lf+ncof
la = lh+iband
call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la),
* wrk(lf),wrk(lh))
do 295 i=1,ncof
q(i,1) = q(i,1)/dmax
295 continue
c add to the sum of squared residuals, the contribution of reducing
c the rank.
fp = fp+sq
300 if(ier.eq.(-2)) fp0 = fp
c test whether the least-squares spline is an acceptable solution.
if(iopt.lt.0) go to 820
fpms = fp-s
if(abs(fpms).le.acc) if(fp) 815,815,820
c test whether we can accept the choice of knots.
if(fpms.lt.0.) go to 430
c test whether we cannot further increase the number of knots.
if(ncof.gt.m) go to 790
ier = 0
c search where to add a new knot.
c find for each interval the sum of squared residuals fpint for the
c data points having the coordinate belonging to that knot interval.
c calculate also coord which is the same sum, weighted by the position
c of the data points considered.
310 do 320 i=1,nrint
fpint(i) = 0.
coord(i) = 0.
320 continue
do 360 num=1,nreg
num1 = num-1
lx = num1/nyy
l1 = lx+1
ly = num1-lx*nyy
l2 = ly+1+nxx
jrot = lx*nk1y+ly
in = index(num)
330 if(in.eq.0) go to 360
store = 0.
i1 = jrot
do 350 i=1,kx1
hxi = spx(in,i)
j1 = i1
do 340 j=1,ky1
j1 = j1+1
store = store+hxi*spy(in,j)*c(j1)
340 continue
i1 = i1+nk1y
350 continue
store = (w(in)*(z(in)-store))**2
fpint(l1) = fpint(l1)+store
coord(l1) = coord(l1)+store*x(in)
fpint(l2) = fpint(l2)+store
coord(l2) = coord(l2)+store*y(in)
in = nummer(in)
go to 330
360 continue
c find the interval for which fpint is maximal on the condition that
c there still can be added a knot.
370 l = 0
fpmax = 0.
l1 = 1
l2 = nrint
if(nx.eq.nxe) l1 = nxx+1
if(ny.eq.nye) l2 = nxx
if(l1.gt.l2) go to 810
do 380 i=l1,l2
if(fpmax.ge.fpint(i)) go to 380
l = i
fpmax = fpint(i)
380 continue
c test whether we cannot further increase the number of knots.
if(l.eq.0) go to 785
c calculate the position of the new knot.
arg = coord(l)/fpint(l)
c test in what direction the new knot is going to be added.
if(l.gt.nxx) go to 400
c addition in the x-direction.
jxy = l+kx1
fpint(l) = 0.
fac1 = tx(jxy)-arg
fac2 = arg-tx(jxy-1)
if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370
j = nx
do 390 i=jxy,nx
tx(j+1) = tx(j)
j = j-1
390 continue
tx(jxy) = arg
nx = nx+1
go to 420
c addition in the y-direction.
400 jxy = l+ky1-nxx
fpint(l) = 0.
fac1 = ty(jxy)-arg
fac2 = arg-ty(jxy-1)
if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370
j = ny
do 410 i=jxy,ny
ty(j+1) = ty(j)
j = j-1
410 continue
ty(jxy) = arg
ny = ny+1
c restart the computations with the new set of knots.
420 continue
c test whether the least-squares polynomial is a solution of our
c approximation problem.
430 if(ier.eq.(-2)) go to 830
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 the observation matrix a is extended by the rows of a matrix, c
c expressing that sp(x,y) must be a polynomial of degree kx in x and c
c ky in y. the corresponding weights of these additional rows are set c
c to 1./p. iteratively we than have to determine the value of p c
c such that f(p)=sum((w(i)*(z(i)-sp(x(i),y(i))))**2) be = s. c
c we already know that the least-squares polynomial corresponds to c
c p=0 and that the least-squares spline corresponds to p=infinity. c
c the iteration process which is proposed here makes use of rational c
c interpolation. since f(p) is a convex and strictly decreasing c
c function of p, it can be approximated by a rational function r(p)= c
c (u*p+v)/(p+w). three values of p(p1,p2,p3) with corresponding values c
c of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the c
c new value of p such that r(p)=s. convergence is guaranteed by taking c
c f1 > 0 and f3 < 0. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
kx2 = kx1+1
c test whether there are interior knots in the x-direction.
if(nk1x.eq.kx1) go to 440
c evaluate the discotinuity jumps of the kx-th order derivative of
c the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1.
call fpdisc(tx,nx,kx2,bx,nmax)
440 ky2 = ky1 + 1
c test whether there are interior knots in the y-direction.
if(nk1y.eq.ky1) go to 450
c evaluate the discontinuity jumps of the ky-th order derivative of
c the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1.
call fpdisc(ty,ny,ky2,by,nmax)
c initial value for p.
450 p1 = 0.
f1 = fp0-s
p3 = -one
f3 = fpms
p = 0.
do 460 i=1,ncof
p = p+a(i,1)
460 continue
rn = ncof
p = rn/p
c find the bandwidth of the extended observation matrix.
iband3 = kx1*nk1y
iband4 = iband3 +1
ich1 = 0
ich3 = 0
c iteration process to find the root of f(p)=s.
do 770 iter=1,maxit
pinv = one/p
c store the triangularized observation matrix into q.
do 480 i=1,ncof
ff(i) = f(i)
do 470 j=1,iband
q(i,j) = a(i,j)
470 continue
ibb = iband+1
do 480 j=ibb,iband4
q(i,j) = 0.
480 continue
if(nk1y.eq.ky1) go to 560
c extend the observation matrix with the rows of a matrix, expressing
c that for x=cst. sp(x,y) must be a polynomial in y of degree ky.
do 550 i=ky2,nk1y
ii = i-ky1
do 550 j=1,nk1x
c initialize the new row.
do 490 l=1,iband
h(l) = 0.
490 continue
c fill in the non-zero elements of the row. jrot records the column
c number of the first non-zero element in the row.
do 500 l=1,ky2
h(l) = by(ii,l)*pinv
500 continue
zi = 0.
jrot = (j-1)*nk1y+ii
c rotate the new row into triangle by givens transformations without
c square roots.
do 540 irot=jrot,ncof
piv = h(1)
i2 = min0(iband1,ncof-irot)
if(piv.eq.0.) if(i2) 550,550,520
c calculate the parameters of the givens transformation.
call fpgivs(piv,q(irot,1),cos,sin)
c apply that givens transformation to the right hand side.
call fprota(cos,sin,zi,ff(irot))
if(i2.eq.0) go to 550
c apply that givens transformation to the left hand side.
do 510 l=1,i2
l1 = l+1
call fprota(cos,sin,h(l1),q(irot,l1))
510 continue
520 do 530 l=1,i2
h(l) = h(l+1)
530 continue
h(i2+1) = 0.
540 continue
550 continue
560 if(nk1x.eq.kx1) go to 640
c extend the observation matrix with the rows of a matrix expressing
c that for y=cst. sp(x,y) must be a polynomial in x of degree kx.
do 630 i=kx2,nk1x
ii = i-kx1
do 630 j=1,nk1y
c initialize the new row
do 570 l=1,iband4
h(l) = 0.
570 continue
c fill in the non-zero elements of the row. jrot records the column
c number of the first non-zero element in the row.
j1 = 1
do 580 l=1,kx2
h(j1) = bx(ii,l)*pinv
j1 = j1+nk1y
580 continue
zi = 0.
jrot = (i-kx2)*nk1y+j
c rotate the new row into triangle by givens transformations .
do 620 irot=jrot,ncof
piv = h(1)
i2 = min0(iband3,ncof-irot)
if(piv.eq.0.) if(i2) 630,630,600
c calculate the parameters of the givens transformation.
call fpgivs(piv,q(irot,1),cos,sin)
c apply that givens transformation to the right hand side.
call fprota(cos,sin,zi,ff(irot))
if(i2.eq.0) go to 630
c apply that givens transformation to the left hand side.
do 590 l=1,i2
l1 = l+1
call fprota(cos,sin,h(l1),q(irot,l1))
590 continue
600 do 610 l=1,i2
h(l) = h(l+1)
610 continue
h(i2+1) = 0.
620 continue
630 continue
c find dmax, the maximum value for the diagonal elements in the
c reduced triangle.
640 dmax = 0.
do 650 i=1,ncof
if(q(i,1).le.dmax) go to 650
dmax = q(i,1)
650 continue
c check whether the matrix is rank deficient.
sigma = eps*dmax
do 660 i=1,ncof
if(q(i,1).le.sigma) go to 670
660 continue
c backward substitution in case of full rank.
call fpback(q,ff,ncof,iband4,c,nc)
rank = ncof
go to 675
c in case of rank deficiency, find the minimum norm solution.
670 lwest = ncof*iband4+ncof+iband4
if(lwrk.lt.lwest) go to 780
lf = 1
lh = lf+ncof
la = lh+iband4
call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la),
* wrk(lf),wrk(lh))
675 do 680 i=1,ncof
q(i,1) = q(i,1)/dmax
680 continue
c compute f(p).
fp = 0.
do 720 num = 1,nreg
num1 = num-1
lx = num1/nyy
ly = num1-lx*nyy
jrot = lx*nk1y+ly
in = index(num)
690 if(in.eq.0) go to 720
store = 0.
i1 = jrot
do 710 i=1,kx1
hxi = spx(in,i)
j1 = i1
do 700 j=1,ky1
j1 = j1+1
store = store+hxi*spy(in,j)*c(j1)
700 continue
i1 = i1+nk1y
710 continue
fp = fp+(w(in)*(z(in)-store))**2
in = nummer(in)
go to 690
720 continue
c test whether the approximation sp(x,y) is an acceptable solution.
fpms = fp-s
if(abs(fpms).le.acc) go to 820
c test whether the maximum allowable number of iterations has been
c reached.
if(iter.eq.maxit) go to 795
c carry out one more step of the iteration process.
p2 = p
f2 = fpms
if(ich3.ne.0) go to 740
if((f2-f3).gt.acc) go to 730
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 770
730 if(f2.lt.0.) ich3 = 1
740 if(ich1.ne.0) go to 760
if((f1-f2).gt.acc) go to 750
c our initial choice of p is too small
p1 = p2
f1 = f2
p = p/con4
if(p3.lt.0.) go to 770
if(p.ge.p3) p = p2*con1 + p3*con9
go to 770
750 if(f2.gt.0.) ich1 = 1
c test whether the iteration process proceeds as theoretically
c expected.
760 if(f2.ge.f1 .or. f2.le.f3) go to 800
c find the new value of p.
p = fprati(p1,f1,p2,f2,p3,f3)
770 continue
c error codes and messages.
780 ier = lwest
go to 830
785 ier = 5
go to 830
790 ier = 4
go to 830
795 ier = 3
go to 830
800 ier = 2
go to 830
810 ier = 1
go to 830
815 ier = -1
fp = 0.
820 if(ncof.ne.rank) ier = -rank
c test whether x and y are in the original order.
830 if(ichang.lt.0) go to 930
c if not, interchange x and y once more.
l1 = 1
do 840 i=1,nk1x
l2 = i
do 840 j=1,nk1y
f(l2) = c(l1)
l1 = l1+1
l2 = l2+nk1x
840 continue
do 850 i=1,ncof
c(i) = f(i)
850 continue
do 860 i=1,m
store = x(i)
x(i) = y(i)
y(i) = store
860 continue
n = min0(nx,ny)
do 870 i=1,n
store = tx(i)
tx(i) = ty(i)
ty(i) = store
870 continue
n1 = n+1
if(nx-ny) 880,920,900
880 do 890 i=n1,ny
tx(i) = ty(i)
890 continue
go to 920
900 do 910 i=n1,nx
ty(i) = tx(i)
910 continue
920 l = nx
nx = ny
ny = l
930 if(iopt.lt.0) go to 940
nx0 = nx
ny0 = ny
940 return
end