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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
|
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)
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
|
