l2appr.f

Go to the documentation of this file.

      subroutine l2appr ( t, n, k, q, diag, bcoef )

c*********************************************************************72
c
cc L2APPR constructs a weighted L2 spline approximation to given data.
c
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
      implicit none

      integer k,n,   i,j,jj,kmax,left,leftmk,ll,mm,ntau,ntmax
      parameter(kmax=20,ntmax=200)
      double precision bcoef(n),diag(n),q(k,n),t(n+k),
     &  biatx(kmax),dw,gtau,tau,weight
      common / i4data / ntau
      common / r8data / tau(ntmax),gtau(ntmax),weight(ntmax)
c
      do 7 j=1,n
         bcoef(j) = 0.0d+00
         do 7 i=1,k
    7       q(i,j) = 0.0d+00
      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