1
subroutine pogrid(iopt,ider,mu,u,mv,v,z,z0,r,s,nuest,nvest,
2
* nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
3
c subroutine pogrid fits a function f(x,y) to a set of data points
4
c z(i,j) given at the nodes (x,y)=(u(i)*cos(v(j)),u(i)*sin(v(j))),
5
c i=1,...,mu ; j=1,...,mv , of a radius-angle grid over a disc
6
c x ** 2 + y ** 2 <= r ** 2 .
8
c this approximation problem is reduced to the determination of a
9
c bicubic spline s(u,v) smoothing the data (u(i),v(j),z(i,j)) on the
10
c rectangle 0<=u<=r, v(1)<=v<=v(1)+2*pi
11
c in order to have continuous partial derivatives
18
c s(u,v)=f(x,y) must satisfy the following conditions
20
c (1) s(0,v) = g(0,0) v(1)<=v<= v(1)+2*pi
23
c (2) -------- = cos(v)*g(1,0)+sin(v)*g(0,1) v(1)<=v<= v(1)+2*pi
26
c moreover, s(u,v) must be periodic in the variable v, i.e.
30
c (3) ---------- = --------- 0 <=u<= r, j=0,1,2 , vb=v(1),
34
c the number of knots of s(u,v) and their position tu(i),i=1,2,...,nu;
35
c tv(j),j=1,2,...,nv, is chosen automatically by the routine. the
36
c smoothness of s(u,v) is achieved by minimalizing the discontinuity
37
c jumps of the derivatives of the spline at the knots. the amount of
38
c smoothness of s(u,v) is determined by the condition that
39
c fp=sumi=1,mu(sumj=1,mv((z(i,j)-s(u(i),v(j)))**2))+(z0-g(0,0))**2<=s,
40
c with s a given non-negative constant.
41
c the fit s(u,v) is given in its b-spline representation and can be
42
c evaluated by means of routine bispev. f(x,y) = s(u,v) can also be
43
c evaluated by means of function program evapol.
46
c call pogrid(iopt,ider,mu,u,mv,v,z,z0,r,s,nuest,nvest,nu,tu,
47
c * ,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
50
c iopt : integer array of dimension 3, specifying different options.
52
c iopt(1):on entry iopt(1) must specify whether a least-squares spline
53
c (iopt(1)=-1) or a smoothing spline (iopt(1)=0 or 1) must be
55
c if iopt(1)=0 the routine will start with an initial set of
56
c knots tu(i)=0,tu(i+4)=r,i=1,...,4;tv(i)=v(1)+(i-4)*2*pi,i=1,.
58
c if iopt(1)=1 the routine will continue with the set of knots
59
c found at the last call of the routine.
60
c attention: a call with iopt(1)=1 must always be immediately
61
c preceded by another call with iopt(1) = 1 or iopt(1) = 0.
62
c iopt(2):on entry iopt(2) must specify the requested order of conti-
63
c nuity for f(x,y) at the origin.
64
c if iopt(2)=0 only condition (1) must be fulfilled and
65
c if iopt(2)=1 conditions (1)+(2) must be fulfilled.
66
c iopt(3):on entry iopt(3) must specify whether (iopt(3)=1) or not
67
c (iopt(3)=0) the approximation f(x,y) must vanish at the
68
c boundary of the approximation domain.
69
c ider : integer array of dimension 2, specifying different options.
71
c ider(1):on entry ider(1) must specify whether (ider(1)=0 or 1) or not
72
c (ider(1)=-1) there is a data value z0 at the origin.
73
c if ider(1)=1, z0 will be considered to be the right function
74
c value, and it will be fitted exactly (g(0,0)=z0=c(1)).
75
c if ider(1)=0, z0 will be considered to be a data value just
76
c like the other data values z(i,j).
77
c ider(2):on entry ider(2) must specify whether (ider(2)=1) or not
78
c (ider(2)=0) f(x,y) must have vanishing partial derivatives
79
c g(1,0) and g(0,1) at the origin. (in case iopt(2)=1)
80
c mu : integer. on entry mu must specify the number of grid points
81
c along the u-axis. unchanged on exit.
82
c mu >= mumin where mumin=4-iopt(3)-ider(2) if ider(1)<0
83
c =3-iopt(3)-ider(2) if ider(1)>=0
84
c u : real array of dimension at least (mu). before entry, u(i)
85
c must be set to the u-co-ordinate of the i-th grid point
86
c along the u-axis, for i=1,2,...,mu. these values must be
87
c positive and supplied in strictly ascending order.
89
c mv : integer. on entry mv must specify the number of grid points
90
c along the v-axis. mv > 3 . unchanged on exit.
91
c v : real array of dimension at least (mv). before entry, v(j)
92
c must be set to the v-co-ordinate of the j-th grid point
93
c along the v-axis, for j=1,2,...,mv. these values must be
94
c supplied in strictly ascending order. unchanged on exit.
95
c -pi <= v(1) < pi , v(mv) < v(1)+2*pi.
96
c z : real array of dimension at least (mu*mv).
97
c before entry, z(mv*(i-1)+j) must be set to the data value at
98
c the grid point (u(i),v(j)) for i=1,...,mu and j=1,...,mv.
100
c z0 : real value. on entry (if ider(1) >=0 ) z0 must specify the
101
c data value at the origin. unchanged on exit.
102
c r : real value. on entry r must specify the radius of the disk.
103
c r>=u(mu) (>u(mu) if iopt(3)=1). unchanged on exit.
104
c s : real. on entry (if iopt(1)>=0) s must specify the smoothing
105
c factor. s >=0. unchanged on exit.
106
c for advice on the choice of s see further comments
107
c nuest : integer. unchanged on exit.
108
c nvest : integer. unchanged on exit.
109
c on entry, nuest and nvest must specify an upper bound for the
110
c number of knots required in the u- and v-directions respect.
111
c these numbers will also determine the storage space needed by
112
c the routine. nuest >= 8, nvest >= 8.
113
c in most practical situation nuest = mu/2, nvest=mv/2, will
114
c be sufficient. always large enough are nuest=mu+5+iopt(2)+
115
c iopt(3), nvest = mv+7, the number of knots needed for
116
c interpolation (s=0). see also further comments.
118
c unless ier=10 (in case iopt(1)>=0), nu will contain the total
119
c number of knots with respect to the u-variable, of the spline
120
c approximation returned. if the computation mode iopt(1)=1 is
121
c used, the value of nu should be left unchanged between sub-
122
c sequent calls. in case iopt(1)=-1, the value of nu should be
123
c specified on entry.
124
c tu : real array of dimension at least (nuest).
125
c on succesful exit, this array will contain the knots of the
126
c spline with respect to the u-variable, i.e. the position of
127
c the interior knots tu(5),...,tu(nu-4) as well as the position
128
c of the additional knots tu(1)=...=tu(4)=0 and tu(nu-3)=...=
129
c tu(nu)=r needed for the b-spline representation.
130
c if the computation mode iopt(1)=1 is used,the values of tu(1)
131
c ...,tu(nu) should be left unchanged between subsequent calls.
132
c if the computation mode iopt(1)=-1 is used, the values tu(5),
133
c ...tu(nu-4) must be supplied by the user, before entry.
134
c see also the restrictions (ier=10).
136
c unless ier=10 (in case iopt(1)>=0), nv will contain the total
137
c number of knots with respect to the v-variable, of the spline
138
c approximation returned. if the computation mode iopt(1)=1 is
139
c used, the value of nv should be left unchanged between sub-
140
c sequent calls. in case iopt(1) = -1, the value of nv should
141
c be specified on entry.
142
c tv : real array of dimension at least (nvest).
143
c on succesful exit, this array will contain the knots of the
144
c spline with respect to the v-variable, i.e. the position of
145
c the interior knots tv(5),...,tv(nv-4) as well as the position
146
c of the additional knots tv(1),...,tv(4) and tv(nv-3),...,
147
c tv(nv) needed for the b-spline representation.
148
c if the computation mode iopt(1)=1 is used,the values of tv(1)
149
c ...,tv(nv) should be left unchanged between subsequent calls.
150
c if the computation mode iopt(1)=-1 is used, the values tv(5),
151
c ...tv(nv-4) must be supplied by the user, before entry.
152
c see also the restrictions (ier=10).
153
c c : real array of dimension at least (nuest-4)*(nvest-4).
154
c on succesful exit, c contains the coefficients of the spline
155
c approximation s(u,v)
156
c fp : real. unless ier=10, fp contains the sum of squared
157
c residuals of the spline approximation returned.
158
c wrk : real array of dimension (lwrk). used as workspace.
159
c if the computation mode iopt(1)=1 is used the values of
160
c wrk(1),...,wrk(8) should be left unchanged between subsequent
162
c lwrk : integer. on entry lwrk must specify the actual dimension of
163
c the array wrk as declared in the calling (sub)program.
164
c lwrk must not be too small.
165
c lwrk >= 8+nuest*(mv+nvest+3)+nvest*21+4*mu+6*mv+q
166
c where q is the larger of (mv+nvest) and nuest.
167
c iwrk : integer array of dimension (kwrk). used as workspace.
168
c if the computation mode iopt(1)=1 is used the values of
169
c iwrk(1),.,iwrk(4) should be left unchanged between subsequent
171
c kwrk : integer. on entry kwrk must specify the actual dimension of
172
c the array iwrk as declared in the calling (sub)program.
173
c kwrk >= 4+mu+mv+nuest+nvest.
174
c ier : integer. unless the routine detects an error, ier contains a
175
c non-positive value on exit, i.e.
176
c ier=0 : normal return. the spline returned has a residual sum of
177
c squares fp such that abs(fp-s)/s <= tol with tol a relat-
178
c ive tolerance set to 0.001 by the program.
179
c ier=-1 : normal return. the spline returned is an interpolating
181
c ier=-2 : normal return. the spline returned is the least-squares
182
c constrained polynomial. in this extreme case fp gives the
183
c upper bound for the smoothing factor s.
184
c ier=1 : error. the required storage space exceeds the available
185
c storage space, as specified by the parameters nuest and
187
c probably causes : nuest or nvest too small. if these param-
188
c eters are already large, it may also indicate that s is
190
c the approximation returned is the least-squares spline
191
c according to the current set of knots. the parameter fp
192
c gives the corresponding sum of squared residuals (fp>s).
193
c ier=2 : error. a theoretically impossible result was found during
194
c the iteration proces for finding a smoothing spline with
195
c fp = s. probably causes : s too small.
196
c there is an approximation returned but the corresponding
197
c sum of squared residuals does not satisfy the condition
199
c ier=3 : error. the maximal number of iterations maxit (set to 20
200
c by the program) allowed for finding a smoothing spline
201
c with fp=s has been reached. probably causes : s too small
202
c there is an approximation returned but the corresponding
203
c sum of squared residuals does not satisfy the condition
205
c ier=10 : error. on entry, the input data are controlled on validity
206
c the following restrictions must be satisfied.
207
c -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1,
208
c -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0.
209
c mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8,
210
c kwrk>=4+mu+mv+nuest+nvest,
211
c lwrk >= 8+nuest*(mv+nvest+3)+nvest*21+4*mu+6*mv+
212
c max(nuest,mv+nvest)
213
c 0< u(i-1)<u(i)<=r,i=2,..,mu, (< r if iopt(3)=1)
214
c -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv
215
c if iopt(1)=-1: 8<=nu<=min(nuest,mu+5+iopt(2)+iopt(3))
216
c 0<tu(5)<tu(6)<...<tu(nu-4)<r
217
c 8<=nv<=min(nvest,mv+7)
218
c v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi
219
c the schoenberg-whitney conditions, i.e. there must
220
c be subset of grid co-ordinates uu(p) and vv(q) such
221
c that tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
222
c (iopt(2)=1 and iopt(3)=1 also count for a uu-value
223
c tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
224
c (vv(q) is either a value v(j) or v(j)+2*pi)
225
c if iopt(1)>=0: s>=0
226
c if s=0: nuest>=mu+5+iopt(2)+iopt(3), nvest>=mv+7
227
c if one of these conditions is found to be violated,control
228
c is immediately repassed to the calling program. in that
229
c case there is no approximation returned.
232
c pogrid does not allow individual weighting of the data-values.
233
c so, if these were determined to widely different accuracies, then
234
c perhaps the general data set routine polar should rather be used
235
c in spite of efficiency.
236
c by means of the parameter s, the user can control the tradeoff
237
c between closeness of fit and smoothness of fit of the approximation.
238
c if s is too large, the spline will be too smooth and signal will be
239
c lost ; if s is too small the spline will pick up too much noise. in
240
c the extreme cases the program will return an interpolating spline if
241
c s=0 and the constrained least-squares polynomial(degrees 3,0)if s is
242
c very large. between these extremes, a properly chosen s will result
243
c in a good compromise between closeness of fit and smoothness of fit.
244
c to decide whether an approximation, corresponding to a certain s is
245
c satisfactory the user is highly recommended to inspect the fits
247
c recommended values for s depend on the accuracy of the data values.
248
c if the user has an idea of the statistical errors on the data, he
249
c can also find a proper estimate for s. for, by assuming that, if he
250
c specifies the right s, pogrid will return a spline s(u,v) which
251
c exactly reproduces the function underlying the data he can evaluate
252
c the sum((z(i,j)-s(u(i),v(j)))**2) to find a good estimate for this s
253
c for example, if he knows that the statistical errors on his z(i,j)-
254
c values is not greater than 0.1, he may expect that a good s should
255
c have a value not larger than mu*mv*(0.1)**2.
256
c if nothing is known about the statistical error in z(i,j), s must
257
c be determined by trial and error, taking account of the comments
258
c above. the best is then to start with a very large value of s (to
259
c determine the least-squares polynomial and the corresponding upper
260
c bound fp0 for s) and then to progressively decrease the value of s
261
c ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
262
c and more carefully as the approximation shows more detail) to
263
c obtain closer fits.
264
c to economize the search for a good s-value the program provides with
265
c different modes of computation. at the first call of the routine, or
266
c whenever he wants to restart with the initial set of knots the user
267
c must set iopt(1)=0.
268
c if iopt(1) = 1 the program will continue with the knots found at
269
c the last call of the routine. this will save a lot of computation
270
c time if pogrid is called repeatedly for different values of s.
271
c the number of knots of the spline returned and their location will
272
c depend on the value of s and on the complexity of the shape of the
273
c function underlying the data. if the computation mode iopt(1) = 1
274
c is used, the knots returned may also depend on the s-values at
275
c previous calls (if these were smaller). therefore, if after a number
276
c of trials with different s-values and iopt(1)=1,the user can finally
277
c accept a fit as satisfactory, it may be worthwhile for him to call
278
c pogrid once more with the chosen value for s but now with iopt(1)=0.
279
c indeed, pogrid may then return an approximation of the same quality
280
c of fit but with fewer knots and therefore better if data reduction
281
c is also an important objective for the user.
282
c the number of knots may also depend on the upper bounds nuest and
283
c nvest. indeed, if at a certain stage in pogrid the number of knots
284
c in one direction (say nu) has reached the value of its upper bound
285
c (nuest), then from that moment on all subsequent knots are added
286
c in the other (v) direction. this may indicate that the value of
287
c nuest is too small. on the other hand, it gives the user the option
288
c of limiting the number of knots the routine locates in any direction
289
c for example, by setting nuest=8 (the lowest allowable value for
290
c nuest), the user can indicate that he wants an approximation which
291
c is a simple cubic polynomial in the variable u.
293
c other subroutines required:
294
c fppogr,fpchec,fpchep,fpknot,fpopdi,fprati,fpgrdi,fpsysy,fpback,
295
c fpbacp,fpbspl,fpcyt1,fpcyt2,fpdisc,fpgivs,fprota
298
c dierckx p. : fast algorithms for smoothing data over a disc or a
299
c sphere using tensor product splines, in "algorithms
300
c for approximation", ed. j.c.mason and m.g.cox,
301
c clarendon press oxford, 1987, pp. 51-65
302
c dierckx p. : fast algorithms for smoothing data over a disc or a
303
c sphere using tensor product splines, report tw73, dept.
304
c computer science,k.u.leuven, 1985.
305
c dierckx p. : curve and surface fitting with splines, monographs on
306
c numerical analysis, oxford university press, 1993.
310
c dept. computer science, k.u. leuven
311
c celestijnenlaan 200a, b-3001 heverlee, belgium.
312
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
314
c creation date : july 1985
315
c latest update : march 1989
318
c ..scalar arguments..
320
integer mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier
321
c ..array arguments..
322
integer iopt(3),ider(2),iwrk(kwrk)
323
real*8 u(mu),v(mv),z(mu*mv),c((nuest-4)*(nvest-4)),tu(nuest),
324
* tv(nvest),wrk(lwrk)
326
real*8 per,pi,tol,uu,ve,zmax,zmin,one,half,rn,zb
327
integer i,i1,i2,j,jwrk,j1,j2,kndu,kndv,knru,knrv,kwest,l,
328
* ldz,lfpu,lfpv,lwest,lww,m,maxit,mumin,muu,nc
329
c ..function references..
332
c ..subroutine references..
333
c fpchec,fpchep,fppogr
338
pi = datan2(0d0,-one)
341
c we set up the parameters tol and maxit.
344
c before starting computations, a data check is made. if the input data
345
c are invalid, control is immediately repassed to the calling program.
347
if(iopt(1).lt.(-1) .or. iopt(1).gt.1) go to 200
348
if(iopt(2).lt.0 .or. iopt(2).gt.1) go to 200
349
if(iopt(3).lt.0 .or. iopt(3).gt.1) go to 200
350
if(ider(1).lt.(-1) .or. ider(1).gt.1) go to 200
351
if(ider(2).lt.0 .or. ider(2).gt.1) go to 200
352
if(ider(2).eq.1 .and. iopt(2).eq.0) go to 200
353
mumin = 4-iopt(3)-ider(2)
354
if(ider(1).ge.0) mumin = mumin-1
355
if(mu.lt.mumin .or. mv.lt.4) go to 200
356
if(nuest.lt.8 .or. nvest.lt.8) go to 200
358
nc = (nuest-4)*(nvest-4)
359
lwest = 8+nuest*(mv+nvest+3)+21*nvest+4*mu+6*mv+
360
* max0(nuest,mv+nvest)
361
kwest = 4+mu+mv+nuest+nvest
362
if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200
363
if(u(1).le.0. .or. u(mu).gt.r) go to 200
364
if(iopt(3).eq.0) go to 10
365
if(u(mu).eq.r) go to 200
366
10 if(mu.eq.1) go to 30
368
if(u(i-1).ge.u(i)) go to 200
370
30 if(v(1).lt. (-pi) .or. v(1).ge.pi ) go to 200
371
if(v(mv).ge.v(1)+per) go to 200
373
if(v(i-1).ge.v(i)) go to 200
375
if(iopt(1).gt.0) go to 140
376
c if not given, we compute an estimate for z0.
377
if(ider(1).lt.0) go to 50
386
c we determine the range of z-values.
390
if(z(i).lt.zmin) zmin = z(i)
391
if(z(i).gt.zmax) zmax = z(i)
398
if(iopt(1).eq.0) go to 140
399
if(nu.lt.8 .or. nu.gt.nuest) go to 200
400
if(nv.lt.11 .or. nv.gt.nvest) go to 200
409
if(iopt(2).eq.0) go to 100
412
if(uu.gt.tu(5)) uu = tu(5)
418
if(iopt(3).eq.0) go to 120
422
call fpchec(wrk(9),muu,tu,nu,3,ier)
423
if(ier.ne.0) go to 200
444
call fpchep(wrk(9),mv+1,tv,nv,3,ier)
445
if (ier.eq.0) go to 150
447
140 if(s.lt.0.) go to 200
448
if(s.eq.0. .and. (nuest.lt.(mu+5+iopt(2)+iopt(3)) .or.
449
* nvest.lt.(mv+7)) ) go to 200
450
c we partition the working space and determine the spline approximation
455
jwrk = lwrk-8-nuest-nvest
460
call fppogr(iopt,ider,u,mu,v,mv,z,m,zb,r,s,nuest,nvest,tol,maxit,
461
* nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),wrk(lfpu),
462
* wrk(lfpv),wrk(ldz),wrk(8),iwrk(1),iwrk(2),iwrk(3),iwrk(4),
463
* iwrk(knru),iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),jwrk,ier)
added
removed
Lib/sandbox/spline/fitpack/pogrid.f
