chapter xiii, example 2. cubic spline interpolant at knot averages
c                        with good knots.
c  from  * a practical guide to splines *  by c. de boor    
calls splint(banfac/slv,bsplvb),newnot,bsplpp(bsplvb*)
      integer i,istep,iter,itermx,n,nhigh,nmax,nmk,nlow
      parameter (nmax=20)
      real algerp,aloger,bcoef(nmax+2),break(nmax),decay,dx,errmax
     *    ,c(4,nmax),g,gtau(nmax),h,pnatx,scrtch(nmax*7),step
     *    ,t(nmax+6),tau(nmax),tnew(nmax)
C     real algerp,aloger,bcoef(22),break(20),decay,dx,errmax,
C    *    c(4,20),g,gtau(20),h,pnatx,scrtch(140),step,t(26),
C    *    tau(20),tnew(20)
c     istep and step = float(istep) specify point density for error det-
c     ermination.
      data step, istep /20., 20/
c                               the function  g  is to be interpolated .
      g(x) = sqrt(x+1.)
      decay = 0.
c      read in the number of iterations to be carried out and the lower
c      and upper limit for the number  n  of data points to be used.
      read 500,itermx,nlow,nhigh
  500 format(3i3)
      print 600, itermx
  600 format(i4,22h cycles through newnot//
     *      28h  n   max.error   decay exp./)
c                                 loop over  n = number of data points .
      do 40 n=nlow,nhigh,2
    4    nmk = n-4
         h = 2./float(nmk+1)
         do 5 i=1,4
            t(i) = -1.
    5       t(n+i) = 1.
         if (nmk .lt. 1)                go to 10
         do 6 i=1,nmk
    6       t(i+4) = float(i)*h - 1.
   10    iter = 1
c               construct cubic spline interpolant. then, itermx  times,
c                determine new knots from it and find a new interpolant.
   11       do 12 i=1,n
               tau(i) = (t(i+1)+t(i+2)+t(i+3))/3.
   12          gtau(i) = g(tau(i))
            call splint ( tau, gtau, t, n, 4, scrtch, bcoef, iflag )
            call bsplpp ( t, bcoef, n, 4, scrtch, break, c, l )
            if (iter .gt. itermx)       go to 19
            iter = iter + 1
            call newnot ( break, c, l, 4, tnew, l, scrtch )
   15       do 18 i=2,l
   18          t(3+i) = tnew(i)
                                        go to 11
c                      estimate max.interpolation error on  (-1,1) .
   19    errmax = 0.
c                           loop over polynomial pieces of interpolant .
         do 30 i=1,l
            dx = (break(i+1)-break(i))/step
c                interpolation error is calculated at  istep  points per
c                    polynomial piece .
            do 30 j=1,istep
               h = float(j)*dx
               pnatx = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)
   30       errmax = amax1(errmax,abs(g(break(i)+h)-pnatx))
c                                             calculate decay exponent .
         aloger = alog(errmax)
         if (n .gt. nlow)  decay =
     *       (aloger - algerp)/alog(float(n)/float(n-2))
         algerp = aloger
c
   40    print 640,n,errmax,decay
  640 format(i3,e12.4,f11.2)
                                        stop
      end