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
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
|
subroutine fppara(iopt,idim,m,u,mx,x,w,ub,ue,k,s,nest,tol,maxit,
* k1,k2,n,t,nc,c,fp,fpint,z,a,b,g,q,nrdata,ier)
c ..
c ..scalar arguments..
real*8 ub,ue,s,tol,fp
integer iopt,idim,m,mx,k,nest,maxit,k1,k2,n,nc,ier
c ..array arguments..
real*8 u(m),x(mx),w(m),t(nest),c(nc),fpint(nest),
* z(nc),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,fac,fpart,fpms,fpold,fp0,f1,f2,f3,
* half,one,p,pinv,piv,p1,p2,p3,rn,sin,store,term,ui,wi
integer i,ich1,ich3,it,iter,i1,i2,i3,j,jj,j1,j2,k3,l,l0,
* mk1,new,nk1,nmax,nmin,nplus,npl1,nrint,n8
c ..local arrays..
real*8 h(7),xi(10)
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.1e+01
con1 = 0.1e0
con9 = 0.9e0
con4 = 0.4e-01
half = 0.5e0
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 curve sinf(u), c
c and the corresponding sum of squared residuals fp=f(p=inf). c
c if iopt=-1 sinf(u) is the requested curve. 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 curve 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 curve 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 polynomial curve 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.) go to 45
c if s=0, s(u) is an interpolating curve.
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) = u(j)
i = i+1
j = j+1
20 continue
go to 60
30 do 40 l=1,mk1
t(i) = (u(j)+u(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 curve which is a spline curve without interior knots.
c if iopt=1 and fp0>s we start computing the least squares spline curve
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.
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(u).
nk1 = n-k1
i = n
do 70 j=1,k1
t(j) = ub
t(i) = ue
i = i-1
70 continue
c compute the b-spline coefficients of the least-squares spline curve
c sinf(u). 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.
c initialize the b-spline coefficients and the observation matrix a.
do 75 i=1,nc
z(i) = 0.
75 continue
do 80 i=1,nk1
do 80 j=1,k1
a(i,j) = 0.
80 continue
l = k1
jj = 0
do 130 it=1,m
c fetch the current data point u(it),x(it).
ui = u(it)
wi = w(it)
do 83 j=1,idim
jj = jj+1
xi(j) = x(jj)*wi
83 continue
c search for knot interval t(l) <= ui < t(l+1).
85 if(ui.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 ui and store them in q.
90 call fpbspl(t,n,k,ui,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.) 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.
j1 = j
do 97 j2 =1,idim
call fprota(cos,sin,xi(j2),z(j1))
j1 = j1+n
97 continue
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 do 125 j2=1,idim
fp = fp+xi(j2)**2
125 continue
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.
j1 = 1
do 135 j2=1,idim
call fpback(a,z(j1),nk1,k1,c(j1),nest)
j1 = j1+n
135 continue
c test whether the approximation sinf(u) 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.) go to 250
c if n = nmax, sinf(u) is an interpolating spline curve.
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 of squared residuals for each knot interval
c t(j+k) <= u(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
fpart = 0.
i = 1
l = k2
new = 0
jj = 0
do 180 it=1,m
if(u(it).lt.t(l) .or. l.gt.nk1) go to 160
new = 1
l = l+1
160 term = 0.
l0 = l-k2
do 175 j2=1,idim
fac = 0.
j1 = l0
do 170 j=1,k1
j1 = j1+1
fac = fac+c(j1)*q(it,j)
170 continue
jj = jj+1
term = term+(w(it)*(fac-x(jj)))**2
l0 = l0+n
175 continue
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(u,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 curve is a
c solution of our approximation problem.
250 if(ier.eq.(-2)) go to 440
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c part 2: determination of the smoothing spline curve sp(u). 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 curve c
c sp(u). the observation matrix a is extended by the rows of matrix c
c b expressing that the kth derivative discontinuities of sp(u) 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 f(p),c
c the sum of squared residuals be = s. we already know that the least c
c squares kth degree polynomial curve corresponds to p=0, and that c
c the least-squares spline curve 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.
f1 = fp0-s
p3 = -one
f3 = fpms
p = 0.
do 252 i=1,nk1
p = p+a(i,1)
252 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 255 i=1,nc
c(i) = z(i)
255 continue
do 260 i=1,nk1
g(i,k2) = 0.
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
do 275 j=1,idim
xi(j) = 0.
275 continue
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.
j1 = j
do 277 j2=1,idim
call fprota(cos,sin,xi(j2),c(j1))
j1 = j1+n
277 continue
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.
290 continue
300 continue
c backward substitution to obtain the b-spline coefficients.
j1 = 1
do 305 j2=1,idim
call fpback(g,c(j1),nk1,k2,c(j1),nest)
j1 =j1+n
305 continue
c computation of f(p).
fp = 0.
l = k2
jj = 0
do 330 it=1,m
if(u(it).lt.t(l) .or. l.gt.nk1) go to 310
l = l+1
310 l0 = l-k2
term = 0.
do 325 j2=1,idim
fac = 0.
j1 = l0
do 320 j=1,k1
j1 = j1+1
fac = fac+c(j1)*q(it,j)
320 continue
jj = jj+1
term = term+(fac-x(jj))**2
l0 = l0+n
325 continue
fp = fp+term*w(it)**2
330 continue
c test whether the approximation sp(u) 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.) 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.) 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
|