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 integer k,n, i,j,jj,kmax,left,leftmk,ll,mm,ntau,ntmax parameter (kmax=20,ntmax=200) real bcoef(n),diag(n),q(k,n),t(n+k), biatx(kmax),dw,gtau,tau,weight common / data / ntau, tau(ntmax),gtau(ntmax),weight(ntmax) c 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