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
|
subroutine l2appr ( t, n, k, q, diag, bcoef )
c from * a practical guide to splines * by c. de boor
c to be called in main program l 2 m a i n .
calls subprograms bsplvb, bchfac/slv
c
constructs the (weighted discrete) l2-approximation by splines of order
c k with knot sequence t(1), ..., t(n+k) to given data points
c ( tau(i), gtau(i) ), i=1,...,ntau. the b-spline coefficients
c b c o e f of the approximating spline are determined from the
c normal equations using cholesky's method.
c
c****** i n p u t ******
c t(1), ..., t(n+k) the knot sequence
c n.....the dimension of the space of splines of order k with knots t.
c k.....the order
c
c w a r n i n g . . . the restriction k .le. kmax (= 20) is impo-
c sed by the arbitrary dimension statement for biatx below, but
c is n o w h e r e c h e c k e d for.
c
c****** w o r k a r r a y s ******
c q....a work array of size (at least) k*n. its first k rows are used
c for the k lower diagonals of the gramian matrix c .
c diag.....a work array of length n used in bchfac .
c
c****** i n p u t via c o m m o n /data/ ******
c ntau.....number of data points
c (tau(i),gtau(i)), i=1,...,ntau are the ntau data points to be
c fitted .
c weight(i), i=1,...,ntau are the corresponding weights .
c
c****** o u t p u t ******
c bcoef(1), ..., bcoef(n) the b-spline coeffs. of the l2-appr.
c
c****** m e t h o d ******
c the b-spline coefficients of the l2-appr. are determined as the sol-
c ution of the normal equations
c sum ( (b(i),b(j))*bcoef(j) : j=1,...,n) = (b(i),g),
c i = 1, ..., n .
c here, b(i) denotes the i-th b-spline, g denotes the function to
c be approximated, and the i n n e r p r o d u c t of two funct-
c ions f and g is given by
c (f,g) := sum ( f(tau(i))*g(tau(i))*weight(i) : i=1,...,ntau) .
c the arrays t a u and w e i g h t are given in common block
c d a t a , as is the array g t a u containing the sequence
c g(tau(i)), i=1,...,ntau.
c the relevant function values of the b-splines b(i), i=1,...,n, are
c supplied by the subprogram b s p l v b .
c the coeff.matrix c , with
c c(i,j) := (b(i), b(j)), i,j=1,...,n,
c of the normal equations is symmetric and (2*k-1)-banded, therefore
c can be specified by giving its k bands at or below the diagonal. for
c i=1,...,n, we store
c (b(i),b(j)) = c(i,j) in q(i-j+1,j), j=i,...,min0(i+k-1,n)
c and the right side
c (b(i), g ) in bcoef(i) .
c since b-spline values are most efficiently generated by finding sim-
c ultaneously the value of e v e r y nonzero b-spline at one point,
c the entries of c (i.e., of q ), are generated by computing, for
c each ll, all the terms involving tau(ll) simultaneously and adding
c them to all relevant entries.
c parameter kmax=20,ntmax=200
integer k,n, i,j,jj,left,leftmk,ll,mm,ntau
c real bcoef(n),diag(n),q(k,n),t(1), biatx(kmax),dw,gtau,tau,weight
real bcoef(n),diag(n),q(k,n),t(1), biatx( 20),dw,gtau,tau,weight
c dimension t(n+k)
current fortran standard makes it impossible to specify the exact dimen-
c sion of t without the introduction of additional otherwise super-
c fluous arguments.
c common / data / ntau, tau(ntmax),gtau(ntmax),weight(ntmax)
common / data / ntau, tau( 200),gtau( 200),weight( 200)
do 7 j=1,n
bcoef(j) = 0.
do 7 i=1,k
7 q(i,j) = 0.
left = k
leftmk = 0
do 20 ll=1,ntau
c locate l e f t s.t. tau(ll) in (t(left),t(left+1))
10 if (left .eq. n) go to 15
if (tau(ll) .lt. t(left+1)) go to 15
left = left+1
leftmk = leftmk + 1
go to 10
15 call bsplvb ( t, k, 1, tau(ll), left, biatx )
c biatx(mm) contains the value of b(left-k+mm) at tau(ll).
c hence, with dw := biatx(mm)*weight(ll), the number dw*gtau(ll)
c is a summand in the inner product
c (b(left-k+mm), g) which goes into bcoef(left-k+mm)
c and the number biatx(jj)*dw is a summand in the inner product
c (b(left-k+jj), b(left-k+mm)), into q(jj-mm+1,left-k+mm)
c since (left-k+jj) - (left-k+mm) + 1 = jj - mm + 1 .
do 20 mm=1,k
dw = biatx(mm)*weight(ll)
j = leftmk + mm
bcoef(j) = dw*gtau(ll) + bcoef(j)
i = 1
do 20 jj=mm,k
q(i,j) = biatx(jj)*dw + q(i,j)
20 i = i + 1
c
c construct cholesky factorization for c in q , then use
c it to solve the normal equations
c c*x = bcoef
c for x , and store x in bcoef .
call bchfac ( q, k, n, diag )
call bchslv ( q, k, n, bcoef )
return
end
|
