1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
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
|
