chapter xiv, additional smoothing example, supplied by C. Reinsch
calls  bsplpp(bsplvb), ppvalu(interv), smooth(setupq,chol1d)
c     perturbed values of  sin , with the opportunity of scaling both
C     x  and  y  arbitrarily  (which does bring ruin to the version of
C     smooth  in the original pgs.
      integer i,j,npoint,nm1
      parameter (npoint=91, nm1 = npoint-1)
      real skalx, skaly, skalesy(6)
      data skalesy /1.e-12, 1., 100., 300., 400., 500. /
      real a(npoint,4),dely,dy(npoint),freq, rands(91)
     *    ,scrtch(npoint,7),sfp,sfpchk,x(npoint),xx,y(npoint)
C    the following random data were generated by matlab
      data rands /
     *-0.5621, -0.9059, 0.3577, 0.3586, 0.8694, -0.2330, 0.0388,
     *0.6619, -0.9309, -0.8931, 0.0594, 0.3423, -0.9846, -0.2332,
     *-0.8663, -0.1650, 0.3735, 0.1780, 0.8609, 0.6923, 0.0539,
     *-0.8161, 0.3078, -0.1680, 0.4024, 0.8206, 0.5244, -0.4751,
     *-0.9051, 0.4722, -0.3435, 0.2653, 0.5128, 0.9821, -0.2693,
     *-0.5059, 0.9651, 0.4453, 0.5067, 0.3030, -0.8546, 0.2633,
     *0.7694, -0.4546, -0.1272, 0.5330, -0.0445, -0.5245, -0.4502,
     *-0.2815, -0.6670, -0.0270, 0.7953, 0.8184, -0.8789, 0.8093,
     *0.0090, 0.0326, -0.3619, 0.9733, -0.0120, -0.4677, -0.8185,
     *0.8955, -0.8525, 0.0014, -0.2317, -0.4458, 0.8276, 0.0595,
     *-0.0711, 0.8820, -0.8998, 0.5230, 0.5404, 0.6556, -0.7493,
     *-0.9683, 0.3769, 0.7365, 0.2591, 0.4724, 0.4508, 0.9989,
     *0.7771, -0.5336, -0.3874, -0.2980, 0.0265, 0.1822, 0.6920/
      freq = 3.14159265/9.
      skalx = 1.
      do 100 j=1,6
      skaly = skalesy(j)
      dely = skaly*.05/sqrt(3.)
      do 10 i=1,npoint
         xx = float(i-1)*.1
         x(i) = skalx*xx
         y(i) = skaly*(sin(xx*freq) + .05*rands(i))
   10    dy(i) = dely
      sfp = smooth ( x, y, dy, npoint, 80., scrtch, a )
c  also compute  sfp from scratch here
      sfpchk = 0.
      do 90 i=1,npoint
         sfpchk = sfpchk + ((y(i) - a(i,1))/dy(i))**2
   90 continue
      print *, 'scales used: (',skalx,', ', skaly, '), sfp = ', sfp, sfpchk
  100 continue
c  also try interpolation, to be sure:
      sfp = smooth ( x, y, dy, npoint, 0., scrtch, a )
      sfpchk = 0.
c  the following check is trivial here since a(:,1) <-- y  in this case.
      do 110 i=1,npoint
         sfpchk = sfpchk + ((y(i) - a(i,1))/dy(i))**2
  110 continue
      print *, 'same scales as last, but S = 0, gives sfp = ', sfp, sfpchk
                                        stop
      end