1
subroutine fpperi(iopt,x,y,w,m,k,s,nest,tol,maxit,k1,k2,n,t,c,
2
* fp,fpint,z,a1,a2,b,g1,g2,q,nrdata,ier)
6
integer iopt,m,k,nest,maxit,k1,k2,n,ier
8
real*8 x(m),y(m),w(m),t(nest),c(nest),fpint(nest),z(nest),
9
* a1(nest,k1),a2(nest,k),b(nest,k2),g1(nest,k2),g2(nest,k1),
13
real*8 acc,cos,c1,d1,fpart,fpms,fpold,fp0,f1,f2,f3,p,per,pinv,piv,
15
* p1,p2,p3,sin,store,term,wi,xi,yi,rn,one,con1,con4,con9,half
16
integer i,ich1,ich3,ij,ik,it,iter,i1,i2,i3,j,jk,jper,j1,j2,kk,
17
* kk1,k3,l,l0,l1,l5,mm,m1,new,nk1,nk2,nmax,nmin,nplus,npl1,
20
real*8 h(6),h1(7),h2(6)
21
c ..function references..
24
c ..subroutine references..
25
c fpbacp,fpbspl,fpgivs,fpdisc,fpknot,fprota
33
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
34
c part 1: determination of the number of knots and their position c
35
c ************************************************************** c
36
c given a set of knots we compute the least-squares periodic spline c
37
c sinf(x). if the sum f(p=inf) <= s we accept the choice of knots. c
38
c the initial choice of knots depends on the value of s and iopt. c
39
c if s=0 we have spline interpolation; in that case the number of c
40
c knots equals nmax = m+2*k. c
42
c iopt=0 we first compute the least-squares polynomial of c
43
c degree k; n = nmin = 2*k+2. since s(x) must be periodic we c
44
c find that s(x) is a constant function. c
45
c iopt=1 we start with the set of knots found at the last c
46
c call of the routine, except for the case that s > fp0; then c
47
c we compute directly the least-squares periodic polynomial. c
48
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
54
c determine the length of the period of s(x).
56
if(iopt.lt.0) go to 50
57
c calculation of acc, the absolute tolerance for the root of f(p)=s.
59
c determine nmax, the number of knots for periodic spline interpolation
61
if(s.gt.0. .or. nmax.eq.nmin) go to 30
62
c if s=0, s(x) is an interpolating spline.
64
c test whether the required storage space exceeds the available one.
65
if(n.gt.nest) go to 620
66
c find the position of the interior knots in case of interpolation.
67
5 if((k/2)*2 .eq. k) go to 20
91
t(j) = (x(i)+x(i-1))*half
94
c if s > 0 our initial choice depends on the value of iopt.
95
c if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
96
c periodic polynomial. (i.e. a constant function).
97
c if iopt=1 and fp0>s we start computing the least-squares periodic
98
c spline according the set of knots found at the last call of the
100
30 if(iopt.eq.0) go to 35
101
if(n.eq.nmin) go to 35
105
if(fp0.gt.s) go to 50
106
c the case that s(x) is a constant function is treated separetely.
107
c find the least-squares constant c1 and compute fp0 at the same time.
114
call fpgivs(wi,d1,cos,sin)
115
call fprota(cos,sin,yi,c1)
119
c test whether that constant function is a solution of our problem.
121
if(fpms.lt.acc .or. nmax.eq.nmin) go to 640
123
c test whether the required storage space exceeds the available one.
124
if(nmin.ge.nest) go to 620
125
c start computing the least-squares periodic spline with one
133
c main loop for the different sets of knots. m is a save upper
134
c bound for the number of trials.
136
c find nrint, the number of knot intervals.
138
c find the position of the additional knots which are needed for
139
c the b-spline representation of s(x). if we take
140
c t(k+1) = x(1), t(n-k) = x(m)
141
c t(k+1-j) = t(n-k-j) - per, j=1,2,...k
142
c t(n-k+j) = t(k+1+j) + per, j=1,2,...k
143
c then s(x) is a periodic spline with period per if the b-spline
144
c coefficients satisfy the following conditions
145
c c(n7+j) = c(j), j=1,...k (**) with n7=n-2*k-1.
158
c compute the b-spline coefficients c(j),j=1,...n7 of the least-squares
159
c periodic spline sinf(x). the observation matrix a is built up row
160
c by row while taking into account condition (**) and is reduced to
161
c triangular form by givens transformations .
162
c at the same time fp=f(p=inf) is computed.
163
c the n7 x n7 triangularised upper matrix a has the form
167
c with a2 a n7 x k matrix and a1 a n10 x n10 upper triangular
168
c matrix of bandwith k+1 ( n10 = n7-k).
181
c fetch the current data point x(it),y(it)
185
c search for knot interval t(l) <= xi < t(l+1).
186
80 if(xi.lt.t(l+1)) go to 85
189
c evaluate the (k+1) non-zero b-splines at xi and store them in q.
190
85 call fpbspl(t,n,k,xi,l,h)
196
c test whether the b-splines nj,k+1(x),j=1+n7,...nk1 are all zero at xi
197
if(l5.lt.n10) go to 285
198
if(jper.ne.0) go to 160
199
c initialize the matrix a2.
208
if(ik.le.0) go to 105
215
c if one of the b-splines nj,k+1(x),j=n7+1,...nk1 is not zero at xi
216
c we take account of condition (**) for setting up the new row
217
c of the observation matrix a. this row is stored in the arrays h1
218
c (the part with respect to a1) and h2 (the part with
230
if(l1.le.0) go to 200
231
if(l1.le.n10) go to 190
236
200 h2(l0) = h2(l0)+h(i)
238
c rotate the new row of the observation matrix into triangle
239
c by givens transformations.
240
if(n10.le.0) go to 250
241
c rotation with the rows 1,2,...n10 of matrix a.
244
if(piv.ne.0.) go to 214
250
c calculate the parameters of the givens transformation.
251
214 call fpgivs(piv,a1(j,1),cos,sin)
252
c transformation to the right hand side.
253
call fprota(cos,sin,yi,z(j))
254
c transformations to the left hand side with respect to a2.
256
call fprota(cos,sin,h2(i),a2(j,i))
258
if(j.eq.n10) go to 250
260
c transformations to the left hand side with respect to a1.
263
call fprota(cos,sin,h1(i1),a1(j,i1))
268
c rotation with the rows n10+1,...n7 of matrix a.
271
if(ij.le.0) go to 270
273
if(piv.eq.0.) go to 270
274
c calculate the parameters of the givens transformation.
275
call fpgivs(piv,a2(ij,j),cos,sin)
276
c transformations to right hand side.
277
call fprota(cos,sin,yi,z(ij))
278
if(j.eq.kk) go to 280
280
c transformations to left hand side.
282
call fprota(cos,sin,h2(i),a2(ij,i))
285
c add contribution of this row to the sum of squares of residual
289
c rotation of the new row of the observation matrix into
290
c triangle in case the b-splines nj,k+1(x),j=n7+1,...n-k-1 are all zero
296
if(piv.eq.0.) go to 140
297
c calculate the parameters of the givens transformation.
298
call fpgivs(piv,a1(j,1),cos,sin)
299
c transformations to right hand side.
300
call fprota(cos,sin,yi,z(j))
301
if(i.eq.kk1) go to 150
304
c transformations to left hand side.
307
call fprota(cos,sin,h(i1),a1(j,i2))
310
c add contribution of this row to the sum of squares of residual
317
c backward substitution to obtain the b-spline coefficients c(j),j=1,.n
318
call fpbacp(a1,a2,z,n7,kk,c,kk1,nest)
319
c calculate from condition (**) the coefficients c(j+n7),j=1,2,...k.
324
if(iopt.lt.0) go to 660
325
c test whether the approximation sinf(x) is an acceptable solution.
327
if(abs(fpms).lt.acc) go to 660
328
c if f(p=inf) < s accept the choice of knots.
329
if(fpms.lt.0.) go to 350
330
c if n=nmax, sinf(x) is an interpolating spline.
331
if(n.eq.nmax) go to 630
332
c increase the number of knots.
333
c if n=nest we cannot increase the number of knots because of the
334
c storage capacity limitation.
335
if(n.eq.nest) go to 620
336
c determine the number of knots nplus we are going to add.
339
if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
340
nplus = min0(nplus*2,max0(npl1,nplus/2,1))
342
c compute the sum(wi*(yi-s(xi))**2) for each knot interval
343
c t(j+k) <= xi <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
348
if(x(it).lt.t(l)) go to 300
355
term = term+c(l0)*q(it,j)
357
term = (w(it)*(term-y(it)))**2
359
if(new.eq.0) go to 320
360
if(l.gt.k2) go to 315
364
315 store = term*half
365
fpint(i) = fpart-store
370
fpint(nrint) = fpint(nrint)+fpart
373
call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
374
c if n=nmax we locate the knots as for interpolation.
375
if(n.eq.nmax) go to 5
376
c test whether we cannot further increase the number of knots.
377
if(n.eq.nest) go to 340
379
c restart the computations with the new set of knots.
381
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
382
c part 2: determination of the smoothing periodic spline sp(x). c
383
c ************************************************************* c
384
c we have determined the number of knots and their position. c
385
c we now compute the b-spline coefficients of the smoothing spline c
386
c sp(x). the observation matrix a is extended by the rows of matrix c
387
c b expressing that the kth derivative discontinuities of sp(x) at c
388
c the interior knots t(k+2),...t(n-k-1) must be zero. the corres- c
389
c ponding weights of these additional rows are set to 1/sqrt(p). c
390
c iteratively we then have to determine the value of p such that c
391
c f(p)=sum(w(i)*(y(i)-sp(x(i)))**2) be = s. we already know that c
392
c the least-squares constant function corresponds to p=0, and that c
393
c the least-squares periodic spline corresponds to p=infinity. the c
394
c iteration process which is proposed here, makes use of rational c
395
c interpolation. since f(p) is a convex and strictly decreasing c
396
c function of p, it can be approximated by a rational function c
397
c r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond- c
398
c ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used c
399
c to calculate the new value of p such that r(p)=s. convergence is c
400
c guaranteed by taking f1>0 and f3<0. c
401
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
402
c evaluate the discontinuity jump of the kth derivative of the
403
c b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
404
350 call fpdisc(t,n,k2,b,nest)
405
c initial value for p.
427
c iteration process to find the root of f(p) = s.
429
c form the matrix g as the matrix a extended by the rows of matrix b.
430
c the rows of matrix b with weight 1/p are rotated into
431
c the triangularised observation matrix a.
432
c after triangularisation our n7 x n7 matrix g takes the form
436
c with g2 a n7 x (k+1) matrix and g1 a n11 x n11 upper triangular
437
c matrix of bandwidth k+2. ( n11 = n7-k-1)
439
c store matrix a into g
456
c fetch a new row of matrix b and store it in the arrays h1 (the part
457
c with respect to g1) and h2 (the part with respect to g2).
464
if(it.gt.n11) go to 420
468
if(l0.eq.n10) go to 400
475
h2(l0) = b(it,l1)*pinv
485
if(l1.le.0) go to 450
486
if(l1.le.n11) go to 440
489
440 h1(l1) = b(it,j)*pinv
491
450 h2(l0) = h2(l0)+b(it,j)*pinv
493
if(n11.le.0) go to 510
494
c rotate this row into triangle by givens transformations without
496
c rotation with the rows l,l+1,...n11.
499
c calculate the parameters of the givens transformation.
500
call fpgivs(piv,g1(j,1),cos,sin)
501
c transformation to right hand side.
502
call fprota(cos,sin,yi,c(j))
503
c transformation to the left hand side with respect to g2.
505
call fprota(cos,sin,h2(i),g2(j,i))
507
if(j.eq.n11) go to 510
509
c transformation to the left hand side with respect to g1.
512
call fprota(cos,sin,h1(i1),g1(j,i1))
517
c rotation with the rows n11+1,...n7
520
if(ij.le.0) go to 530
522
c calculate the parameters of the givens transformation
523
call fpgivs(piv,g2(ij,j),cos,sin)
524
c transformation to the right hand side.
525
call fprota(cos,sin,yi,c(ij))
526
if(j.eq.k1) go to 540
528
c transformation to the left hand side.
530
call fprota(cos,sin,h2(i),g2(ij,i))
534
c backward substitution to obtain the b-spline coefficients
535
c c(j),j=1,2,...n7 of sp(x).
536
call fpbacp(g1,g2,c,n7,k1,c,k2,nest)
537
c calculate from condition (**) the b-spline coefficients c(n7+j),j=1,.
542
c computation of f(p).
546
if(x(it).lt.t(l)) go to 550
552
term = term+c(l0)*q(it,j)
554
fp = fp+(w(it)*(term-y(it)))**2
556
c test whether the approximation sp(x) is an acceptable solution.
558
if(abs(fpms).lt.acc) go to 660
559
c test whether the maximal number of iterations is reached.
560
if(iter.eq.maxit) go to 600
561
c carry out one more step of the iteration process.
564
if(ich3.ne.0) go to 580
565
if((f2-f3) .gt. acc) go to 575
566
c our initial choice of p is too large.
570
if(p.le.p1) p = p1*con9 +p2*con1
572
575 if(f2.lt.0.) ich3 = 1
573
580 if(ich1.ne.0) go to 590
574
if((f1-f2) .gt. acc) go to 585
575
c our initial choice of p is too small
579
if(p3.lt.0.) go to 595
580
if(p.ge.p3) p = p2*con1 +p3*con9
582
585 if(f2.gt.0.) ich1 = 1
583
c test whether the iteration process proceeds as theoretically
585
590 if(f2.ge.f1 .or. f2.le.f3) go to 610
586
c find the new value for p.
587
p = fprati(p1,f1,p2,f2,p3,f3)
589
c error codes and messages.
599
c the least-squares constant function c1 is a solution of our problem.
600
c a constant function is a spline of degree k with all b-spline
601
c coefficients equal to that constant c1.
added
removed
Lib/sandbox/spline/fitpack/fpperi.f
