1      subroutine parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest,
2     * nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
3c  given the set of ordered points f(i,j) in the idim-dimensional space,
4c  corresponding to grid values (u(i),v(j)) ,i=1,...,mu ; j=1,...,mv,
5c  parsur determines a smooth approximating spline surface s(u,v) , i.e.
6c    f1 = s1(u,v)
7c      ...                u(1) <= u <= u(mu) ; v(1) <= v <= v(mv)
8c    fidim = sidim(u,v)
9c  with sl(u,v), l=1,2,...,idim bicubic spline functions with common
10c  knots tu(i),i=1,...,nu in the u-variable and tv(j),j=1,...,nv in the
11c  v-variable.
12c  in addition, these splines will be periodic in the variable u if
13c  ipar(1) = 1 and periodic in the variable v if ipar(2) = 1.
14c  if iopt=-1, parsur determines the least-squares bicubic spline
15c  surface according to a given set of knots.
16c  if iopt>=0, the number of knots of s(u,v) and their position
17c  is chosen automatically by the routine. the smoothness of s(u,v) is
18c  achieved by minimalizing the discontinuity jumps of the derivatives
19c  of the splines at the knots. the amount of smoothness of s(u,v) is
20c  determined by the condition that
21c  fp=sumi=1,mu(sumj=1,mv(dist(f(i,j)-s(u(i),v(j)))**2))<=s,
22c  with s a given non-negative constant.
23c  the fit s(u,v) is given in its b-spline representation and can be
24c  evaluated by means of routine surev.
25c
26c calling sequence:
27c     call parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest,nu,tu,
28c    *  nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
29c
30c parameters:
31c  iopt  : integer flag. unchanged on exit.
32c          on entry iopt must specify whether a least-squares surface
33c          (iopt=-1) or a smoothing surface (iopt=0 or 1)must be
34c          determined.
35c          if iopt=0 the routine will start with the initial set of
36c          knots needed for determining the least-squares polynomial
37c          surface.
38c          if iopt=1 the routine will continue with the set of knots
39c          found at the last call of the routine.
40c          attention: a call with iopt=1 must always be immediately
41c          preceded by another call with iopt = 1 or iopt = 0.
42c  ipar  : integer array of dimension 2. unchanged on exit.
43c          on entry ipar(1) must specify whether (ipar(1)=1) or not
44c          (ipar(1)=0) the splines must be periodic in the variable u.
45c          on entry ipar(2) must specify whether (ipar(2)=1) or not
46c          (ipar(2)=0) the splines must be periodic in the variable v.
47c  idim  : integer. on entry idim must specify the dimension of the
48c          surface. 1 <= idim <= 3. unchanged on exit.
49c  mu    : integer. on entry mu must specify the number of grid points
50c          along the u-axis. unchanged on exit.
51c          mu >= mumin where mumin=4-2*ipar(1)
52c  u     : real array of dimension at least (mu). before entry, u(i)
53c          must be set to the u-co-ordinate of the i-th grid point
54c          along the u-axis, for i=1,2,...,mu. these values must be
55c          supplied in strictly ascending order. unchanged on exit.
56c  mv    : integer. on entry mv must specify the number of grid points
57c          along the v-axis. unchanged on exit.
58c          mv >= mvmin where mvmin=4-2*ipar(2)
59c  v     : real array of dimension at least (mv). before entry, v(j)
60c          must be set to the v-co-ordinate of the j-th grid point
61c          along the v-axis, for j=1,2,...,mv. these values must be
62c          supplied in strictly ascending order. unchanged on exit.
63c  f     : real array of dimension at least (mu*mv*idim).
64c          before entry, f(mu*mv*(l-1)+mv*(i-1)+j) must be set to the
65c          l-th co-ordinate of the data point corresponding to the
66c          the grid point (u(i),v(j)) for l=1,...,idim ,i=1,...,mu
67c          and j=1,...,mv. unchanged on exit.
68c          if ipar(1)=1 it is expected that f(mu*mv*(l-1)+mv*(mu-1)+j)
69c          = f(mu*mv*(l-1)+j), l=1,...,idim ; j=1,...,mv
70c          if ipar(2)=1 it is expected that f(mu*mv*(l-1)+mv*(i-1)+mv)
71c          = f(mu*mv*(l-1)+mv*(i-1)+1), l=1,...,idim ; i=1,...,mu
72c  s     : real. on entry (if iopt>=0) s must specify the smoothing
73c          factor. s >=0. unchanged on exit.
74c          for advice on the choice of s see further comments
75c  nuest : integer. unchanged on exit.
76c  nvest : integer. unchanged on exit.
77c          on entry, nuest and nvest must specify an upper bound for the
78c          number of knots required in the u- and v-directions respect.
79c          these numbers will also determine the storage space needed by
80c          the routine. nuest >= 8, nvest >= 8.
81c          in most practical situation nuest = mu/2, nvest=mv/2, will
82c          be sufficient. always large enough are nuest=mu+4+2*ipar(1),
83c          nvest = mv+4+2*ipar(2), the number of knots needed for
84c          interpolation (s=0). see also further comments.
85c  nu    : integer.
86c          unless ier=10 (in case iopt>=0), nu will contain the total
87c          number of knots with respect to the u-variable, of the spline
88c          surface returned. if the computation mode iopt=1 is used,
89c          the value of nu should be left unchanged between subsequent
90c          calls. in case iopt=-1, the value of nu should be specified
91c          on entry.
92c  tu    : real array of dimension at least (nuest).
93c          on successful exit, this array will contain the knots of the
94c          splines with respect to the u-variable, i.e. the position of
95c          the interior knots tu(5),...,tu(nu-4) as well as the position
96c          of the additional knots tu(1),...,tu(4) and tu(nu-3),...,
97c          tu(nu) needed for the b-spline representation.
98c          if the computation mode iopt=1 is used,the values of tu(1)
99c          ...,tu(nu) should be left unchanged between subsequent calls.
100c          if the computation mode iopt=-1 is used, the values tu(5),
101c          ...tu(nu-4) must be supplied by the user, before entry.
102c          see also the restrictions (ier=10).
103c  nv    : integer.
104c          unless ier=10 (in case iopt>=0), nv will contain the total
105c          number of knots with respect to the v-variable, of the spline
106c          surface returned. if the computation mode iopt=1 is used,
107c          the value of nv should be left unchanged between subsequent
108c          calls. in case iopt=-1, the value of nv should be specified
109c          on entry.
110c  tv    : real array of dimension at least (nvest).
111c          on successful exit, this array will contain the knots of the
112c          splines with respect to the v-variable, i.e. the position of
113c          the interior knots tv(5),...,tv(nv-4) as well as the position
114c          of the additional knots tv(1),...,tv(4) and tv(nv-3),...,
115c          tv(nv) needed for the b-spline representation.
116c          if the computation mode iopt=1 is used,the values of tv(1)
117c          ...,tv(nv) should be left unchanged between subsequent calls.
118c          if the computation mode iopt=-1 is used, the values tv(5),
119c          ...tv(nv-4) must be supplied by the user, before entry.
120c          see also the restrictions (ier=10).
121c  c     : real array of dimension at least (nuest-4)*(nvest-4)*idim.
122c          on successful exit, c contains the coefficients of the spline
123c          approximation s(u,v)
124c  fp    : real. unless ier=10, fp contains the sum of squared
125c          residuals of the spline surface returned.
126c  wrk   : real array of dimension (lwrk). used as workspace.
127c          if the computation mode iopt=1 is used the values of
128c          wrk(1),...,wrk(4) should be left unchanged between subsequent
129c          calls.
130c  lwrk  : integer. on entry lwrk must specify the actual dimension of
131c          the array wrk as declared in the calling (sub)program.
132c          lwrk must not be too small.
133c           lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+
134c           4*(mu+mv)+q*idim where q is the larger of mv and nuest.
135c  iwrk  : integer array of dimension (kwrk). used as workspace.
136c          if the computation mode iopt=1 is used the values of
137c          iwrk(1),.,iwrk(3) should be left unchanged between subsequent
138c          calls.
139c  kwrk  : integer. on entry kwrk must specify the actual dimension of
140c          the array iwrk as declared in the calling (sub)program.
141c          kwrk >= 3+mu+mv+nuest+nvest.
142c  ier   : integer. unless the routine detects an error, ier contains a
143c          non-positive value on exit, i.e.
144c   ier=0  : normal return. the surface returned has a residual sum of
145c            squares fp such that abs(fp-s)/s <= tol with tol a relat-
146c            ive tolerance set to 0.001 by the program.
147c   ier=-1 : normal return. the spline surface returned is an
148c            interpolating surface (fp=0).
149c   ier=-2 : normal return. the surface returned is the least-squares
150c            polynomial surface. in this extreme case fp gives the
151c            upper bound for the smoothing factor s.
152c   ier=1  : error. the required storage space exceeds the available
153c            storage space, as specified by the parameters nuest and
154c            nvest.
155c            probably causes : nuest or nvest too small. if these param-
156c            eters are already large, it may also indicate that s is
157c            too small
158c            the approximation returned is the least-squares surface
159c            according to the current set of knots. the parameter fp
160c            gives the corresponding sum of squared residuals (fp>s).
161c   ier=2  : error. a theoretically impossible result was found during
162c            the iteration process for finding a smoothing surface with
163c            fp = s. probably causes : s too small.
164c            there is an approximation returned but the corresponding
165c            sum of squared residuals does not satisfy the condition
166c            abs(fp-s)/s < tol.
167c   ier=3  : error. the maximal number of iterations maxit (set to 20
168c            by the program) allowed for finding a smoothing surface
169c            with fp=s has been reached. probably causes : s too small
170c            there is an approximation returned but the corresponding
171c            sum of squared residuals does not satisfy the condition
172c            abs(fp-s)/s < tol.
173c   ier=10 : error. on entry, the input data are controlled on validity
174c            the following restrictions must be satisfied.
175c            -1<=iopt<=1, 0<=ipar(1)<=1, 0<=ipar(2)<=1, 1 <=idim<=3
176c            mu >= 4-2*ipar(1),mv >= 4-2*ipar(2), nuest >=8, nvest >= 8,
177c            kwrk>=3+mu+mv+nuest+nvest,
178c            lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))
179c             +4*(mu+mv)+max(nuest,mv)*idim
180c            u(i-1)<u(i),i=2,..,mu, v(i-1)<v(i),i=2,...,mv
181c            if iopt=-1: 8<=nu<=min(nuest,mu+4+2*ipar(1))
182c                        u(1)<tu(5)<tu(6)<...<tu(nu-4)<u(mu)
183c                        8<=nv<=min(nvest,mv+4+2*ipar(2))
184c                        v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(mv)
185c                    the schoenberg-whitney conditions, i.e. there must
186c                    be subset of grid co-ordinates uu(p) and vv(q) such
187c                    that   tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
188c                           tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
189c                     (see fpchec or fpchep)
190c            if iopt>=0: s>=0
191c                       if s=0: nuest>=mu+4+2*ipar(1)
192c                               nvest>=mv+4+2*ipar(2)
193c            if one of these conditions is found to be violated,control
194c            is immediately repassed to the calling program. in that
195c            case there is no approximation returned.
196c
197c further comments:
198c   by means of the parameter s, the user can control the tradeoff
199c   between closeness of fit and smoothness of fit of the approximation.
200c   if s is too large, the surface will be too smooth and signal will be
201c   lost ; if s is too small the surface will pick up too much noise. in
202c   the extreme cases the program will return an interpolating surface
203c   if s=0 and the constrained least-squares polynomial surface if s is
204c   very large. between these extremes, a properly chosen s will result
205c   in a good compromise between closeness of fit and smoothness of fit.
206c   to decide whether an approximation, corresponding to a certain s is
207c   satisfactory the user is highly recommended to inspect the fits
208c   graphically.
209c   recommended values for s depend on the accuracy of the data values.
210c   if the user has an idea of the statistical errors on the data, he
211c   can also find a proper estimate for s. for, by assuming that, if he
212c   specifies the right s, parsur will return a surface s(u,v) which
213c   exactly reproduces the surface underlying the data he can evaluate
214c   the sum(dist(f(i,j)-s(u(i),v(j)))**2) to find a good estimate for s.
215c   for example, if he knows that the statistical errors on his f(i,j)-
216c   values is not greater than 0.1, he may expect that a good s should
217c   have a value not larger than mu*mv*(0.1)**2.
218c   if nothing is known about the statistical error in f(i,j), s must
219c   be determined by trial and error, taking account of the comments
220c   above. the best is then to start with a very large value of s (to
221c   determine the le-sq polynomial surface and the corresponding upper
222c   bound fp0 for s) and then to progressively decrease the value of s
223c   ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
224c   and more carefully as the approximation shows more detail) to
225c   obtain closer fits.
226c   to economize the search for a good s-value the program provides with
227c   different modes of computation. at the first call of the routine, or
228c   whenever he wants to restart with the initial set of knots the user
229c   must set iopt=0.
230c   if iopt = 1 the program will continue with the knots found at
231c   the last call of the routine. this will save a lot of computation
232c   time if parsur is called repeatedly for different values of s.
233c   the number of knots of the surface returned and their location will
234c   depend on the value of s and on the complexity of the shape of the
235c   surface underlying the data. if the computation mode iopt = 1
236c   is used, the knots returned may also depend on the s-values at
237c   previous calls (if these were smaller). therefore, if after a number
238c   of trials with different s-values and iopt=1,the user can finally
239c   accept a fit as satisfactory, it may be worthwhile for him to call
240c   parsur once more with the chosen value for s but now with iopt=0.
241c   indeed, parsur may then return an approximation of the same quality
242c   of fit but with fewer knots and therefore better if data reduction
243c   is also an important objective for the user.
244c   the number of knots may also depend on the upper bounds nuest and
245c   nvest. indeed, if at a certain stage in parsur the number of knots
246c   in one direction (say nu) has reached the value of its upper bound
247c   (nuest), then from that moment on all subsequent knots are added
248c   in the other (v) direction. this may indicate that the value of
249c   nuest is too small. on the other hand, it gives the user the option
250c   of limiting the number of knots the routine locates in any direction
251c   for example, by setting nuest=8 (the lowest allowable value for
252c   nuest), the user can indicate that he wants an approximation with
253c   splines which are simple cubic polynomials in the variable u.
254c
255c  other subroutines required:
256c    fppasu,fpchec,fpchep,fpknot,fprati,fpgrpa,fptrnp,fpback,
257c    fpbacp,fpbspl,fptrpe,fpdisc,fpgivs,fprota
258c
259c  author:
260c    p.dierckx
261c    dept. computer science, k.u. leuven
262c    celestijnenlaan 200a, b-3001 heverlee, belgium.
263c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
264c
265c  latest update : march 1989
266c
267c  ..
268c  ..scalar arguments..
269      real*8 s,fp
270      integer iopt,idim,mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier
271c  ..array arguments..
272      real*8 u(mu),v(mv),f(mu*mv*idim),tu(nuest),tv(nvest),
273     * c((nuest-4)*(nvest-4)*idim),wrk(lwrk)
274      integer ipar(2),iwrk(kwrk)
275c  ..local scalars..
276      real*8 tol,ub,ue,vb,ve,peru,perv
277      integer i,j,jwrk,kndu,kndv,knru,knrv,kwest,l1,l2,l3,l4,
278     * lfpu,lfpv,lwest,lww,maxit,nc,mf,mumin,mvmin
279c  ..function references..
280      integer max0
281c  ..subroutine references..
282c    fppasu,fpchec,fpchep
283c  ..
284c  we set up the parameters tol and maxit.
285      maxit = 20
286      tol = 0.1e-02
287c  before starting computations a data check is made. if the input data
288c  are invalid, control is immediately repassed to the calling program.
289      ier = 10
290      if(iopt.lt.(-1) .or. iopt.gt.1) go to 200
291      if(ipar(1).lt.0 .or. ipar(1).gt.1) go to 200
292      if(ipar(2).lt.0 .or. ipar(2).gt.1) go to 200
293      if(idim.le.0 .or. idim.gt.3) go to 200
294      mumin = 4-2*ipar(1)
295      if(mu.lt.mumin .or. nuest.lt.8) go to 200
296      mvmin = 4-2*ipar(2)
297      if(mv.lt.mvmin .or. nvest.lt.8) go to 200
298      mf = mu*mv
299      nc = (nuest-4)*(nvest-4)
300      lwest = 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+
301     * 4*(mu+mv)+max0(nuest,mv)*idim
302      kwest = 3+mu+mv+nuest+nvest
303      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200
304      do 10 i=2,mu
305        if(u(i-1).ge.u(i)) go to 200
306  10  continue
307      do 20 i=2,mv
308        if(v(i-1).ge.v(i)) go to 200
309  20  continue
310      if(iopt.ge.0) go to 100
311      if(nu.lt.8 .or. nu.gt.nuest) go to 200
312      ub = u(1)
313      ue = u(mu)
314      if (ipar(1).ne.0) go to 40
315      j = nu
316      do 30 i=1,4
317        tu(i) = ub
318        tu(j) = ue
319        j = j-1
320  30  continue
321      call fpchec(u,mu,tu,nu,3,ier)
322      if(ier.ne.0) go to 200
323      go to 60
324  40  l1 = 4
325      l2 = l1
326      l3 = nu-3
327      l4 = l3
328      peru = ue-ub
329      tu(l2) = ub
330      tu(l3) = ue
331      do 50 j=1,3
332        l1 = l1+1
333        l2 = l2-1
334        l3 = l3+1
335        l4 = l4-1
336        tu(l2) = tu(l4)-peru
337        tu(l3) = tu(l1)+peru
338  50  continue
339      call fpchep(u,mu,tu,nu,3,ier)
340      if(ier.ne.0) go to 200
341  60  if(nv.lt.8 .or. nv.gt.nvest) go to 200
342      vb = v(1)
343      ve = v(mv)
344      if (ipar(2).ne.0) go to 80
345      j = nv
346      do 70 i=1,4
347        tv(i) = vb
348        tv(j) = ve
349        j = j-1
350  70  continue
351      call fpchec(v,mv,tv,nv,3,ier)
352      if(ier.ne.0) go to 200
353      go to 150
354  80  l1 = 4
355      l2 = l1
356      l3 = nv-3
357      l4 = l3
358      perv = ve-vb
359      tv(l2) = vb
360      tv(l3) = ve
361      do 90 j=1,3
362        l1 = l1+1
363        l2 = l2-1
364        l3 = l3+1
365        l4 = l4-1
366        tv(l2) = tv(l4)-perv
367        tv(l3) = tv(l1)+perv
368  90  continue
369      call fpchep(v,mv,tv,nv,3,ier)
370      if (ier.eq.0) go to 150
371      go to 200
372 100  if(s.lt.0.) go to 200
373      if(s.eq.0. .and. (nuest.lt.(mu+4+2*ipar(1)) .or.
374     * nvest.lt.(mv+4+2*ipar(2))) )go to 200
375      ier = 0
376c  we partition the working space and determine the spline approximation
377 150  lfpu = 5
378      lfpv = lfpu+nuest
379      lww = lfpv+nvest
380      jwrk = lwrk-4-nuest-nvest
381      knru = 4
382      knrv = knru+mu
383      kndu = knrv+mv
384      kndv = kndu+nuest
385      call fppasu(iopt,ipar,idim,u,mu,v,mv,f,mf,s,nuest,nvest,
386     * tol,maxit,nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),
387     * wrk(lfpu),wrk(lfpv),iwrk(1),iwrk(2),iwrk(3),iwrk(knru),
388     * iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),jwrk,ier)
389 200  return
390      end
391
392