Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
!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..
double precision x(m),t(nest),fpint(nest)
integer nrdata(nest)
!c ..local scalars..
double precision an,am,fpmax
integer ihalf,j,jbegin,jj,jk,jpoint,k,maxbeg,maxpt
integer next,nrx,number
!c ..
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
fpmax = fpint(j)
number = j
maxpt = jpoint
maxbeg = jbegin
10 jbegin = jbegin+jpoint+1
20 continue
!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)
n = n+1
nrint = nrint+1
return
end