real function smooth ( x, y, dy, npoint, s, v, a ) c from * a practical guide to splines * by c. de boor calls setupq, chol1d c c constructs the cubic smoothing spline f to given data (x(i),y(i)), c i=1,...,npoint, which has as small a second derivative as possible c while c s(f) = sum( ((y(i)-f(x(i)))/dy(i))**2 , i=1,...,npoint ) .le. s . c c****** i n p u t ****** c x(1),...,x(npoint) data abscissae, a s s u m e d to be strictly c increasing . c y(1),...,y(npoint) corresponding data ordinates . c dy(1),...,dy(npoint) estimate of uncertainty in data, a s s u m- c e d to be positive . c npoint.....number of data points, a s s u m e d .gt. 1 c s.....upper bound on the discrete weighted mean square distance of c the approximation f from the data . c c****** w o r k a r r a y s ***** c v.....of size (npoint,7) c a.....of size (npoint,4) c c***** o u t p u t ***** c a(.,1).....contains the sequence of smoothed ordinates . c a(i,j) = f^(j-1)(x(i)), j=2,3,4, i=1,...,npoint-1 , i.e., the c first three derivatives of the smoothing spline f at the c left end of each of the data intervals . c w a r n i n g . . . a would have to be transposed before it c could be used in ppvalu . c c****** m e t h o d ****** c The matrices Q-transp*d and Q-transp*D**2*Q are constructed in c s e t u p q from x and dy , as is the vector qty = Q-transp*y . c Then, for given p , the vector u is determined in c h o l 1 d as c the solution of the linear system c (6(1-p)Q-transp*D**2*Q + p*Q)u = qty . c From u , the smoothing spline f (for this choice of smoothing par- c ameter p ) is obtained in the sense that c f(x(.)) = y - 6(1-p)D**2*Q*u and c f''(x(.)) = 6*p*u . c The smoothing parameter p is found (if possible) so that c sf(p) = s , c with sf(p) = s(f) , where f is the smoothing spline as it depends c on p . if s = 0, then p = 1 . if sf(0) .le. s , then p = 0 . c Otherwise, the secant method is used to locate an appropriate p in c the open interval (0,1) . However, straightforward application of c the secant method, as done in the original version of this program, c can be very slow and is influenced by the units in which x and y c are measured, as C. Reinsch has pointed out. Instead, on recommend- c ation from C. Reinsch, the secant method is applied to the function c g:q |--> 1/sqrt{sfq(q)} - 1/sqrt{s} , c with sfq(q) := sf(q/(1+q)), since 1/sqrt{sfq} is monotone increasing c and close to linear for larger q . One starts at q = 0 with a c Newton step, i.e., c q_0 = 0, q_1 = -g(0)/g'(0) c with g'(0) = -(1/2) sfq(0)^{-3/2} dsfq, where dsfq = -12*u-transp*r*u , c and u as obtained for p = 0 . Iteration terminates as soon as c abs(sf - s) .le. .01*s . c c logical test c parameter (test = .true.) c integer itercnt integer npoint, i,npm1 real a(npoint,4),dy(npoint),s,v(npoint,7),x(npoint),y(npoint) * ,change,ooss,oosf,p,prevsf,prevq,q,sfq,sixp,six1mp,utru call setupq(x,dy,y,npoint,v,a(1,4)) if (s .gt. 0.) go to 20 10 p = 1. call chol1d(p,v,a(1,4),npoint,1,a(1,3),a(1,1)) sfq = 0. go to 60 20 p = 0. call chol1d(p,v,a(1,4),npoint,1,a(1,3),a(1,1)) sfq = 0. do 21 i=1,npoint 21 sfq = sfq + (a(i,1)*dy(i))**2 sfq = sfq*36. if (sfq .le. s) go to 60 utru = 0. do 25 i=2,npoint 25 utru = utru + v(i-1,4)*(a(i-1,3)*(a(i-1,3)+a(i,3))+a(i,3)**2) ooss = 1./sqrt(s) oosf = 1./sqrt(sfq) q = -(oosf-ooss)*sfq/(6.*utru*oosf) c secant iteration for the determination of p starts here. c itercnt = 0 prevq = 0. prevsf = oosf 30 call chol1d(q/(1.+q),v,a(1,4),npoint,1,a(1,3),a(1,1)) sfq = 0. do 35 i=1,npoint 35 sfq = sfq + (a(i,1)*dy(i))**2 sfq = sfq*36./(1.+q)**2 if (abs(sfq-s) .le. .01*s) go to 59 oosf = 1./sqrt(sfq) change = (q-prevq)/(oosf-prevsf)*(oosf-ooss) prevq = q q = q - change prevsf = oosf c itercnt = itercnt + 1 go to 30 59 p = q/(1.+q) correct value of p has been found. compute pol.coefficients from Q*u (in a(.,1)). 60 smooth = sfq c if (test) then c print *, 'number of iterations = ', itercnt c end if six1mp = 6./(1.+q) do 61 i=1,npoint 61 a(i,1) = y(i) - six1mp*dy(i)**2*a(i,1) sixp = 6.*p do 62 i=1,npoint 62 a(i,3) = a(i,3)*sixp npm1 = npoint - 1 do 63 i=1,npm1 a(i,4) = (a(i+1,3)-a(i,3))/v(i,4) 63 a(i,2) = (a(i+1,1)-a(i,1))/v(i,4) * - (a(i,3)+a(i,4)/3.*v(i,4))/2.*v(i,4) return end