diff options
Diffstat (limited to 'math/deboor/l2appr.f')
| -rw-r--r-- | math/deboor/l2appr.f | 109 |
1 files changed, 109 insertions, 0 deletions
diff --git a/math/deboor/l2appr.f b/math/deboor/l2appr.f new file mode 100644 index 00000000..6f61150c --- /dev/null +++ b/math/deboor/l2appr.f @@ -0,0 +1,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 |