l2appr Subroutine

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

  !*****************************************************************************80
  !
  !! L2APPR constructs a weighted L2 spline approximation to given data.
  !
  !  Discussion:
  !
  !    The routine constructs the weighted discrete L2-approximation by 
  !    splines of order K with knot sequence T(1:N+K) to 
  !    given data points ( TAU(1:NTAU), GTAU(1:NTAU) ).  
  !
  !    The B-spline coefficients BCOEF of the approximating spline are 
  !    determined from the normal equations using Cholesky's method.
  !
  !  Method:
  !
  !    The B-spline coefficients of the L2-approximation are determined as the 
  !    solution of the normal equations, for 1 <= I <= N:
  !
  !      sum ( 1 <= J <= N ) ( B(I), B(J) ) * BCOEF(J) = ( B(I), G ).
  !
  !    Here, B(I) denotes the I-th B-spline, G denotes the function to
  !    be approximated, and the inner product of two functions F and G 
  !    is given by
  !
  !      ( F, G ) = sum ( 1 <= I <= NTAU ) WEIGHT(I) * F(TAU(I)) * G(TAU(I)).
  !
  !    The arrays TAU and WEIGHT are given in common block DATA, as is the 
  !    array GTAU(1:NTAU) = G(TAU(1:NTAU)).
  !
  !    The values of the B-splines B(1:N) are supplied by BSPLVB.
  !
  !    The coefficient matrix C, with
  !
  !       C(I,J) = ( B(I), B(J) )
  !
  !    of the normal equations is symmetric and (2*K-1)-banded, therefore
  !    can be specified by giving its K bands at or below the diagonal. 
  !
  !    For I = 1:N and J = I:min(I+K-1,N), we store
  !
  !      ( B(I), B(J) ) = C(I,J) 
  !
  !    in
  !
  !      Q(I-J+1,J), 
  !
  !    and the right hand side
  !
  !      ( B(I), G )  
  !
  !    in
  !
  !      BCOEF(I).
  !
  !    Since B-spline values are most efficiently generated by finding
  !    simultaneously the value of every nonzero B-spline at one point,
  !    the entries of C (that is, of Q), are generated by computing, for
  !    each LL, all the terms involving TAU(LL) simultaneously and adding
  !    them to all relevant entries.
  !
  !  Modified:
  !
  !    14 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) T(N+K), the knot sequence.
  !
  !    Input, integer ( kind = 4 ) N, the dimension of the space of splines 
  !    of order K with knots T.
  !
  !    Input, integer ( kind = 4 ) K, the order of the splines.
  !
  !    Workspace, real ( kind = 8 ) Q(K,N), used to store the K lower 
  !    diagonals of the Gramian matrix C.
  !
  !    Workspace, real ( kind = 8 ) DIAG(N), used in BCHFAC.
  !
  !    Output, real ( kind = 8 ) BCOEF(N), the B-spline coefficients of
  !    the L2 approximation to the data.
  !
  implicit none

  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  integer ( kind = 4 ), parameter :: ntmax = 200

  real ( kind = 8 ) bcoef(n)
  real ( kind = 8 ) biatx(k)
  real ( kind = 8 ) diag(n)
  real ( kind = 8 ) dw
  real ( kind = 8 ) gtau
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jj
  integer ( kind = 4 ) left
  integer ( kind = 4 ) leftmk
  integer ( kind = 4 ) ll
  integer ( kind = 4 ) mm
  integer ( kind = 4 ) ntau
  real ( kind = 8 ) q(k,n)
  real ( kind = 8 ) t(n+k)
  real ( kind = 8 ) tau
  real ( kind = 8 ) totalw
  real ( kind = 8 ) weight

  save / i4data /
  save / r8data /

  common / i4data / ntau
  common / r8data / tau(ntmax), gtau(ntmax), weight(ntmax), totalw

  bcoef(1:n) = 0.0D+00
  q(1:k,1:n) = 0.0D+00

  left = k
  leftmk = 0

  do ll = 1, ntau
     !
     !  Locate LEFT such that TAU(LL) is in ( T(LEFT), T(LEFT+1) ).
     !
     do

        if ( left == n ) then
           exit
        end if

        if ( tau(ll) < t(left+1) ) then
           exit
        end if

        left = left + 1
        leftmk = leftmk + 1

     end do

     call bsplvb ( t, k, 1, tau(ll), left, biatx )
     !
     !  BIATX(MM) contains the value of B(LEFT-K+MM) at TAU(LL).
     !
     !  Hence, with DW = BIATX(MM) * WEIGHT(LL), the number DW * GTAU(LL)
     !  is a summand in the inner product
     !
     !    ( B(LEFT-K+MM), G)
     !
     !  which goes into  BCOEF(LEFT-K+MM)
     !  and the number BIATX(JJ)*DW is a summand in the inner product
     !    (B(LEFT-K+JJ), B(LEFT-K+MM)), into  Q(JJ-MM+1,LEFT-K+MM)
     !  since  (LEFT-K+JJ)-(LEFT-K+MM)+1 = JJ - MM + 1.
     !
     do mm = 1, k

        dw = biatx(mm) * weight(ll)
        j = leftmk + mm
        bcoef(j) = dw * gtau(ll) + bcoef(j)
        i = 1

        do jj = mm, k
           q(i,j) = biatx(jj) * dw + q(i,j)
           i = i + 1
        end do

     end do

  end do
  !
  !  Construct the Cholesky factorization for C in Q, then 
  !  use it to solve the normal equations
  !
  !    C * X = BCOEF
  !
  !  for X, and store X in BCOEF.
  !
  call bchfac ( q, k, n, diag )

  call bchslv ( q, k, n, bcoef )

  return
end subroutine l2appr