c main program for least-squares approximation by splines
c from * a practical guide to splines * by c. de Boor (7 may 92)
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 = 'on' gives printout of b-spline coeffs. of approxim.
c p r p c o = 'on' gives printout of pp repr. of approximation.
c p r f u n = 'on' 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
integer i,icount,ii,j,k,l,lbegin,lnew,ll,lpkmax,ltkmax,n,nt,ntau
* ,ntimes,ntmax,on,prbco,prfun,prpco
parameter (lpkmax=100,ntmax=200,ltkmax=2000)
real addbrk,bcoef(lpkmax),break,coef,gtau,q(ltkmax),scrtch(ntmax)
* ,t(ntmax),tau,totalw,weight
common / data / ntau, tau(ntmax),gtau(ntmax),weight(ntmax),totalw
C real addbrk,bcoef(100),break,coef,gtau,q(2000),scrtch(200)
C * ,t(200),tau,totalw,weight
C 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.
common /approx/ break(lpkmax),coef(ltkmax),l,k
C common /approx/ break(100),coef(2000),l,k
c common /approx/ also occurs in setdat and l2err.
data on /'ON'/
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(' 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 = 0
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(//' b-spline coefficients'/(4e20.10))
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(//' approximation by splines of order',i3,' on ',
* i3,' intervals. breakpoints -'/(4e20.10))
if (prpco .ne. on) go to 15
print 611
611 format(/' 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,4e20.10/(11x,4e20.10))
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