1
subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,
3
c given the set of data points (x(i),y(i)) and the set of positive
4
c numbers w(i),i=1,2,...,m-1, subroutine percur determines a smooth
5
c periodic spline approximation of degree k with period per=x(m)-x(1).
6
c if iopt=-1 percur calculates the weighted least-squares periodic
7
c spline according to a given set of knots.
8
c if iopt>=0 the number of knots of the spline s(x) and the position
9
c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
10
c ness of s(x) is then achieved by minimalizing the discontinuity
11
c jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
12
c n-k-1. the amount of smoothness is determined by the condition that
13
c f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
14
c negative constant, called the smoothing factor.
15
c the fit s(x) is given in the b-spline representation (b-spline coef-
16
c ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
20
c call percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,
24
c iopt : integer flag. on entry iopt must specify whether a weighted
25
c least-squares spline (iopt=-1) or a smoothing spline (iopt=
26
c 0 or 1) must be determined. if iopt=0 the routine will start
27
c with an initial set of knots t(i)=x(1)+(x(m)-x(1))*(i-k-1),
28
c i=1,2,...,2*k+2. if iopt=1 the routine will continue with
29
c the knots found at the last call of the routine.
30
c attention: a call with iopt=1 must always be immediately
31
c preceded by another call with iopt=1 or iopt=0.
33
c m : integer. on entry m must specify the number of data points.
34
c m > 1. unchanged on exit.
35
c x : real array of dimension at least (m). before entry, x(i)
36
c must be set to the i-th value of the independent variable x,
37
c for i=1,2,...,m. these values must be supplied in strictly
38
c ascending order. x(m) only indicates the length of the
39
c period of the spline, i.e per=x(m)-x(1).
41
c y : real array of dimension at least (m). before entry, y(i)
42
c must be set to the i-th value of the dependent variable y,
43
c for i=1,2,...,m-1. the element y(m) is not used.
45
c w : real array of dimension at least (m). before entry, w(i)
46
c must be set to the i-th value in the set of weights. the
47
c w(i) must be strictly positive. w(m) is not used.
48
c see also further comments. unchanged on exit.
49
c k : integer. on entry k must specify the degree of the spline.
50
c 1<=k<=5. it is recommended to use cubic splines (k=3).
51
c the user is strongly dissuaded from choosing k even,together
52
c with a small s-value. unchanged on exit.
53
c s : real.on entry (in case iopt>=0) s must specify the smoothing
54
c factor. s >=0. unchanged on exit.
55
c for advice on the choice of s see further comments.
56
c nest : integer. on entry nest must contain an over-estimate of the
57
c total number of knots of the spline returned, to indicate
58
c the storage space available to the routine. nest >=2*k+2.
59
c in most practical situation nest=m/2 will be sufficient.
60
c always large enough is nest=m+2*k,the number of knots needed
61
c for interpolation (s=0). unchanged on exit.
63
c unless ier = 10 (in case iopt >=0), n will contain the
64
c total number of knots of the spline approximation returned.
65
c if the computation mode iopt=1 is used this value of n
66
c should be left unchanged between subsequent calls.
67
c in case iopt=-1, the value of n must be specified on entry.
68
c t : real array of dimension at least (nest).
69
c on succesful exit, this array will contain the knots of the
70
c spline,i.e. the position of the interior knots t(k+2),t(k+3)
71
c ...,t(n-k-1) as well as the position of the additional knots
72
c t(1),t(2),...,t(k+1)=x(1) and t(n-k)=x(m),..,t(n) needed for
73
c the b-spline representation.
74
c if the computation mode iopt=1 is used, the values of t(1),
75
c t(2),...,t(n) should be left unchanged between subsequent
76
c calls. if the computation mode iopt=-1 is used, the values
77
c t(k+2),...,t(n-k-1) must be supplied by the user, before
78
c entry. see also the restrictions (ier=10).
79
c c : real array of dimension at least (nest).
80
c on succesful exit, this array will contain the coefficients
81
c c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
82
c fp : real. unless ier = 10, fp contains the weighted sum of
83
c squared residuals of the spline approximation returned.
84
c wrk : real array of dimension at least (m*(k+1)+nest*(8+5*k)).
85
c used as working space. if the computation mode iopt=1 is
86
c used, the values wrk(1),...,wrk(n) should be left unchanged
87
c between subsequent calls.
88
c lwrk : integer. on entry,lwrk must specify the actual dimension of
89
c the array wrk as declared in the calling (sub)program. lwrk
90
c must not be too small (see wrk). unchanged on exit.
91
c iwrk : integer array of dimension at least (nest).
92
c used as working space. if the computation mode iopt=1 is
93
c used,the values iwrk(1),...,iwrk(n) should be left unchanged
94
c between subsequent calls.
95
c ier : integer. unless the routine detects an error, ier contains a
96
c non-positive value on exit, i.e.
97
c ier=0 : normal return. the spline returned has a residual sum of
98
c squares fp such that abs(fp-s)/s <= tol with tol a relat-
99
c ive tolerance set to 0.001 by the program.
100
c ier=-1 : normal return. the spline returned is an interpolating
101
c periodic spline (fp=0).
102
c ier=-2 : normal return. the spline returned is the weighted least-
103
c squares constant. in this extreme case fp gives the upper
104
c bound fp0 for the smoothing factor s.
105
c ier=1 : error. the required storage space exceeds the available
106
c storage space, as specified by the parameter nest.
107
c probably causes : nest too small. if nest is already
108
c large (say nest > m/2), it may also indicate that s is
110
c the approximation returned is the least-squares periodic
111
c spline according to the knots t(1),t(2),...,t(n). (n=nest)
112
c the parameter fp gives the corresponding weighted sum of
113
c squared residuals (fp>s).
114
c ier=2 : error. a theoretically impossible result was found during
115
c the iteration proces for finding a smoothing spline with
116
c fp = s. probably causes : s too small.
117
c there is an approximation returned but the corresponding
118
c weighted sum of squared residuals does not satisfy the
119
c condition abs(fp-s)/s < tol.
120
c ier=3 : error. the maximal number of iterations maxit (set to 20
121
c by the program) allowed for finding a smoothing spline
122
c with fp=s has been reached. probably causes : s too small
123
c there is an approximation returned but the corresponding
124
c weighted sum of squared residuals does not satisfy the
125
c condition abs(fp-s)/s < tol.
126
c ier=10 : error. on entry, the input data are controlled on validity
127
c the following restrictions must be satisfied.
128
c -1<=iopt<=1, 1<=k<=5, m>1, nest>2*k+2, w(i)>0,i=1,...,m-1
129
c x(1)<x(2)<...<x(m), lwrk>=(k+1)*m+nest*(8+5*k)
130
c if iopt=-1: 2*k+2<=n<=min(nest,m+2*k)
131
c x(1)<t(k+2)<t(k+3)<...<t(n-k-1)<x(m)
132
c the schoenberg-whitney conditions, i.e. there
133
c must be a subset of data points xx(j) with
134
c xx(j) = x(i) or x(i)+(x(m)-x(1)) such that
135
c t(j) < xx(j) < t(j+k+1), j=k+1,...,n-k-1
137
c if s=0 : nest >= m+2*k
138
c if one of these conditions is found to be violated,control
139
c is immediately repassed to the calling program. in that
140
c case there is no approximation returned.
143
c by means of the parameter s, the user can control the tradeoff
144
c between closeness of fit and smoothness of fit of the approximation.
145
c if s is too large, the spline will be too smooth and signal will be
146
c lost ; if s is too small the spline will pick up too much noise. in
147
c the extreme cases the program will return an interpolating periodic
148
c spline if s=0 and the weighted least-squares constant if s is very
149
c large. between these extremes, a properly chosen s will result in
150
c a good compromise between closeness of fit and smoothness of fit.
151
c to decide whether an approximation, corresponding to a certain s is
152
c satisfactory the user is highly recommended to inspect the fits
154
c recommended values for s depend on the weights w(i). if these are
155
c taken as 1/d(i) with d(i) an estimate of the standard deviation of
156
c y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
157
c sqrt(2*m)). if nothing is known about the statistical error in y(i)
158
c each w(i) can be set equal to one and s determined by trial and
159
c error, taking account of the comments above. the best is then to
160
c start with a very large value of s ( to determine the least-squares
161
c constant and the corresponding upper bound fp0 for s) and then to
162
c progressively decrease the value of s ( say by a factor 10 in the
163
c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
164
c approximation shows more detail) to obtain closer fits.
165
c to economize the search for a good s-value the program provides with
166
c different modes of computation. at the first call of the routine, or
167
c whenever he wants to restart with the initial set of knots the user
169
c if iopt=1 the program will continue with the set of knots found at
170
c the last call of the routine. this will save a lot of computation
171
c time if percur is called repeatedly for different values of s.
172
c the number of knots of the spline returned and their location will
173
c depend on the value of s and on the complexity of the shape of the
174
c function underlying the data. but, if the computation mode iopt=1
175
c is used, the knots returned may also depend on the s-values at
176
c previous calls (if these were smaller). therefore, if after a number
177
c of trials with different s-values and iopt=1, the user can finally
178
c accept a fit as satisfactory, it may be worthwhile for him to call
179
c percur once more with the selected value for s but now with iopt=0.
180
c indeed, percur may then return an approximation of the same quality
181
c of fit but with fewer knots and therefore better if data reduction
182
c is also an important objective for the user.
184
c other subroutines required:
185
c fpbacp,fpbspl,fpchep,fpperi,fpdisc,fpgivs,fpknot,fprati,fprota
188
c dierckx p. : algorithms for smoothing data with periodic and
189
c parametric splines, computer graphics and image
190
c processing 20 (1982) 171-184.
191
c dierckx p. : algorithms for smoothing data with periodic and param-
192
c etric splines, report tw55, dept. computer science,
194
c dierckx p. : curve and surface fitting with splines, monographs on
195
c numerical analysis, oxford university press, 1993.
199
c dept. computer science, k.u. leuven
200
c celestijnenlaan 200a, b-3001 heverlee, belgium.
201
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
203
c creation date : may 1979
204
c latest update : march 1987
207
c ..scalar arguments..
209
integer iopt,m,k,nest,n,lwrk,ier
210
c ..array arguments..
211
real*8 x(m),y(m),w(m),t(nest),c(nest),wrk(lwrk)
215
integer i,ia1,ia2,ib,ifp,ig1,ig2,iq,iz,i1,i2,j1,j2,k1,k2,lwest,
217
c ..subroutine references..
220
c we set up the parameters tol and maxit
223
c before starting computations a data check is made. if the input data
224
c are invalid, control is immediately repassed to the calling program.
226
if(k.le.0 .or. k.gt.5) go to 50
229
if(iopt.lt.(-1) .or. iopt.gt.1) go to 50
231
if(m.lt.2 .or. nest.lt.nmin) go to 50
232
lwest = m*k1+nest*(8+5*k)
233
if(lwrk.lt.lwest) go to 50
236
if(x(i).ge.x(i+1) .or. w(i).le.0.) go to 50
238
if(iopt.ge.0) go to 30
239
if(n.le.nmin .or. n.gt.nest) go to 50
255
call fpchep(x,m,t,n,k,ier)
256
if (ier.eq.0) go to 40
258
30 if(s.lt.0.) go to 50
259
if(s.eq.0. .and. nest.lt.(m+2*k)) go to 50
261
c we partition the working space and determine the spline approximation.
270
call fpperi(iopt,x,y,w,m,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,
271
* wrk(ifp),wrk(iz),wrk(ia1),wrk(ia2),wrk(ib),wrk(ig1),wrk(ig2),
added
removed
Lib/sandbox/spline/fitpack/percur.f
