Timur A. Fatkhullin
5279d1c41a
start developing of FITPACK C++ bindings mount_server.cpp: fix compilation error with GCC15
361 lines
13 KiB
Fortran
361 lines
13 KiB
Fortran
recursive subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,
|
|
* maxit,k1,k2,n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
|
|
implicit none
|
|
c ..
|
|
c ..scalar arguments..
|
|
real*8 xb,xe,s,tol,fp
|
|
integer iopt,m,k,nest,maxit,k1,k2,n,ier
|
|
c ..array arguments..
|
|
real*8 x(m),y(m),w(m),t(nest),c(nest),fpint(nest),
|
|
* z(nest),a(nest,k1),b(nest,k2),g(nest,k2),q(m,k1)
|
|
integer nrdata(nest)
|
|
c ..local scalars..
|
|
real*8 acc,con1,con4,con9,cos,half,fpart,fpms,fpold,fp0,f1,f2,f3,
|
|
* one,p,pinv,piv,p1,p2,p3,rn,sin,store,term,wi,xi,yi
|
|
integer i,ich1,ich3,it,iter,i1,i2,i3,j,k3,l,l0,
|
|
* mk1,new,nk1,nmax,nmin,nplus,npl1,nrint,n8
|
|
c ..local arrays..
|
|
real*8 h(7)
|
|
c ..function references
|
|
real*8 abs,fprati
|
|
integer max0,min0
|
|
c ..subroutine references..
|
|
c fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
|
|
c ..
|
|
c set constants
|
|
one = 0.1d+01
|
|
con1 = 0.1d0
|
|
con9 = 0.9d0
|
|
con4 = 0.4d-01
|
|
half = 0.5d0
|
|
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), c
|
|
c and the corresponding sum of squared residuals fp=f(p=inf). c
|
|
c if iopt=-1 sinf(x) 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 nmax = m+k+1. c
|
|
c if s > 0 and c
|
|
c iopt=0 we first compute the least-squares polynomial of c
|
|
c degree k; n = nmin = 2*k+2 c
|
|
c iopt=1 we start with the set of knots found at the last c
|
|
c call of the routine, except for the case that s > fp0; then c
|
|
c we compute directly the least-squares polynomial of degree k. c
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
c determine nmin, the number of knots for polynomial approximation.
|
|
nmin = 2*k1
|
|
if(iopt.lt.0) go to 60
|
|
c calculation of acc, the absolute tolerance for the root of f(p)=s.
|
|
acc = tol*s
|
|
c determine nmax, the number of knots for spline interpolation.
|
|
nmax = m+k1
|
|
if(s.gt.0.0d0) go to 45
|
|
c if s=0, s(x) is an interpolating spline.
|
|
c test whether the required storage space exceeds the available one.
|
|
n = nmax
|
|
if(nmax.gt.nest) go to 420
|
|
c find the position of the interior knots in case of interpolation.
|
|
10 mk1 = m-k1
|
|
if(mk1.eq.0) go to 60
|
|
k3 = k/2
|
|
i = k2
|
|
j = k3+2
|
|
if(k3*2.eq.k) go to 30
|
|
do 20 l=1,mk1
|
|
t(i) = x(j)
|
|
i = i+1
|
|
j = j+1
|
|
20 continue
|
|
go to 60
|
|
30 do 40 l=1,mk1
|
|
t(i) = (x(j)+x(j-1))*half
|
|
i = i+1
|
|
j = j+1
|
|
40 continue
|
|
go to 60
|
|
c if s>0 our initial choice of knots depends on the value of iopt.
|
|
c if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
|
|
c polynomial of degree k which is a spline without interior knots.
|
|
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.
|
|
45 if(iopt.eq.0) go to 50
|
|
if(n.eq.nmin) go to 50
|
|
fp0 = fpint(n)
|
|
fpold = fpint(n-1)
|
|
nplus = nrdata(n)
|
|
if(fp0.gt.s) go to 60
|
|
50 n = nmin
|
|
fpold = 0.0d0
|
|
nplus = 0
|
|
nrdata(1) = m-2
|
|
c main loop for the different sets of knots. m is a save upper bound
|
|
c for the number of trials.
|
|
60 do 200 iter = 1,m
|
|
if(n.eq.nmin) ier = -2
|
|
c find nrint, tne number of knot intervals.
|
|
nrint = n-nmin+1
|
|
c find the position of the additional knots which are needed for
|
|
c the b-spline representation of s(x).
|
|
nk1 = n-k1
|
|
i = n
|
|
do 70 j=1,k1
|
|
t(j) = xb
|
|
t(i) = xe
|
|
i = i-1
|
|
70 continue
|
|
c compute the b-spline coefficients of the least-squares spline
|
|
c sinf(x). the observation matrix a is built up row by row and
|
|
c reduced to upper triangular form by givens transformations.
|
|
c at the same time fp=f(p=inf) is computed.
|
|
fp = 0.0d0
|
|
c initialize the observation matrix a.
|
|
do 80 i=1,nk1
|
|
z(i) = 0.0d0
|
|
do 80 j=1,k1
|
|
a(i,j) = 0.0d0
|
|
80 continue
|
|
l = k1
|
|
do 130 it=1,m
|
|
c fetch the current data point x(it),y(it).
|
|
xi = x(it)
|
|
wi = w(it)
|
|
yi = y(it)*wi
|
|
c search for knot interval t(l) <= xi < t(l+1).
|
|
85 if(xi.lt.t(l+1) .or. l.eq.nk1) go to 90
|
|
l = l+1
|
|
go to 85
|
|
c evaluate the (k+1) non-zero b-splines at xi and store them in q.
|
|
90 call fpbspl(t,n,k,xi,l,h)
|
|
do 95 i=1,k1
|
|
q(it,i) = h(i)
|
|
h(i) = h(i)*wi
|
|
95 continue
|
|
c rotate the new row of the observation matrix into triangle.
|
|
j = l-k1
|
|
do 110 i=1,k1
|
|
j = j+1
|
|
piv = h(i)
|
|
if(piv.eq.0.0d0) go to 110
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,a(j,1),cos,sin)
|
|
c transformations to right hand side.
|
|
call fprota(cos,sin,yi,z(j))
|
|
if(i.eq.k1) go to 120
|
|
i2 = 1
|
|
i3 = i+1
|
|
do 100 i1 = i3,k1
|
|
i2 = i2+1
|
|
c transformations to left hand side.
|
|
call fprota(cos,sin,h(i1),a(j,i2))
|
|
100 continue
|
|
110 continue
|
|
c add contribution of this row to the sum of squares of residual
|
|
c right hand sides.
|
|
120 fp = fp+yi*yi
|
|
130 continue
|
|
if(ier.eq.(-2)) fp0 = fp
|
|
fpint(n) = fp0
|
|
fpint(n-1) = fpold
|
|
nrdata(n) = nplus
|
|
c backward substitution to obtain the b-spline coefficients.
|
|
call fpback(a,z,nk1,k1,c,nest)
|
|
c test whether the approximation sinf(x) 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 accept the choice of knots.
|
|
if(fpms.lt.0.0d0) go to 250
|
|
c if n = nmax, sinf(x) is an interpolating spline.
|
|
if(n.eq.nmax) go to 430
|
|
c increase the number of knots.
|
|
c if n=nest we cannot increase the number of knots because of
|
|
c the storage capacity limitation.
|
|
if(n.eq.nest) go to 420
|
|
c determine the number of knots nplus we are going to add.
|
|
if(ier.eq.0) go to 140
|
|
nplus = 1
|
|
ier = 0
|
|
go to 150
|
|
140 npl1 = nplus*2
|
|
rn = nplus
|
|
if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
|
|
nplus = min0(nplus*2,max0(npl1,nplus/2,1))
|
|
150 fpold = fp
|
|
c compute the sum((w(i)*(y(i)-s(x(i))))**2) for each knot interval
|
|
c t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
|
|
fpart = 0.0d0
|
|
i = 1
|
|
l = k2
|
|
new = 0
|
|
do 180 it=1,m
|
|
if(x(it).lt.t(l) .or. l.gt.nk1) go to 160
|
|
new = 1
|
|
l = l+1
|
|
160 term = 0.0d0
|
|
l0 = l-k2
|
|
do 170 j=1,k1
|
|
l0 = l0+1
|
|
term = term+c(l0)*q(it,j)
|
|
170 continue
|
|
term = (w(it)*(term-y(it)))**2
|
|
fpart = fpart+term
|
|
if(new.eq.0) go to 180
|
|
store = term*half
|
|
fpint(i) = fpart-store
|
|
i = i+1
|
|
fpart = store
|
|
new = 0
|
|
180 continue
|
|
fpint(nrint) = fpart
|
|
do 190 l=1,nplus
|
|
c add a new knot.
|
|
call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
|
|
c if n=nmax we locate the knots as for interpolation.
|
|
if(n.eq.nmax) go to 10
|
|
c test whether we cannot further increase the number of knots.
|
|
if(n.eq.nest) go to 200
|
|
190 continue
|
|
c restart the computations with the new set of knots.
|
|
200 continue
|
|
c test whether the least-squares kth degree polynomial is a solution
|
|
c of our approximation problem.
|
|
250 if(ier.eq.(-2)) go to 440
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
c part 2: determination of the smoothing spline sp(x). c
|
|
c *************************************************** c
|
|
c we have determined the number of knots and their position. c
|
|
c we now compute the b-spline coefficients of the smoothing spline c
|
|
c sp(x). the observation matrix a is extended by the rows of matrix c
|
|
c b expressing that the kth derivative discontinuities of sp(x) at c
|
|
c the interior knots t(k+2),...t(n-k-1) must be zero. the corres- c
|
|
c ponding weights of these additional rows are set to 1/p. c
|
|
c iteratively we then have to determine the value of p such that c
|
|
c f(p)=sum((w(i)*(y(i)-sp(x(i))))**2) be = s. we already know that c
|
|
c the least-squares kth degree polynomial corresponds to p=0, and c
|
|
c that the least-squares spline corresponds to p=infinity. the c
|
|
c 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 c
|
|
c r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond- c
|
|
c ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used c
|
|
c to calculate the new value of p such that r(p)=s. convergence is c
|
|
c guaranteed by taking f1>0 and f3<0. c
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
c evaluate the discontinuity jump of the kth derivative of the
|
|
c b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
|
|
call fpdisc(t,n,k2,b,nest)
|
|
c initial value for p.
|
|
p1 = 0.0d0
|
|
f1 = fp0-s
|
|
p3 = -one
|
|
f3 = fpms
|
|
p = 0.
|
|
do 255 i=1,nk1
|
|
p = p+a(i,1)
|
|
255 continue
|
|
rn = nk1
|
|
p = rn/p
|
|
ich1 = 0
|
|
ich3 = 0
|
|
n8 = n-nmin
|
|
c iteration process to find the root of f(p) = s.
|
|
do 360 iter=1,maxit
|
|
c the rows of matrix b with weight 1/p are rotated into the
|
|
c triangularised observation matrix a which is stored in g.
|
|
pinv = one/p
|
|
do 260 i=1,nk1
|
|
c(i) = z(i)
|
|
g(i,k2) = 0.0d0
|
|
do 260 j=1,k1
|
|
g(i,j) = a(i,j)
|
|
260 continue
|
|
do 300 it=1,n8
|
|
c the row of matrix b is rotated into triangle by givens transformation
|
|
do 270 i=1,k2
|
|
h(i) = b(it,i)*pinv
|
|
270 continue
|
|
yi = 0.0d0
|
|
do 290 j=it,nk1
|
|
piv = h(1)
|
|
c calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,g(j,1),cos,sin)
|
|
c transformations to right hand side.
|
|
call fprota(cos,sin,yi,c(j))
|
|
if(j.eq.nk1) go to 300
|
|
i2 = k1
|
|
if(j.gt.n8) i2 = nk1-j
|
|
do 280 i=1,i2
|
|
c transformations to left hand side.
|
|
i1 = i+1
|
|
call fprota(cos,sin,h(i1),g(j,i1))
|
|
h(i) = h(i1)
|
|
280 continue
|
|
h(i2+1) = 0.0d0
|
|
290 continue
|
|
300 continue
|
|
c backward substitution to obtain the b-spline coefficients.
|
|
call fpback(g,c,nk1,k2,c,nest)
|
|
c computation of f(p).
|
|
fp = 0.0d0
|
|
l = k2
|
|
do 330 it=1,m
|
|
if(x(it).lt.t(l) .or. l.gt.nk1) go to 310
|
|
l = l+1
|
|
310 l0 = l-k2
|
|
term = 0.0d0
|
|
do 320 j=1,k1
|
|
l0 = l0+1
|
|
term = term+c(l0)*q(it,j)
|
|
320 continue
|
|
fp = fp+(w(it)*(term-y(it)))**2
|
|
330 continue
|
|
c test whether the approximation sp(x) is an acceptable solution.
|
|
fpms = fp-s
|
|
if(abs(fpms).lt.acc) go to 440
|
|
c test whether the maximal number of iterations is 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 340
|
|
if((f2-f3).gt.acc) go to 335
|
|
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 360
|
|
335 if(f2.lt.0.0d0) ich3=1
|
|
340 if(ich1.ne.0) go to 350
|
|
if((f1-f2).gt.acc) go to 345
|
|
c our initial choice of p is too small
|
|
p1 = p2
|
|
f1 = f2
|
|
p = p/con4
|
|
if(p3.lt.0.) go to 360
|
|
if(p.ge.p3) p = p2*con1 + p3*con9
|
|
go to 360
|
|
345 if(f2.gt.0.0d0) ich1=1
|
|
c test whether the iteration process proceeds as theoretically
|
|
c expected.
|
|
350 if(f2.ge.f1 .or. f2.le.f3) go to 410
|
|
c find the new value for p.
|
|
p = fprati(p1,f1,p2,f2,p3,f3)
|
|
360 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
|
|
440 return
|
|
end
|