1
subroutine fppara(iopt,idim,m,u,mx,x,w,ub,ue,k,s,nest,tol,maxit,
2
* k1,k2,n,t,nc,c,fp,fpint,z,a,b,g,q,nrdata,ier)
6
integer iopt,idim,m,mx,k,nest,maxit,k1,k2,n,nc,ier
8
real*8 u(m),x(mx),w(m),t(nest),c(nc),fpint(nest),
9
* z(nc),a(nest,k1),b(nest,k2),g(nest,k2),q(m,k1)
12
real*8 acc,con1,con4,con9,cos,fac,fpart,fpms,fpold,fp0,f1,f2,f3,
13
* half,one,p,pinv,piv,p1,p2,p3,rn,sin,store,term,ui,wi
14
integer i,ich1,ich3,it,iter,i1,i2,i3,j,jj,j1,j2,k3,l,l0,
15
* mk1,new,nk1,nmax,nmin,nplus,npl1,nrint,n8
18
c ..function references
21
c ..subroutine references..
22
c fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
30
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
31
c part 1: determination of the number of knots and their position c
32
c ************************************************************** c
33
c given a set of knots we compute the least-squares curve sinf(u), c
34
c and the corresponding sum of squared residuals fp=f(p=inf). c
35
c if iopt=-1 sinf(u) is the requested curve. c
36
c if iopt=0 or iopt=1 we check whether we can accept the knots: c
37
c if fp <=s we will continue with the current set of knots. c
38
c if fp > s we will increase the number of knots and compute the c
39
c corresponding least-squares curve until finally fp<=s. c
40
c the initial choice of knots depends on the value of s and iopt. c
41
c if s=0 we have spline interpolation; in that case the number of c
42
c knots equals nmax = m+k+1. c
44
c iopt=0 we first compute the least-squares polynomial curve of c
45
c degree k; n = nmin = 2*k+2 c
46
c iopt=1 we start with the set of knots found at the last c
47
c call of the routine, except for the case that s > fp0; then c
48
c we compute directly the polynomial curve of degree k. c
49
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
50
c determine nmin, the number of knots for polynomial approximation.
52
if(iopt.lt.0) go to 60
53
c calculation of acc, the absolute tolerance for the root of f(p)=s.
55
c determine nmax, the number of knots for spline interpolation.
58
c if s=0, s(u) is an interpolating curve.
59
c test whether the required storage space exceeds the available one.
61
if(nmax.gt.nest) go to 420
62
c find the position of the interior knots in case of interpolation.
68
if(k3*2.eq.k) go to 30
76
t(i) = (u(j)+u(j-1))*half
81
c if s>0 our initial choice of knots depends on the value of iopt.
82
c if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
83
c polynomial curve which is a spline curve without interior knots.
84
c if iopt=1 and fp0>s we start computing the least squares spline curve
85
c according to the set of knots found at the last call of the routine.
86
45 if(iopt.eq.0) go to 50
87
if(n.eq.nmin) go to 50
96
c main loop for the different sets of knots. m is a save upper bound
97
c for the number of trials.
99
if(n.eq.nmin) ier = -2
100
c find nrint, tne number of knot intervals.
102
c find the position of the additional knots which are needed for
103
c the b-spline representation of s(u).
111
c compute the b-spline coefficients of the least-squares spline curve
112
c sinf(u). the observation matrix a is built up row by row and
113
c reduced to upper triangular form by givens transformations.
114
c at the same time fp=f(p=inf) is computed.
116
c initialize the b-spline coefficients and the observation matrix a.
127
c fetch the current data point u(it),x(it).
134
c search for knot interval t(l) <= ui < t(l+1).
135
85 if(ui.lt.t(l+1) .or. l.eq.nk1) go to 90
138
c evaluate the (k+1) non-zero b-splines at ui and store them in q.
139
90 call fpbspl(t,n,k,ui,l,h)
144
c rotate the new row of the observation matrix into triangle.
149
if(piv.eq.0.) go to 110
150
c calculate the parameters of the givens transformation.
151
call fpgivs(piv,a(j,1),cos,sin)
152
c transformations to right hand side.
155
call fprota(cos,sin,xi(j2),z(j1))
158
if(i.eq.k1) go to 120
163
c transformations to left hand side.
164
call fprota(cos,sin,h(i1),a(j,i2))
167
c add contribution of this row to the sum of squares of residual
173
if(ier.eq.(-2)) fp0 = fp
177
c backward substitution to obtain the b-spline coefficients.
180
call fpback(a,z(j1),nk1,k1,c(j1),nest)
183
c test whether the approximation sinf(u) is an acceptable solution.
184
if(iopt.lt.0) go to 440
186
if(abs(fpms).lt.acc) go to 440
187
c if f(p=inf) < s accept the choice of knots.
188
if(fpms.lt.0.) go to 250
189
c if n = nmax, sinf(u) is an interpolating spline curve.
190
if(n.eq.nmax) go to 430
191
c increase the number of knots.
192
c if n=nest we cannot increase the number of knots because of
193
c the storage capacity limitation.
194
if(n.eq.nest) go to 420
195
c determine the number of knots nplus we are going to add.
196
if(ier.eq.0) go to 140
202
if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
203
nplus = min0(nplus*2,max0(npl1,nplus/2,1))
205
c compute the sum of squared residuals for each knot interval
206
c t(j+k) <= u(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
213
if(u(it).lt.t(l) .or. l.gt.nk1) go to 160
223
fac = fac+c(j1)*q(it,j)
226
term = term+(w(it)*(fac-x(jj)))**2
230
if(new.eq.0) go to 180
232
fpint(i) = fpart-store
240
call fpknot(u,m,t,n,fpint,nrdata,nrint,nest,1)
241
c if n=nmax we locate the knots as for interpolation
242
if(n.eq.nmax) go to 10
243
c test whether we cannot further increase the number of knots.
244
if(n.eq.nest) go to 200
246
c restart the computations with the new set of knots.
248
c test whether the least-squares kth degree polynomial curve is a
249
c solution of our approximation problem.
250
250 if(ier.eq.(-2)) go to 440
251
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
252
c part 2: determination of the smoothing spline curve sp(u). c
253
c ********************************************************** c
254
c we have determined the number of knots and their position. c
255
c we now compute the b-spline coefficients of the smoothing curve c
256
c sp(u). the observation matrix a is extended by the rows of matrix c
257
c b expressing that the kth derivative discontinuities of sp(u) at c
258
c the interior knots t(k+2),...t(n-k-1) must be zero. the corres- c
259
c ponding weights of these additional rows are set to 1/p. c
260
c iteratively we then have to determine the value of p such that f(p),c
261
c the sum of squared residuals be = s. we already know that the least c
262
c squares kth degree polynomial curve corresponds to p=0, and that c
263
c the least-squares spline curve corresponds to p=infinity. the c
264
c iteration process which is proposed here, makes use of rational c
265
c interpolation. since f(p) is a convex and strictly decreasing c
266
c function of p, it can be approximated by a rational function c
267
c r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond- c
268
c ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used c
269
c to calculate the new value of p such that r(p)=s. convergence is c
270
c guaranteed by taking f1>0 and f3<0. c
271
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
272
c evaluate the discontinuity jump of the kth derivative of the
273
c b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
274
call fpdisc(t,n,k2,b,nest)
275
c initial value for p.
289
c iteration process to find the root of f(p) = s.
291
c the rows of matrix b with weight 1/p are rotated into the
292
c triangularised observation matrix a which is stored in g.
303
c the row of matrix b is rotated into triangle by givens transformation
312
c calculate the parameters of the givens transformation.
313
call fpgivs(piv,g(j,1),cos,sin)
314
c transformations to right hand side.
317
call fprota(cos,sin,xi(j2),c(j1))
320
if(j.eq.nk1) go to 300
322
if(j.gt.n8) i2 = nk1-j
324
c transformations to left hand side.
326
call fprota(cos,sin,h(i1),g(j,i1))
332
c backward substitution to obtain the b-spline coefficients.
335
call fpback(g,c(j1),nk1,k2,c(j1),nest)
338
c computation of f(p).
343
if(u(it).lt.t(l) .or. l.gt.nk1) go to 310
352
fac = fac+c(j1)*q(it,j)
355
term = term+(fac-x(jj))**2
358
fp = fp+term*w(it)**2
360
c test whether the approximation sp(u) is an acceptable solution.
362
if(abs(fpms).lt.acc) go to 440
363
c test whether the maximal number of iterations is reached.
364
if(iter.eq.maxit) go to 400
365
c carry out one more step of the iteration process.
368
if(ich3.ne.0) go to 340
369
if((f2-f3).gt.acc) go to 335
370
c our initial choice of p is too large.
374
if(p.le.p1) p=p1*con9 + p2*con1
376
335 if(f2.lt.0.) ich3=1
377
340 if(ich1.ne.0) go to 350
378
if((f1-f2).gt.acc) go to 345
379
c our initial choice of p is too small
383
if(p3.lt.0.) go to 360
384
if(p.ge.p3) p = p2*con1 + p3*con9
386
345 if(f2.gt.0.) ich1=1
387
c test whether the iteration process proceeds as theoretically
389
350 if(f2.ge.f1 .or. f2.le.f3) go to 410
390
c find the new value for p.
391
p = fprati(p1,f1,p2,f2,p3,f3)
393
c error codes and messages.
added
removed
Lib/sandbox/spline/fitpack/fppara.f
