aboutsummaryrefslogtreecommitdiff
path: root/math/deboor/progs
diff options
context:
space:
mode:
Diffstat (limited to 'math/deboor/progs')
-rw-r--r--math/deboor/progs/prog1.f45
-rw-r--r--math/deboor/progs/prog10.f76
-rw-r--r--math/deboor/progs/prog11.f69
-rw-r--r--math/deboor/progs/prog12.f70
-rw-r--r--math/deboor/progs/prog13.f77
-rw-r--r--math/deboor/progs/prog14.f37
-rw-r--r--math/deboor/progs/prog15.f48
-rw-r--r--math/deboor/progs/prog16.f115
-rw-r--r--math/deboor/progs/prog17.f10
-rw-r--r--math/deboor/progs/prog18.f35
-rw-r--r--math/deboor/progs/prog19.f58
-rw-r--r--math/deboor/progs/prog2.f48
-rw-r--r--math/deboor/progs/prog20.f62
-rw-r--r--math/deboor/progs/prog21.f74
-rw-r--r--math/deboor/progs/prog3.f44
-rw-r--r--math/deboor/progs/prog4.f35
-rw-r--r--math/deboor/progs/prog5.f27
-rw-r--r--math/deboor/progs/prog6.f23
-rw-r--r--math/deboor/progs/prog7.f17
-rw-r--r--math/deboor/progs/prog8.f63
-rw-r--r--math/deboor/progs/prog9.f38
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
generated by cgit (git 2.34.1) at 2025-09-20 08:46:17 -0400