Initial commit

21 files changed, 1071 insertions, 0 deletions

diff --git a/math/deboor/progs/prog1.f b/math/deboor/progs/prog1.f
new file mode 100644
index 00000000..2a0b869c
--- /dev/null
+++ b/math/deboor/progs/prog1.f
@@ -0,0 +1,45 @@
+chapter ii. runge example
+c  from  * a practical guide to splines *  by c. de boor
+      integer i,istep,j,k,n,nmk,nm1
+      real aloger,algerp,d(20),decay,dx,errmax,g,h,pnatx,step,tau(20),x
+      data step, istep /20., 20/
+      g(x) = 1./(1.+(5.*x)**2)
+      print 600
+  600 format(28h  n   max.error   decay exp.//)
+      decay = 0.
+      do 40 n=2,20,2
+c	 choose interpolation points  tau(1), ..., tau(n) , equally
+c	 spaced in (-1,1), and set  d(i) = g(tau(i)), i=1,...,n.
+	 nm1 = n-1
+	 h = 2./float(nm1)
+	 do 10 i=1,n
+	    tau(i) = float(i-1)*h - 1.
+   10	    d(i) = g(tau(i))
+c	 calculate the divided differences for the newton form.
+c
+	 do 20 k=1,nm1
+	    nmk = n-k
+	    do 20 i=1,nmk
+   20	       d(i) = (d(i+1)-d(i))/(tau(i+k)-tau(i))
+c
+c	 estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+	 do 30 i=2,n
+	    dx = (tau(i)-tau(i-1))/step
+	    do 30 j=1,istep
+	       x = tau(i-1) + float(j)*dx
+c	       evaluate interp.pol. by nested multiplication
+c
+	       pnatx = d(1)
+	       do 29 k=2,n
+   29		  pnatx = d(k) + (x-tau(k))*pnatx
+c
+   30	       errmax = amax1(errmax,abs(g(x)-pnatx))
+	 aloger = alog(errmax)
+	 if (n .gt. 2)	decay =
+     *	     (aloger - algerp)/alog(float(n)/float(n-2))
+	 algerp = aloger
+   40	 print 640,n,errmax,decay
+  640 format(i3,e12.4,f11.2)
+					stop
+      end
diff --git a/math/deboor/progs/prog10.f b/math/deboor/progs/prog10.f
new file mode 100644
index 00000000..dc580c99
--- /dev/null
+++ b/math/deboor/progs/prog10.f
@@ -0,0 +1,76 @@
+chapter xii, example 4. quasi-interpolant with good knots.
+c  from  * a practical guide to splines *  by c. de boor
+calls bsplpp(bsplvb)
+c
+      integer i,irate,istep,j,l,m,mmk,n,nlow,nhigh,nm1
+      real algerp,aloger,bcoef(22),break(20),c(4,20),decay,dg,ddg,x
+     *	  ,dtip1,dtip2,dx,errmax,g,h,pnatx,scrtch(4,4),step,t(26),taui
+c     istep and step = float(istep) specify point density for error det-
+c     ermination.
+      data step, istep /20., 20/
+c	 g  is the function to be approximated,  dg  is its first, and
+c	 ddg  its second derivative .
+      g(x) = sqrt(x+1.)
+      dg(x) = .5/g(x)
+      ddg(x) = -.5*dg(x)/(x+1.)
+      decay = 0.
+c      read in the exponent  irate  for the knot distribution and the
+c      lower and upper limit for the number  n	.
+      read 500,irate,nlow,nhigh
+  500 format(3i3)
+      print 600
+  600 format(28h  n   max.error   decay exp./)
+c			       loop over  n = dim( spline(4,t) ) - 2 .
+c	       n  is chosen as the parameter in order to afford compar-
+c	       ison with examples 2 and 3  in which cubic spline interp-
+c	       olation at  n  data points was used .
+      do 40 n=nlow,nhigh,2
+	 nm1 = n-1
+	 h = 1./float(nm1)
+	 m = n+2
+	 mmk = m-4
+	 do 5 i=1,4
+	    t(i) = -1.
+    5	    t(m+i) = 1.
+c		  interior knots are equidistributed with respect to the
+c		  function  (x + 1)**(1/irate) .
+	 do 6 i=1,mmk
+    6	    t(i+4) = 2.*(float(i)*h)**irate - 1.
+c					    construct quasi-interpolant.
+c						 bcoef(1) = g(-1.) = 0.
+	 bcoef(1) = 0.
+	 dtip2 = t(5) - t(4)
+	 taui = t(5)
+c			    special choice of  tau(2)  to avoid infinite
+c			    derivatives of  g  at left endpoint .
+	 bcoef(2) = g(taui) - 2.*dtip2*dg(taui)/3.
+     *		   + dtip2**2*ddg(taui)/6.
+	 do 15 i=3,m
+	    taui = t(i+2)
+	    dtip1 = dtip2
+	    dtip2 = t(i+3) - t(i+2)
+c				       formula xii(30) of text is used .
+   15	    bcoef(i) = g(taui) + (dtip2-dtip1)*dg(taui)/3.
+     *		     - dtip1*dtip2*ddg(taui)/6.
+c					  convert to pp-representation .
+	 call bsplpp(t,bcoef,m,4,scrtch,break,c,l)
+c			     estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+c					      loop over cubic pieces ...
+	 do 30 i=1,l
+	    dx = (break(i+1)-break(i))/step
+c			error is calculated at	istep  points per 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
diff --git a/math/deboor/progs/prog11.f b/math/deboor/progs/prog11.f
new file mode 100644
index 00000000..2deea951
--- /dev/null
+++ b/math/deboor/progs/prog11.f
@@ -0,0 +1,69 @@
+chapter xiii, example 1. a large norm amplifies noise.
+c  from  * a practical guide to splines *  by c. de boor
+calls splint(banfac/slv,bsplvb),bsplpp(bsplvb*),ppvalu(interv),round
+c     an initially uniform data point distribution of  n  points is
+c  changed  i t e r m x  times by moving the  j c l o s e -th data point
+c  toward its left neighbor, cutting the distance between the two by a
+c  factor of  r a t e  each time . to demonstrate the corresponding
+c  increase in the norm of cubic spline interpolation at these data
+c  points, the data are taken from a cubic polynomial (which would be
+c  interpolated exactly) but with noise of size  s i z e  added. the re-
+c  sulting noise or error in the interpolant (compared to the cubic)
+c  gives the norm or noise amplification factor and is printed out
+c  together with the diminishing distance  h  between the two data
+c  points.
+c     parameter nmax=200,nmaxt4=800,nmaxt7=1400,nmaxp4=204
+      integer i,iflag,istep,iter,itermx,j,jclose,l,n,nm1
+c     real amax,bcoef(nmax),break(nmax),coef(nmaxt4),dx,fltnm1,fx
+c    * ,gtau(nmax),h,rate,scrtch(nmaxt7),size,step,t(nmaxp4),tau(nmax),x
+      real amax,bcoef(200),break(200),coef(800),dx,fltnm1,fx
+     * ,gtau(200),h,rate,scrtch(1400),size,step,t(204),tau(200),x
+      real ppvalu,round,g
+      common /rount/ size
+      data step, istep / 20., 20 /
+c					   function to be interpolated .
+      g(x) = 1.+x*(1.+x*(1.+x))
+      read 500,n,itermx,jclose,size,rate
+  500 format(3i3/e10.3/e10.3)
+      print 600,size
+  600 format(16h size of noise =,e8.2//
+     *	    25h    h	       max.error)
+c					start with uniform data points .
+      nm1 = n - 1
+      fltnm1 = float(nm1)
+      do 10 i=1,n
+   10	 tau(i) = float(i-1)/fltnm1
+c		     set up knot sequence for not-a-knot end condition .
+      do 11 i=1,4
+	 t(i) = tau(1)
+   11	 t(n+i) = tau(n)
+      do 12 i=5,n
+   12	 t(i) = tau(i-2)
+c
+      do 100 iter=1,itermx
+	 do 21 i=1,n
+   21	    gtau(i) = round(g(tau(i)))
+	 call splint ( tau, gtau, t, n, 4, scrtch, bcoef, iflag )
+					go to (24,23),iflag
+   23	 print 623
+  623	 format(27h something wrong in splint.)
+					stop
+   24	 call bsplpp ( t, bcoef, n, 4, scrtch, break, coef, l )
+c				     calculate max.interpolation error .
+	 amax = 0.
+	 do 30 i=4,n
+	    dx = (break(i-2) - break(i-3))/step
+	    do 25 j=2,istep
+	       x = break(i-2) - dx*float(j-1)
+	       fx = ppvalu(break,coef,l,4,x,0)
+   25	       amax = amax1(amax,abs(fx-g(x)))
+   30	    continue
+	 h = tau(jclose) - tau(jclose-1)
+	 print 630,h,amax
+  630	 format(e9.2,e15.3)
+c		  move tau(jclose) toward its left neighbor so as to cut
+c		  their distance by a factor of  rate .
+	 tau(jclose) = (tau(jclose) + (rate-1.)*tau(jclose-1))/rate
+  100	 continue
+					stop
+      end
diff --git a/math/deboor/progs/prog12.f b/math/deboor/progs/prog12.f
new file mode 100644
index 00000000..6efe8614
--- /dev/null
+++ b/math/deboor/progs/prog12.f
@@ -0,0 +1,70 @@
+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*)
+c     parameter nmax=20,nmaxp6=26,nmaxp2=22,nmaxt7=140
+      integer i,istep,iter,itermx,n,nhigh,nmk,nlow
+c     real algerp,aloger,bcoef(nmaxp2),break(nmax),decay,dx,errmax,
+c    *	  c(4,nmax),g,gtau(nmax),h,pnatx,scrtch(nmaxt7),step,t(nmaxp6),
+c    *	  tau(nmax),tnew(nmax)
+      real algerp,aloger,bcoef(22),break(20),decay,dx,errmax,
+     *	  c(4,20),g,gtau(20),h,pnatx,scrtch(140),step,t(26),
+     *	  tau(20),tnew(20),x
+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
diff --git a/math/deboor/progs/prog13.f b/math/deboor/progs/prog13.f
new file mode 100644
index 00000000..45e6362f
--- /dev/null
+++ b/math/deboor/progs/prog13.f
@@ -0,0 +1,77 @@
+chapter xiii, example 2m. cubic spline interpolant at knot averages
+c			 with good knots. modified around label 4.
+c  from  * a practical guide to splines *  by c. de boor
+calls splint(banfac/slv,bsplvb),newnot,bsplpp(bsplvb*)
+c     parameter nmax=20,nmaxp6=26,nmaxp2=22,nmaxt7=140
+      integer i,istep,iter,itermx,n,nhigh,nmk,nlow
+c     real algerp,aloger,bcoef(nmaxp2),break(nmax),decay,dx,errmax,
+c    *	  c(4,nmax),g,gtau(nmax),h,pnatx,scrtch(nmaxt7),step,t(nmaxp6),
+c    *	  tau(nmax),tnew(nmax)
+      real algerp,aloger,bcoef(22),break(20),decay,dx,errmax,
+     *	  c(4,20),g,gtau(20),h,pnatx,scrtch(140),step,t(26),
+     *	  tau(20),tnew(20),x
+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
+	 if (n .eq. nlow)		go to 4
+	 call newnot ( break, c, l, 4, tnew, l+2, scrtch )
+	 l = l + 2
+	 t(5+l) = 1.
+	 t(6+l) = 1.
+	 iter = 1
+					go to 15
+    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
diff --git a/math/deboor/progs/prog14.f b/math/deboor/progs/prog14.f
new file mode 100644
index 00000000..03804f47
--- /dev/null
+++ b/math/deboor/progs/prog14.f
@@ -0,0 +1,37 @@
+chapter xiii, example 3 . test of optimal spline interpolation routine
+c			  on titanium heat data .
+c  from  * a practical guide to splines *  by c. de boor
+calls titand,splopt(bsplvb,banfac/slv),splint(*),bvalue(interv)
+c
+c     lenscr = (n-k)(2k+3)+5k+3  is the length of scrtch required in
+c				 splopt .
+c     parameter n=12,ntitan=49,k=5,npk=17,lenscr=119
+c     integer i,iflag,ipick(n),ipicki,lx,nmk
+c     real a(n),gtitan(ntitan),gtau(ntitan),scrtch(lenscr),t(npk),tau(n)
+c    *	  ,x(ntitan)
+      real bvalue
+      integer i,iflag,ipick(12),ipicki,lx,nmk
+      real a(12),gtitan(49),gtau(49),scrtch(119),t(17),tau(12)
+     *	  ,x(49)
+      data n,k /12,5/
+      data ipick /1,5,11,21,27,29,31,33,35,40,45,49/
+      call titand ( x, gtitan, lx )
+      do 10 i=1,n
+	 ipicki = ipick(i)
+	 tau(i) = x(ipicki)
+   10	 gtau(i) = gtitan(ipicki)
+      call splopt ( tau, n, k, scrtch, t, iflag )
+      if (iflag .gt. 1) 		stop
+      call splint ( tau, gtau, t, n, k, scrtch, a, iflag )
+      if (iflag .gt. 1) 		stop
+      do 20 i=1,lx
+	 gtau(i) = bvalue ( t, a, n, k, x(i), 0 )
+   20	 scrtch(i) = gtitan(i) - gtau(i)
+      print 620,(i,x(i),gtitan(i),gtau(i),scrtch(i),i=1,lx)
+  620 format(41h  i, data point, data, interpolant, error//
+     2	       (i3,f8.0,f10.4,f9.4,e11.3))
+      nmk = n-k
+      print 621,(i,t(k+i),i=1,nmk)
+  621 format(///16h optimal knots =/(i5,f15.9))
+					stop
+      end
diff --git a/math/deboor/progs/prog15.f b/math/deboor/progs/prog15.f
new file mode 100644
index 00000000..2158931e
--- /dev/null
+++ b/math/deboor/progs/prog15.f
@@ -0,0 +1,48 @@
+chapter xiv, example 1. cubic smoothing spline
+c  from  * a practical guide to splines *  by c. de boor
+calls  bsplpp(bsplvb), ppvalu(interv), smooth(setupq,chol1d)
+c     values from a cubic b-spline are rounded to  n d i g i t	places
+c  after the decimal point, then smoothed via  s m o o t h  for
+c  various values of the control parameter  s .
+      integer i,is,j,l,lsm,ndigit,npoint,ns
+      real a(61,4),bcoef(7),break(5),coef(4,4),coefsm(4,60),dely,dy(61)
+     *	  ,s(20),scrtch(427),sfp,t(11),tenton,x(61),y(61)
+      real ppvalu,smooth
+      equivalence (scrtch,coefsm)
+      data t /4*0.,1.,3.,4.,4*6./
+      data bcoef /3*0.,1.,3*0./
+      call bsplpp(t,bcoef,7,4,scrtch,break,coef,l)
+      npoint = 61
+      read 500,ndigit,ns,(s(i),i=1,ns)
+  500 format(2i3/(e10.4))
+      print 600,ndigit
+  600 format(24h exact values rounded to,i2,21h digits after decimal
+     *	     ,7h point./)
+      tenton = 10.**ndigit
+      dely = .5/tenton
+      do 10 i=1,npoint
+	 x(i) = .1*float(i-1)
+	 y(i) = ppvalu(break,coef,l,4,x(i),0)
+	 y(i) = float(int(y(i)*tenton + .5))/tenton
+   10	 dy(i) = dely
+      do 15 i=1,npoint,5
+	 do 15 j=1,4
+   15	    a(i,j) = ppvalu(break,coef,l,4,x(i),j-1)
+      print 615,(i,(a(i,j),j=1,4),i=1,npoint,5)
+  615 format(52h value and derivatives of noisefree function at some
+     *	    ,7h points/(i4,4e15.7))
+      do 20 is=1,ns
+	 sfp = smooth ( x, y, dy, npoint, s(is), scrtch, a )
+	 lsm = npoint - 1
+	 do 16 i=1,lsm
+	    do 16 j=1,4
+   16	       coefsm(j,i) = a(i,j)
+	 do 18 i=1,npoint,5
+	    do 18 j=1,4
+   18	       a(i,j) = ppvalu(x,coefsm,lsm,4,x(i),j-1)
+   20	 print 620, s(is), sfp, (i,(a(i,j),j=1,4),i=1,npoint,5)
+  620 format(15h prescribed s =,e10.3,23h, s(smoothing spline) =,e10.3/
+     *	     54h value and derivatives of smoothing spline at corresp.
+     *	    ,7h points/(i4,4e15.7))
+					stop
+      end
diff --git a/math/deboor/progs/prog16.f b/math/deboor/progs/prog16.f
new file mode 100644
index 00000000..4dd39911
--- /dev/null
+++ b/math/deboor/progs/prog16.f
@@ -0,0 +1,115 @@
+c  main program for least-squares approximation by splines
+c  from  * a practical guide to splines *  by c. de boor
+calls setdat,l2knts,l2appr(bsplvb,bchfac,bchslv),bsplpp(bsplvb*)
+c     ,l2err(ppvalu(interv)),ppvalu*,newnot
+c
+c  the program, though ostensibly written for l2-approximation, is typ-
+c  ical for programs constructing a pp approximation to a function gi-
+c  ven in some sense. the subprogram  l 2 a p p r , for instance, could
+c  easily be replaced by one carrying out interpolation or some other
+c  form of approximation.
+c
+c******  i n p u t  ******
+c  is expected in  s e t d a t	(quo vide), specifying both the data to
+c  be approximated and the order and breakpoint sequence of the pp ap-
+c  proximating function to be used. further,  s e t d a t  is expected
+c  to  t e r m i n a t e  the run (for lack of further input or because
+c   i c o u n t  has reached a critical value).
+c     the number  n t i m e s  is read in in the main program. it speci
+c  fies the number of passes through the knot improvement algorithm in
+c  n e w n o t	to be made. also,  a d d b r k	is read in to specify
+c  that, on the average, addbrk knots are to be added per pass through
+c  newnot. for example,  addbrk = .34  would cause a knot to be added
+c  every third pass (as long as  ntimes .lt. 50).
+c
+c******  p r i n t e d	o u t p u t  ******
+c  is governed by the three print control hollerith strings
+c  p r b c o  = 2hon  gives printout of b-spline coeffs. of approxim.
+c  p r p c o  = 2hon  gives printout of pp repr. of approximation.
+c  p r f u n  = 2hon  gives printout of approximation and error at
+c		      every data point.
+c  the order  k , the number of pieces	l, and the interior breakpoints
+c  are always printed out as are (in l2err) the mean, mean square, and
+c  maximum errors in the approximation.
+c
+c     parameter lpkmax=100,ntmax=200,ltkmax=2000
+      integer i,icount,ii,j,k,l,lbegin,lnew,ll,n,nt,ntimes,ntau
+     *	     ,on,prbco,prfun,prpco
+c     real addbrk,bcoef(lpkmax),break,coef,gtau,q(ltkmax),scrtch(ntmax)
+c    *	  ,t(ntmax),tau,totalw,weight
+c     common / data / ntau, tau(ntmax),gtau(ntmax),weight(ntmax),totalw
+      real addbrk,bcoef(100),break,coef,gtau,q(2000),scrtch(200)
+     *	  ,t(200),tau,totalw,weight
+      common / data / ntau, tau(200),gtau(200),weight(200),totalw
+c     common /data/ also occurs in setdat, l2appr and l2err. it is ment-
+c     ioned here only because it might otherwise become undefined be-
+c     tween calls to those subroutines.
+c     common /approx/ break(lpkmax),coef(ltkmax),l,k
+      common /approx/ break(100),coef(2000),l,k
+c     common /approx/ also occurs in setdat and l2err.
+      data on /2hon /
+c
+      icount = 0
+c	 i c o u n t  provides communication with the data-input-and-
+c     termination routine  s e t d a t . it is initialized to  0  to
+c     signal to setdat when it is being called for the first time. after
+c     that, it is up to setdat to use icount for keeping track of the
+c     passes through setdat .
+c
+c     information about the function to be approximated and order and
+c     breakpoint sequence of the approximating pp functions is gathered
+c     by a
+    1 call setdat(icount)
+c
+c     breakpoints are translated into knots, and the number  n	of
+c     b-splines to be used is obtained by a
+      call l2knts ( break, l, k, t, n )
+c
+c     the integer  n t i m e s	and the real  a d d b r k  are requested
+c     as well as the print controls  p r b c o ,  p r p c o  and
+c     p r f u n .  ntimes  passes  are made through the subroutine new-
+c     not, with an increase of	addbrk	knots for every pass .
+      print 600
+  600 format(52h ntimes,addbrk , prbco,prpco,prfun =? (i3,f10.5/3a2))
+      read 500,ntimes,addbrk,prbco,prpco,prfun
+  500 format(i3,f10.5/3a2)
+c
+      lbegin = l
+      nt = 1
+c	 the b-spline coeffs.  b c o e f  of the l2-approx. are obtain-
+c	 ed by a
+   10	 call l2appr ( t, n, k, q, scrtch, bcoef )
+	 if (prbco .eq. on)  print 609, (bcoef(i),i=1,n)
+  609	 format(//22h b-spline coefficients/(5e16.9))
+c
+c	 convert the b-repr. of the approximation to pp repr.
+	 call bsplpp ( t, bcoef, n, k, q, break, coef, l )
+	 print 610, k, l, (break(ll),ll=2,l)
+  610	 format(//34h approximation by splines of order,i3,4h on ,
+     *	       i3,25h intervals. breakpoints -/(5e16.9))
+	 if (prpco .ne. on)		go to 15
+	 print 611
+  611	 format(/36h pp-representation for approximation)
+	 do 12 i=1,l
+	    ii = (i-1)*k
+   12	    print 613,break(i),(coef(ii+j),j=1,k)
+  613	 format(f9.3,5e16.9/(11x,5e16.9))
+c
+c	 compute and print out various error norms by a
+   15	 call l2err ( prfun, scrtch, q )
+c
+c	 if newnot has been applied less than  n t i m e s  times, try
+c	 it again to obtain, from the current approx. a possibly improv-
+c	 ed sequence of breakpoints with  addbrk  more breakpoints (on
+c	 the average) than the current approximation has.
+c	    if only an increase in breakpoints is wanted, without the
+c	 adjustment that newnot provides, a fake newnot routine could be
+c	 used here which merely returns the breakpoints for  l n e w
+c	 equal intervals .
+	 if (nt .ge. ntimes)		go to 1
+	 lnew = lbegin + float(nt)*addbrk
+	 call newnot (break, coef, l, k, scrtch, lnew, t )
+	 call l2knts ( scrtch, lnew, k, t, n )
+	 nt = nt + 1
+					go to 10
+      end
diff --git a/math/deboor/progs/prog17.f b/math/deboor/progs/prog17.f
new file mode 100644
index 00000000..1008c3a1
--- /dev/null
+++ b/math/deboor/progs/prog17.f
@@ -0,0 +1,10 @@
+c  main program for example in chapter xv.
+c  from  * a practical guide to splines *  by c. de boor
+c  solution of a second order nonlinear two point boundary value
+c  problem  on (0., 1.) , by collocation with pp functions having  4
+c  pieces of order  6 .  2  passes through newnot are to be made,
+c  without any knots being added. newton iteration is to be stopped
+c  when two iterates agree to  6  decimal places.
+       call colloc(0.,1.,4,6,2,0.,1.e-6)
+					stop
+       end
diff --git a/math/deboor/progs/prog18.f b/math/deboor/progs/prog18.f
new file mode 100644
index 00000000..fad18c62
--- /dev/null
+++ b/math/deboor/progs/prog18.f
@@ -0,0 +1,35 @@
+chapter xvi, example 2, taut spline interpolation to the
+c  titanium heat data of example xiii.3.
+c  from  * a practical guide to splines *  by c. de boor
+calls titand,tautsp,ppvalu(interv)
+      integer i,iflag,ipick(12),ipicki,k,l,lx,n,npoint
+      real break(122),coef(4,22),gamma,gtau(49),gtitan(49),plotf(201)
+     *	  ,plott(201),plotts(201),scrtch(119),step,tau(12),x(49)
+      real ppvalu
+      data n,ipick /12,1,5,11,21,27,29,31,33,35,40,45,49/
+      data k,npoint /4,201/
+      call titand ( x, gtitan, lx )
+      do 10 i=1,n
+	 ipicki = ipick(i)
+	 tau(i) = x(ipicki)
+   10	 gtau(i) = gtitan(ipicki)
+      call tautsp(tau,gtau,n,0.,scrtch,break,coef,l,k,iflag)
+      if (iflag .gt. 1) 		stop
+      step = (tau(n) - tau(1))/float(npoint-1)
+      do 20 i=1,npoint
+	 plott(i) = tau(1) + step*float(i-1)
+   20	 plotf(i) = ppvalu(break,coef,l,k,plott(i),0)
+    1 print 601
+  601 format(18h gamma = ? (f10.3))
+      read 500,gamma
+  500 format(f10.3)
+      if (gamma .lt. 0.)		stop
+      call tautsp(tau,gtau,n,gamma,scrtch,break,coef,l,k,iflag)
+      if (iflag .gt. 1) 		stop
+      do 30 i=1,npoint
+   30	 plotts(i) = ppvalu(break,coef,l,k,plott(i),0)
+      print 602,gamma,(plott(i),plotf(i),plotts(i),i=1,npoint)
+  602 format(/42h cubic spline vs. taut spline with gamma =,f6.3//
+     *	    (f15.7,2e20.8))
+					go to 1
+      end
diff --git a/math/deboor/progs/prog19.f b/math/deboor/progs/prog19.f
new file mode 100644
index 00000000..6595b76c
--- /dev/null
+++ b/math/deboor/progs/prog19.f
@@ -0,0 +1,58 @@
+chapter xvi, example 3. two parametrizations of some data
+c  from  * a practical guide to splines *  by c. de boor
+calls  splint(bsplvb,banfac/slv),bsplpp(bsplvb*),banslv*,ppvalu(interv)
+c     parameter k=4,kpkm1=7, n=8,npk=12, npiece=6,npoint=21
+c     integer i,icount,iflag,kp1,l
+c     real bcoef(n),break(npiece),ds,q(n,kpkm1),s(n),scrtch(k,k)
+c    *	  ,ss,t(npk),x(n),xcoef(k,npiece),xx(npoint),y(n)
+c    *	  ,ycoef(k,npiece),yy(npoint)
+      integer i,icount,iflag,k,kpkm1,kp1,l,n,npoint
+      real bcoef(8),break(6),ds,q(8,7),s(8),scrtch(4,4)
+     *	  ,ss,t(12),x(8),xcoef(4,6),xx(21),y(8)
+     *	  ,ycoef(4,6),yy(21)
+      real ppvalu
+      data k,kpkm1,n,npoint /4,7,8,21/
+      data x /0.,.1,.2,.3,.301,.4,.5,.6/
+c  *** compute y-component and set 'natural' parametrization
+      do 1 i=1,n
+	 y(i) = (x(i)-.3)**2
+    1	 s(i) = x(i)
+      print 601
+  601 format(26h 'natural' parametrization/6x,1hx,11x,1hy)
+      icount = 1
+c  *** convert data abscissae to knots. note that second and second
+c     last data abscissae are not knots.
+    5 do 6 i=1,k
+	 t(i) = s(1)
+    6	 t(n+i) = s(n)
+      kp1 = k+1
+      do 7 i=kp1,n
+    7	 t(i) = s(i+2-k)
+c  *** interpolate to x-component
+      call splint(s,x,t,n,k,q,bcoef,iflag)
+      call bsplpp(t,bcoef,n,k,scrtch,break,xcoef,l)
+c  *** interpolate to y-component. since data abscissae and knots are
+c     the same for both components, we only need to use backsubstitution
+      do 10 i=1,n
+   10	 bcoef(i) = y(i)
+      call banslv(q,kpkm1,n,k-1,k-1,bcoef)
+      call bsplpp(t,bcoef,n,k,scrtch,break,ycoef,l)
+c  *** evaluate curve at some points near the potential trouble spot,
+c     the fourth and fifth data points.
+      ss = s(3)
+      ds = (s(6)-s(3))/float(npoint-1)
+      do 20 i=1,npoint
+	 xx(i) = ppvalu(break,xcoef,l,k,ss,0)
+	 yy(i) = ppvalu(break,ycoef,l,k,ss,0)
+   20	 ss = ss + ds
+      print 620,(xx(i),yy(i),i=1,npoint)
+  620 format(2f12.7)
+      if (icount .ge. 2)		stop
+c  *** now repeat the whole process with uniform parametrization
+      icount = icount + 1
+      do 30 i=1,n
+   30	 s(i) = float(i)
+      print 630
+  630 format(/26h 'uniform' parametrization/6x,1hx,11x,1hy)
+					go to 5
+      end
diff --git a/math/deboor/progs/prog2.f b/math/deboor/progs/prog2.f
new file mode 100644
index 00000000..1da20d52
--- /dev/null
+++ b/math/deboor/progs/prog2.f
@@ -0,0 +1,48 @@
+chapter iv. runge example, with cubic hermite interpolation
+c  from  * a practical guide to splines *  by c. de boor
+      integer i,istep,j,n,nm1
+      real aloger,algerp,c(4,20),decay,divdf1,divdf3,dtau,dx,errmax,g,h
+     *	  ,pnatx,step,tau(20),x
+      data step, istep /20., 20/
+      g(x) = 1./(1.+(5.*x)**2)
+      print 600
+  600 format(28h  n   max.error   decay exp.//)
+      decay = 0.
+      do 40 n=2,20,2
+c	 choose interpolation points  tau(1), ..., tau(n) , equally
+c	 spaced in (-1,1), and set  c(1,i) = g(tau(i)), c(2,i) =
+c	 gprime(tau(i)) = -50.*tau(i)*g(tau(i))**2, i=1,...,n.
+	 nm1 = n-1
+	 h = 2./float(nm1)
+	 do 10 i=1,n
+	    tau(i) = float(i-1)*h - 1.
+	    c(1,i) = g(tau(i))
+   10	    c(2,i) = -50.*tau(i)*c(1,i)**2
+c	 calculate the coefficients of the polynomial pieces
+c
+	 do 20 i=1,nm1
+	    dtau = tau(i+1) - tau(i)
+	    divdf1 = (c(1,i+1) - c(1,i))/dtau
+	    divdf3 = c(2,i) + c(2,i+1) - 2.*divdf1
+	    c(3,i) = (divdf1 - c(2,i) - divdf3)/dtau
+   20	    c(4,i) = (divdf3/dtau)/dtau
+c
+c	 estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+	 do 30 i=2,n
+	    dx = (tau(i)-tau(i-1))/step
+	    do 30 j=1,istep
+	       h = float(j)*dx
+c	       evaluate (i-1)st cubic piece
+c
+	       pnatx = c(1,i-1)+h*(c(2,i-1)+h*(c(3,i-1)+h*c(4,i-1)))
+c
+   30	       errmax = amax1(errmax,abs(g(tau(i-1)+h)-pnatx))
+	 aloger = alog(errmax)
+	 if (n .gt. 2)	decay =
+     *	     (aloger - algerp)/alog(float(n)/float(n-2))
+	 algerp = aloger
+   40	 print 640,n,errmax,decay
+  640 format(i3,e12.4,f11.2)
+					stop
+      end
diff --git a/math/deboor/progs/prog20.f b/math/deboor/progs/prog20.f
new file mode 100644
index 00000000..78143095
--- /dev/null
+++ b/math/deboor/progs/prog20.f
@@ -0,0 +1,62 @@
+chapter xvii, example 2. bivariate spline interpolation
+c  from  * a practical guide to splines *  by c. de boor
+calls  spli2d(bsplvb,banfac/slv),interv,bvalue(bsplvb*,interv*)
+c     parameter nx=7,kx=3,ny=6,ky=4
+c     integer i,iflag,j,jj,lefty,mflag,
+c     real bcoef(nx,ny),taux(nx),tauy(ny),tx(10),ty(10)
+c    *	       ,work1(nx,ny),work2(nx),work3(42)
+      integer i,iflag,j,jj,kp1,kx,ky,lefty,mflag,nx,ny
+      real bcoef(7,6),taux(7),tauy(6),tx(10),ty(10)
+     *	       ,work1(7,6),work2(7),work3(42)
+      real bvalue,g,x,y
+      data nx,kx,ny,ky /7,3,6,4/
+      g(x,y) = amax1(x-3.5,0.)**2 + amax1(y-3.,0.)**3
+c *** set up data points and knots
+c     in x, interpolate between knots by parabolic splines, using
+c     not-a-knot end condition
+      do 10 i=1,nx
+   10	 taux(i) = float(i)
+      do 11 i=1,kx
+	 tx(i) = taux(1)
+   11	 tx(nx+i) = taux(nx)
+      kp1 = kx+1
+      do 12 i=kp1,nx
+   12	 tx(i) = (taux(i-kx+1) + taux(i-kx+2))/2.
+c     in y, interpolate at knots by cubic splines, using not-a-knot
+c     end condition
+      do 20 i=1,ny
+   20	 tauy(i) = float(i)
+      do 21 i=1,ky
+	 ty(i) = tauy(1)
+   21	 ty(ny+i) = tauy(ny)
+      kp1 = ky+1
+      do 22 i=kp1,ny
+   22	 ty(i) = tauy(i-ky+2)
+c  *** generate and print out function values
+      print 620,(tauy(i),i=1,ny)
+  620 format(11h given data//6f12.1)
+      do 32 i=1,nx
+	 do 31 j=1,ny
+   31	    bcoef(i,j) = g(taux(i),tauy(j))
+   32	 print 632,taux(i),(bcoef(i,j),j=1,ny)
+  632 format(f5.1,6e12.5)
+c
+c  *** construct b-coefficients of interpolant
+c
+      call spli2d(taux,bcoef,tx,nx,kx,ny,work2,work3,work1,iflag)
+      call spli2d(tauy,work1,ty,ny,ky,nx,work2,work3,bcoef,iflag)
+c
+c  *** evaluate interpolation error at mesh points and print out
+      print 640,(tauy(j),j=1,ny)
+  640 format(//20h interpolation error//6f12.1)
+      do 45 j=1,ny
+	 call interv(ty,ny,tauy(j),lefty,mflag)
+	 do 45 i=1,nx
+	    do 41 jj=1,ky
+   41	       work2(jj)=bvalue(tx,bcoef(1,lefty-ky+jj),nx,kx,taux(i),0)
+   45	    work1(i,j) = g(taux(i),tauy(j)) -
+     *			bvalue(ty(lefty-ky+1),work2,ky,ky,tauy(j),0)
+      do 46 i=1,nx
+   46	 print 632,taux(i),(work1(i,j),j=1,ny)
+					stop
+      end
diff --git a/math/deboor/progs/prog21.f b/math/deboor/progs/prog21.f
new file mode 100644
index 00000000..6d8ca62a
--- /dev/null
+++ b/math/deboor/progs/prog21.f
@@ -0,0 +1,74 @@
+chapter xvii, example 3. bivariate spline interpolation
+c  followed by conversion to pp-repr and evaluation
+c  from  * a practical guide to splines *  by c. de boor
+calls  spli2d(bsplvb,banfac/slv),interv,bspp2d(bsplvb*),ppvalu(interv*)
+c     parameter nx=7,kx=3,ny=6,ky=4
+c     integer i,iflag,j,jj,lefty,lx,ly,mflag
+c     real bcoef(nx,ny),breakx(6),breaky(4),coef(kx,5,ky,3),taux(nx)
+c    *	  ,tauy(ny),tx(10),ty(10)
+c    *	       ,work1(nx,ny),work2(nx),work3(42)
+c    *	       ,work4(kx,kx,ny),work5(ny,kx,5),work6(ky,ky,21)
+      integer i,iflag,j,jj,kp1,kx,ky,lefty,lx,ly,mflag,nx,ny
+      real bcoef(7,6),breakx(6),breaky(4),coef(3,5,4,3),taux(7)
+     *	  ,tauy(6),tx(10),ty(10)
+     *	       ,work1(7,6),work2(7),work3(42)
+     *	       ,work4(3,3,6),work5(6,3,5),work6(4,4,21)
+      real ppvalu,g,x,y
+      data nx,kx,ny,ky / 7,3,6,4 /
+c     note that , with the above parameters, lx=5, ly = 3
+      equivalence (work4,work6)
+      g(x,y) = amax1(x-3.5,0.)**2 + amax1(y-3.,0.)**3
+c *** set up data points and knots
+c     in x, interpolate between knots by parabolic splines, using
+c     not-a-knot end condition
+      do 10 i=1,nx
+   10	 taux(i) = float(i)
+      do 11 i=1,kx
+	 tx(i) = taux(1)
+   11	 tx(nx+i) = taux(nx)
+      kp1 = kx+1
+      do 12 i=kp1,nx
+   12	 tx(i) = (taux(i-kx+1) + taux(i-kx+2))/2.
+c     in y, interpolate at knots by cubic splines, using not-a-knot
+c     end condition
+      do 20 i=1,ny
+   20	 tauy(i) = float(i)
+      do 21 i=1,ky
+	 ty(i) = tauy(1)
+   21	 ty(ny+i) = tauy(ny)
+      kp1 = ky+1
+      do 22 i=kp1,ny
+   22	 ty(i) = tauy(i-ky+2)
+c  *** generate and print out function values
+      print 620,(tauy(i),i=1,ny)
+  620 format(11h given data//6f12.1)
+      do 32 i=1,nx
+	 do 31 j=1,ny
+   31	    bcoef(i,j) = g(taux(i),tauy(j))
+   32	 print 632,taux(i),(bcoef(i,j),j=1,ny)
+  632 format(f5.1,6e12.5)
+c
+c  *** construct b-coefficients of interpolant
+c
+      call spli2d(taux,bcoef,tx,nx,kx,ny,work2,work3,work1,iflag)
+      call spli2d(tauy,work1,ty,ny,ky,nx,work2,work3,bcoef,iflag)
+c
+c  *** convert to pp-representation
+c
+      call bspp2d(tx,bcoef,nx,kx,ny,work4,breakx,work5,lx)
+      call bspp2d(ty,work5,ny,ky,lx*kx,work6,breaky,coef,ly)
+c
+c  *** evaluate interpolation error at mesh points and print out
+      print 640,(tauy(j),j=1,ny)
+  640 format(//20h interpolation error//6f12.1)
+      do 45 j=1,ny
+	 call interv(breaky,ly,tauy(j),lefty,mflag)
+	 do 45 i=1,nx
+	    do 41 jj=1,ky
+   41	     work2(jj)=ppvalu(breakx,coef(1,1,jj,lefty),lx,kx,taux(i),0)
+   45	    work1(i,j) = g(taux(i),tauy(j)) -
+     *	    ppvalu(breaky(lefty),work2,1,ky,tauy(j),0)
+      do 46 i=1,nx
+   46	 print 632,taux(i),(work1(i,j),j=1,ny)
+					stop
+      end
diff --git a/math/deboor/progs/prog3.f b/math/deboor/progs/prog3.f
new file mode 100644
index 00000000..337fe7ef
--- /dev/null
+++ b/math/deboor/progs/prog3.f
@@ -0,0 +1,44 @@
+chapter ix. example comparing the b-representation of a cubic  f  with
+c  its values at knot averages.
+c  from  * a practical guide to splines *  by c. de boor
+c
+      integer i,id,j,jj,n,nm4
+      real bcoef(23),d(4),d0(4),dtip1,dtip2,f(23),t(27),tave(23),x
+c	     the taylor coefficients at  0  for the polynomial	f  are
+      data d0 /-162.,99.,-18.,1./
+c
+c		  set up knot sequence in the array  t .
+      n = 13
+      do 5 i=1,4
+	 t(i) = 0.
+    5	 t(n+i) = 10.
+      nm4 = n-4
+      do 6 i=1,nm4
+    6	 t(i+4) = float(i)
+c
+      do 50 i=1,n
+c	 use nested multiplication to get taylor coefficients  d  at
+c		t(i+2)	from those at  0 .
+	 do 20 j=1,4
+   20	    d(j) = d0(j)
+	 do 21 j=1,3
+	    id = 4
+	    do 21 jj=j,3
+	       id = id-1
+   21	       d(id) = d(id) + d(id+1)*t(i+2)
+c
+c		 compute b-spline coefficients by formula (9).
+	 dtip1 = t(i+2) - t(i+1)
+	 dtip2 = t(i+3) - t(i+2)
+	 bcoef(i) = d(1) + (d(2)*(dtip2-dtip1)-d(3)*dtip1*dtip2)/3.
+c
+c		   evaluate  f	at corresp. knot average.
+	 tave(i) = (t(i+1) + t(i+2) + t(i+3))/3.
+	 x = tave(i)
+   50	 f(i) = d0(1) + x*(d0(2) + x*(d0(3) + x*d0(4)))
+c
+      print 650, (i,tave(i), f(i), bcoef(i),i=1,n)
+  650 format(45h  i   tave(i)	   f at tave(i)      bcoef(i)//
+     *	     (i3,f10.5,2f16.5))
+					stop
+      end
diff --git a/math/deboor/progs/prog4.f b/math/deboor/progs/prog4.f
new file mode 100644
index 00000000..8fc73304
--- /dev/null
+++ b/math/deboor/progs/prog4.f
@@ -0,0 +1,35 @@
+chapter x. example 1. plotting some b-splines
+c  from  * a practical guide to splines *  by c. de boor
+calls  bsplvb, interv
+      integer i,j,k,left,leftmk,mflag,n,npoint
+      real dx,t(10),values(7),x,xl
+c  dimension, order and knot sequence for spline space are specified...
+      data n,k /7,3/, t /3*0.,2*1.,3.,4.,3*6./
+c  b-spline values are initialized to 0., number of evaluation points...
+      data values /7*0./, npoint /31/
+c  set leftmost evaluation point  xl , and spacing  dx	to be used...
+      xl = t(k)
+      dx = (t(n+1)-t(k))/float(npoint-1)
+c
+      print 600,(i,i=1,5)
+  600 format(4h1  x,8x,5(1hb,i1,3h(x),7x))
+c
+      do 10 i=1,npoint
+	 x = xl + float(i-1)*dx
+c		      locate  x  with respect to knot array  t .
+	 call interv ( t, n, x, left, mflag )
+	 leftmk = left - k
+c	    get b(i,k)(x)  in  values(i) , i=1,...,n .	k  of these,
+c	 viz.  b(left-k+1,k)(x), ..., b(left,k)(x),  are supplied by
+c	 bsplvb . all others are known to be zero a priori.
+	 call bsplvb ( t, k, 1, x, left, values(leftmk+1) )
+c
+	 print 610, x, (values(j),j=3,7)
+  610	 format(f7.3,5f12.7)
+c
+c	    zero out the values just computed in preparation for next
+c	 evalulation point .
+	 do 10 j=1,k
+   10	    values(leftmk+j) = 0.
+					stop
+      end
diff --git a/math/deboor/progs/prog5.f b/math/deboor/progs/prog5.f
new file mode 100644
index 00000000..b1fbb927
--- /dev/null
+++ b/math/deboor/progs/prog5.f
@@ -0,0 +1,27 @@
+chapter x. example 2. plotting the pol,s which make up a b-spline
+c  from  * a practical guide to splines *  by c. de boor
+calls  bsplvb
+c
+      integer ia,left
+      real biatx(4),t(11),values(4),x
+c		  knot sequence set here....
+      data t / 4*0.,1.,3.,4.,4*6. /
+      do 20 ia=1,40
+	 x = float(ia)*.2 - 1.
+	 do 10 left=4,7
+	    call bsplvb ( t, 4, 1, x, left, biatx )
+c
+c	    according to  bsplvb  listing,  biatx(.) now contains value
+c	    at	x  of polynomial which agrees on the interval  (t(left),
+c	    t(left+1) )  with the b-spline  b(left-4 + . ,4,t) . hence,
+c	    biatx(8-left)  now contains value of that polynomial for
+c	    b(left-4 +(8-left) ,4,t) = b(4,4,t) . since this b-spline
+c	    has support  (t(4),t(8)), it makes sense to run  left = 4,
+c	    ...,7, storing the resulting values in  values(1),...,
+c	    values(4)  for later printing.
+c
+   10	    values(left-3) = biatx(8-left)
+   20	 print 620, x, values
+  620 format(f10.1,4f20.8)
+					stop
+      end
diff --git a/math/deboor/progs/prog6.f b/math/deboor/progs/prog6.f
new file mode 100644
index 00000000..10c55c1f
--- /dev/null
+++ b/math/deboor/progs/prog6.f
@@ -0,0 +1,23 @@
+chapter x. example 3. construction and evaluation of the pp-representat-
+c		      ion of a b-spline.
+c  from  * a practical guide to splines *  by c. de boor
+calls bsplpp(bsplvb),ppvalu(interv)
+c
+      integer ia,l
+      real bcoef(7),break(5),coef(4,4),scrtch(4,4),t(11),value,x
+      real ppvalu
+c		 set knot sequence  t  and b-coeffs for  b(4,4,t)  ....
+      data t / 4*0.,1.,3.,4.,4*6. /, bcoef / 3*0.,1.,3*0. /
+c			   construct pp-representation	....
+      call bsplpp ( t, bcoef, 7, 4, scrtch, break, coef, l )
+c
+c     as a check, evaluate  b(4,4,t)  from its pp-repr. on a fine mesh.
+c	      the values should agree with (some of) those generated in
+c			 example 2 .
+      do 20 ia=1,40
+	 x = float(ia)*.2 - 1.
+	 value = ppvalu ( break, coef, l, 4, x, 0 )
+   20	 print 620, x, value
+  620 format(f10.1,f20.8)
+					stop
+      end
diff --git a/math/deboor/progs/prog7.f b/math/deboor/progs/prog7.f
new file mode 100644
index 00000000..e08ab120
--- /dev/null
+++ b/math/deboor/progs/prog7.f
@@ -0,0 +1,17 @@
+chapter x.  example 4. construction of a b-spline via  bvalue
+c  from  * a practical guide to splines *  by c. de boor
+calls  bvalue(interv)
+      integer ia
+      real bcoef(1),t(5),value,x
+      real bvalue
+c		   set knot sequence  t  and  b-coeffs for  b(1,4,t)
+      data t / 0.,1.,3.,4.,6. / , bcoef / 1. /
+c	evaluate  b(1,4,t)  on a fine mesh.  on (0,6), the values should
+c		coincide with those obtained in example 3 .
+      do 20 ia=1,40
+	 x = float(ia)*.2 - 1.
+	 value = bvalue ( t, bcoef, 1, 4, x, 0 )
+   20	 print 620, x, value
+  620 format(f10.1,f20.8)
+					stop
+      end
diff --git a/math/deboor/progs/prog8.f b/math/deboor/progs/prog8.f
new file mode 100644
index 00000000..197f56f3
--- /dev/null
+++ b/math/deboor/progs/prog8.f
@@ -0,0 +1,63 @@
+chapter xii, example 2. cubic spline interpolation with good knots
+c  from  * a practical guide to splines *  by c. de boor
+calls  cubspl, newnot
+c     parameter nmax=20
+      integer i,istep,iter,itermx,j,n,nhigh,nlow,nm1
+c     real algerp,aloger,decay,dx,errmax,c(4,nmax),g,h,pnatx
+c    *	  ,scrtch(2,nmax),step,tau(nmax),taunew(nmax)
+      real algerp,aloger,decay,dx,errmax,c(4,20),g,h,pnatx
+     *	  ,scrtch(2,20),step,tau(20),taunew(20),x
+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
+c					 knots are initially equispaced.
+	 nm1 = n-1
+	 h = 2./float(nm1)
+	 do 10 i=1,n
+   10	    tau(i) = float(i-1)*h - 1.
+	 iter = 1
+c		construct cubic spline interpolant. then, itermx  times,
+c		 determine new knots from it and find a new interpolant.
+   11	    do 15 i=1,n
+   15	       c(1,i) = g(tau(i))
+	    call cubspl ( tau, c, n, 0, 0 )
+	    if (iter .gt. itermx)	go to 19
+	    iter = iter+1
+	    call newnot(tau,c,nm1,4,taunew,nm1,scrtch)
+	    do 18 i=1,n
+   18	       tau(i) = taunew(i)
+					go to 11
+   19	 continue
+c			     estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+c			    loop over polynomial pieces of interpolant .
+	 do 30 i=1,nm1
+	    dx = (tau(i+1)-tau(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(tau(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
diff --git a/math/deboor/progs/prog9.f b/math/deboor/progs/prog9.f
new file mode 100644
index 00000000..a0dcfd60
--- /dev/null
+++ b/math/deboor/progs/prog9.f
@@ -0,0 +1,38 @@
+chapter xii, example 3. cubic spline interpolation with good knots
+c  from  * a practical guide to splines *  by c. de boor
+calls cubspl
+      integer i,irate,istep,j,n,nhigh,nlow,nm1
+      real algerp,aloger,c(4,20),decay,dx,errmax,g,h,pnatx,step,tau(20)
+     *	  ,x
+      data step, istep /20., 20/
+      g(x) = sqrt(x+1.)
+      decay = 0.
+      read 500,irate,nlow,nhigh
+  500 format(3i3)
+      print 600
+  600 format(28h  n   max.error   decay exp./)
+      do 40 n=nlow,nhigh,2
+	 nm1 = n-1
+	 h = 1./float(nm1)
+	 do 10 i=1,n
+	    tau(i) = 2.*(float(i-1)*h)**irate - 1.
+   10	    c(1,i) = g(tau(i))
+c	 construct cubic spline interpolant.
+	 call cubspl ( tau, c, n, 0, 0 )
+c	 estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+	 do 30 i=1,nm1
+	    dx = (tau(i+1)-tau(i))/step
+	    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.)
+	       x = tau(i) + h
+   30	       errmax = amax1(errmax,abs(g(x)-pnatx))
+	 aloger = alog(errmax)
+	 if (n .gt. nlow)  decay =
+     *	     (aloger - algerp)/alog(float(n)/float(n-2))
+	 algerp = aloger
+   40	 print 640,n,errmax,decay
+  640 format(i3,e12.4,f11.2)
+					stop
+      end