74 lines
2.2 KiB
Fortran
74 lines
2.2 KiB
Fortran
recursive subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,
|
|
* istart)
|
|
implicit none
|
|
c subroutine fpknot locates an additional knot for a spline of degree
|
|
c k and adjusts the corresponding parameters,i.e.
|
|
c t : the position of the knots.
|
|
c n : the number of knots.
|
|
c nrint : the number of knotintervals.
|
|
c fpint : the sum of squares of residual right hand sides
|
|
c for each knot interval.
|
|
c nrdata: the number of data points inside each knot interval.
|
|
c istart indicates that the smallest data point at which the new knot
|
|
c may be added is x(istart+1)
|
|
c ..
|
|
c ..scalar arguments..
|
|
integer m,n,nrint,nest,istart
|
|
c ..array arguments..
|
|
real*8 x(m),t(nest),fpint(nest)
|
|
integer nrdata(nest)
|
|
c ..local scalars..
|
|
real*8 an,am,fpmax
|
|
integer ihalf,j,jbegin,jj,jk,jpoint,k,maxbeg,maxpt,
|
|
* next,nrx,number
|
|
|
|
c note: do not initialize on the same line to avoid saving between calls
|
|
logical iserr
|
|
iserr = .TRUE.
|
|
|
|
k = (n-nrint-1)/2
|
|
c search for knot interval t(number+k) <= x <= t(number+k+1) where
|
|
c fpint(number) is maximal on the condition that nrdata(number)
|
|
c not equals zero.
|
|
fpmax = 0.
|
|
jbegin = istart
|
|
do 20 j=1,nrint
|
|
jpoint = nrdata(j)
|
|
if(fpmax.ge.fpint(j) .or. jpoint.eq.0) go to 10
|
|
iserr = .FALSE.
|
|
fpmax = fpint(j)
|
|
number = j
|
|
maxpt = jpoint
|
|
maxbeg = jbegin
|
|
10 jbegin = jbegin+jpoint+1
|
|
20 continue
|
|
c error condition detected, go to exit
|
|
if(iserr) go to 50
|
|
c let coincide the new knot t(number+k+1) with a data point x(nrx)
|
|
c inside the old knot interval t(number+k) <= x <= t(number+k+1).
|
|
ihalf = maxpt/2+1
|
|
nrx = maxbeg+ihalf
|
|
next = number+1
|
|
if(next.gt.nrint) go to 40
|
|
c adjust the different parameters.
|
|
do 30 j=next,nrint
|
|
jj = next+nrint-j
|
|
fpint(jj+1) = fpint(jj)
|
|
nrdata(jj+1) = nrdata(jj)
|
|
jk = jj+k
|
|
t(jk+1) = t(jk)
|
|
30 continue
|
|
40 nrdata(number) = ihalf-1
|
|
nrdata(next) = maxpt-ihalf
|
|
am = maxpt
|
|
an = nrdata(number)
|
|
fpint(number) = fpmax*an/am
|
|
an = nrdata(next)
|
|
fpint(next) = fpmax*an/am
|
|
jk = next+k
|
|
t(jk) = x(nrx)
|
|
50 n = n+1
|
|
nrint = nrint+1
|
|
return
|
|
end
|