1      subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,
2     * n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
3c  ..
4c  ..scalar arguments..
5      real*8 xb,xe,s,tol,fp
6      integer iopt,m,k,nest,maxit,k1,k2,n,ier
7c  ..array arguments..
8      real*8 x(m),y(m),w(m),t(nest),c(nest),fpint(nest),
9     * z(nest),a(nest,k1),b(nest,k2),g(nest,k2),q(m,k1)
10      integer nrdata(nest)
11c  ..local scalars..
12      real*8 acc,con1,con4,con9,cos,half,fpart,fpms,fpold,fp0,f1,f2,f3,
13     * one,p,pinv,piv,p1,p2,p3,rn,sin,store,term,wi,xi,yi
14      integer i,ich1,ich3,it,iter,i1,i2,i3,j,k3,l,l0,
15     * mk1,new,nk1,nmax,nmin,nplus,npl1,nrint,n8
16c  ..local arrays..
17      real*8 h(7)
18c  ..function references
19      real*8 abs,fprati
20      integer max0,min0
21c  ..subroutine references..
22c    fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
23c  ..
24c  set constants
25      one = 0.1d+01
26      con1 = 0.1d0
27      con9 = 0.9d0
28      con4 = 0.4d-01
29      half = 0.5d0
30cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
31c  part 1: determination of the number of knots and their position     c
32c  **************************************************************      c
33c  given a set of knots we compute the least-squares spline sinf(x),   c
34c  and the corresponding sum of squared residuals fp=f(p=inf).         c
35c  if iopt=-1 sinf(x) is the requested approximation.                  c
36c  if iopt=0 or iopt=1 we check whether we can accept the knots:       c
37c    if fp <=s we will continue with the current set of knots.         c
38c    if fp > s we will increase the number of knots and compute the    c
39c       corresponding least-squares spline until finally fp<=s.        c
40c    the initial choice of knots depends on the value of s and iopt.   c
41c    if s=0 we have spline interpolation; in that case the number of   c
42c    knots equals nmax = m+k+1.                                        c
43c    if s > 0 and                                                      c
44c      iopt=0 we first compute the least-squares polynomial of         c
45c      degree k; n = nmin = 2*k+2                                      c
46c      iopt=1 we start with the set of knots found at the last         c
47c      call of the routine, except for the case that s > fp0; then     c
48c      we compute directly the least-squares polynomial of degree k.   c
49cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
50c  determine nmin, the number of knots for polynomial approximation.
51      nmin = 2*k1
52      if(iopt.lt.0) go to 60
53c  calculation of acc, the absolute tolerance for the root of f(p)=s.
54      acc = tol*s
55c  determine nmax, the number of knots for spline interpolation.
56      nmax = m+k1
57      if(s.gt.0.0d0) go to 45
58c  if s=0, s(x) is an interpolating spline.
59c  test whether the required storage space exceeds the available one.
60      n = nmax
61      if(nmax.gt.nest) go to 420
62c  find the position of the interior knots in case of interpolation.
63  10  mk1 = m-k1
64      if(mk1.eq.0) go to 60
65      k3 = k/2
66      i = k2
67      j = k3+2
68      if(k3*2.eq.k) go to 30
69      do 20 l=1,mk1
70        t(i) = x(j)
71        i = i+1
72        j = j+1
73  20  continue
74      go to 60
75  30  do 40 l=1,mk1
76        t(i) = (x(j)+x(j-1))*half
77        i = i+1
78        j = j+1
79  40  continue
80      go to 60
81c  if s>0 our initial choice of knots depends on the value of iopt.
82c  if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
83c  polynomial of degree k which is a spline without interior knots.
84c  if iopt=1 and fp0>s we start computing the least squares spline
85c  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
88      fp0 = fpint(n)
89      fpold = fpint(n-1)
90      nplus = nrdata(n)
91      if(fp0.gt.s) go to 60
92  50  n = nmin
93      fpold = 0.0d0
94      nplus = 0
95      nrdata(1) = m-2
96c  main loop for the different sets of knots. m is a save upper bound
97c  for the number of trials.
98  60  do 200 iter = 1,m
99        if(n.eq.nmin) ier = -2
100c  find nrint, tne number of knot intervals.
101        nrint = n-nmin+1
102c  find the position of the additional knots which are needed for
103c  the b-spline representation of s(x).
104        nk1 = n-k1
105        i = n
106        do 70 j=1,k1
107          t(j) = xb
108          t(i) = xe
109          i = i-1
110  70    continue
111c  compute the b-spline coefficients of the least-squares spline
112c  sinf(x). the observation matrix a is built up row by row and
113c  reduced to upper triangular form by givens transformations.
114c  at the same time fp=f(p=inf) is computed.
115        fp = 0.0d0
116c  initialize the observation matrix a.
117        do 80 i=1,nk1
118          z(i) = 0.0d0
119          do 80 j=1,k1
120            a(i,j) = 0.0d0
121  80    continue
122        l = k1
123        do 130 it=1,m
124c  fetch the current data point x(it),y(it).
125          xi = x(it)
126          wi = w(it)
127          yi = y(it)*wi
128c  search for knot interval t(l) <= xi < t(l+1).
129  85      if(xi.lt.t(l+1) .or. l.eq.nk1) go to 90
130          l = l+1
131          go to 85
132c  evaluate the (k+1) non-zero b-splines at xi and store them in q.
133  90      call fpbspl(t,n,k,xi,l,h)
134          do 95 i=1,k1
135            q(it,i) = h(i)
136            h(i) = h(i)*wi
137  95      continue
138c  rotate the new row of the observation matrix into triangle.
139          j = l-k1
140          do 110 i=1,k1
141            j = j+1
142            piv = h(i)
143            if(piv.eq.0.0d0) go to 110
144c  calculate the parameters of the givens transformation.
145            call fpgivs(piv,a(j,1),cos,sin)
146c  transformations to right hand side.
147            call fprota(cos,sin,yi,z(j))
148            if(i.eq.k1) go to 120
149            i2 = 1
150            i3 = i+1
151            do 100 i1 = i3,k1
152              i2 = i2+1
153c  transformations to left hand side.
154              call fprota(cos,sin,h(i1),a(j,i2))
155 100        continue
156 110      continue
157c  add contribution of this row to the sum of squares of residual
158c  right hand sides.
159 120      fp = fp+yi*yi
160 130    continue
161        if(ier.eq.(-2)) fp0 = fp
162        fpint(n) = fp0
163        fpint(n-1) = fpold
164        nrdata(n) = nplus
165c  backward substitution to obtain the b-spline coefficients.
166        call fpback(a,z,nk1,k1,c,nest)
167c  test whether the approximation sinf(x) is an acceptable solution.
168        if(iopt.lt.0) go to 440
169        fpms = fp-s
170        if(abs(fpms).lt.acc) go to 440
171c  if f(p=inf) < s accept the choice of knots.
172        if(fpms.lt.0.0d0) go to 250
173c  if n = nmax, sinf(x) is an interpolating spline.
174        if(n.eq.nmax) go to 430
175c  increase the number of knots.
176c  if n=nest we cannot increase the number of knots because of
177c  the storage capacity limitation.
178        if(n.eq.nest) go to 420
179c  determine the number of knots nplus we are going to add.
180        if(ier.eq.0) go to 140
181        nplus = 1
182        ier = 0
183        go to 150
184 140    npl1 = nplus*2
185        rn = nplus
186        if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
187        nplus = min0(nplus*2,max0(npl1,nplus/2,1))
188 150    fpold = fp
189c  compute the sum((w(i)*(y(i)-s(x(i))))**2) for each knot interval
190c  t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
191        fpart = 0.0d0
192        i = 1
193        l = k2
194        new = 0
195        do 180 it=1,m
196          if(x(it).lt.t(l) .or. l.gt.nk1) go to 160
197          new = 1
198          l = l+1
199 160      term = 0.0d0
200          l0 = l-k2
201          do 170 j=1,k1
202            l0 = l0+1
203            term = term+c(l0)*q(it,j)
204 170      continue
205          term = (w(it)*(term-y(it)))**2
206          fpart = fpart+term
207          if(new.eq.0) go to 180
208          store = term*half
209          fpint(i) = fpart-store
210          i = i+1
211          fpart = store
212          new = 0
213 180    continue
214        fpint(nrint) = fpart
215        do 190 l=1,nplus
216c  add a new knot.
217          call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
218c  if n=nmax we locate the knots as for interpolation.
219          if(n.eq.nmax) go to 10
220c  test whether we cannot further increase the number of knots.
221          if(n.eq.nest) go to 200
222 190    continue
223c  restart the computations with the new set of knots.
224 200  continue
225c  test whether the least-squares kth degree polynomial is a solution
226c  of our approximation problem.
227 250  if(ier.eq.(-2)) go to 440
228cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
229c  part 2: determination of the smoothing spline sp(x).                c
230c  ***************************************************                 c
231c  we have determined the number of knots and their position.          c
232c  we now compute the b-spline coefficients of the smoothing spline    c
233c  sp(x). the observation matrix a is extended by the rows of matrix   c
234c  b expressing that the kth derivative discontinuities of sp(x) at    c
235c  the interior knots t(k+2),...t(n-k-1) must be zero. the corres-     c
236c  ponding weights of these additional rows are set to 1/p.            c
237c  iteratively we then have to determine the value of p such that      c
238c  f(p)=sum((w(i)*(y(i)-sp(x(i))))**2) be = s. we already know that    c
239c  the least-squares kth degree polynomial corresponds to p=0, and     c
240c  that the least-squares spline corresponds to p=infinity. the        c
241c  iteration process which is proposed here, makes use of rational     c
242c  interpolation. since f(p) is a convex and strictly decreasing       c
243c  function of p, it can be approximated by a rational function        c
244c  r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond-  c
245c  ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used      c
246c  to calculate the new value of p such that r(p)=s. convergence is    c
247c  guaranteed by taking f1>0 and f3<0.                                 c
248cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
249c  evaluate the discontinuity jump of the kth derivative of the
250c  b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
251      call fpdisc(t,n,k2,b,nest)
252c  initial value for p.
253      p1 = 0.0d0
254      f1 = fp0-s
255      p3 = -one
256      f3 = fpms
257      p = 0.
258      do 255 i=1,nk1
259         p = p+a(i,1)
260 255  continue
261      rn = nk1
262      p = rn/p
263      ich1 = 0
264      ich3 = 0
265      n8 = n-nmin
266c  iteration process to find the root of f(p) = s.
267      do 360 iter=1,maxit
268c  the rows of matrix b with weight 1/p are rotated into the
269c  triangularised observation matrix a which is stored in g.
270        pinv = one/p
271        do 260 i=1,nk1
272          c(i) = z(i)
273          g(i,k2) = 0.0d0
274          do 260 j=1,k1
275            g(i,j) = a(i,j)
276 260    continue
277        do 300 it=1,n8
278c  the row of matrix b is rotated into triangle by givens transformation
279          do 270 i=1,k2
280            h(i) = b(it,i)*pinv
281 270      continue
282          yi = 0.0d0
283          do 290 j=it,nk1
284            piv = h(1)
285c  calculate the parameters of the givens transformation.
286            call fpgivs(piv,g(j,1),cos,sin)
287c  transformations to right hand side.
288            call fprota(cos,sin,yi,c(j))
289            if(j.eq.nk1) go to 300
290            i2 = k1
291            if(j.gt.n8) i2 = nk1-j
292            do 280 i=1,i2
293c  transformations to left hand side.
294              i1 = i+1
295              call fprota(cos,sin,h(i1),g(j,i1))
296              h(i) = h(i1)
297 280        continue
298            h(i2+1) = 0.0d0
299 290      continue
300 300    continue
301c  backward substitution to obtain the b-spline coefficients.
302        call fpback(g,c,nk1,k2,c,nest)
303c  computation of f(p).
304        fp = 0.0d0
305        l = k2
306        do 330 it=1,m
307          if(x(it).lt.t(l) .or. l.gt.nk1) go to 310
308          l = l+1
309 310      l0 = l-k2
310          term = 0.0d0
311          do 320 j=1,k1
312            l0 = l0+1
313            term = term+c(l0)*q(it,j)
314 320      continue
315          fp = fp+(w(it)*(term-y(it)))**2
316 330    continue
317c  test whether the approximation sp(x) is an acceptable solution.
318        fpms = fp-s
319        if(abs(fpms).lt.acc) go to 440
320c  test whether the maximal number of iterations is reached.
321        if(iter.eq.maxit) go to 400
322c  carry out one more step of the iteration process.
323        p2 = p
324        f2 = fpms
325        if(ich3.ne.0) go to 340
326        if((f2-f3).gt.acc) go to 335
327c  our initial choice of p is too large.
328        p3 = p2
329        f3 = f2
330        p = p*con4
331        if(p.le.p1) p=p1*con9 + p2*con1
332        go to 360
333 335    if(f2.lt.0.0d0) ich3=1
334 340    if(ich1.ne.0) go to 350
335        if((f1-f2).gt.acc) go to 345
336c  our initial choice of p is too small
337        p1 = p2
338        f1 = f2
339        p = p/con4
340        if(p3.lt.0.) go to 360
341        if(p.ge.p3) p = p2*con1 + p3*con9
342        go to 360
343 345    if(f2.gt.0.0d0) ich1=1
344c  test whether the iteration process proceeds as theoretically
345c  expected.
346 350    if(f2.ge.f1 .or. f2.le.f3) go to 410
347c  find the new value for p.
348        p = fprati(p1,f1,p2,f2,p3,f3)
349 360  continue
350c  error codes and messages.
351 400  ier = 3
352      go to 440
353 410  ier = 2
354      go to 440
355 420  ier = 1
356      go to 440
357 430  ier = -1
358 440  return
359      end
360