bspp2d Subroutine

subroutine bspp2d ( t, bcoef, n, k, m, scrtch, break, coef, l )

  !*****************************************************************************80
  !
  !! BSPP2D converts from B-spline to piecewise polynomial representation.
  !
  !  Discussion:
  !
  !    The B-spline representation
  ! 
  !      T, BCOEF(.,J), N, K 
  !
  !    is converted to its piecewise polynomial representation 
  !
  !      BREAK, COEF(J,.,.), L, K, J=1, ..., M.
  !
  !    This is an extended version of BSPLPP for use with tensor products.
  !
  !    For each breakpoint interval, the K relevant B-spline 
  !    coefficients of the spline are found and then differenced 
  !    repeatedly to get the B-spline coefficients of all the 
  !    derivatives of the spline on that interval. 
  !
  !    The spline and its first K-1 derivatives are then evaluated 
  !    at the left endpoint of that interval, using BSPLVB 
  !    repeatedly to obtain the values of all B-splines of the 
  !    appropriate order at that point.
  !
  !  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, real ( kind = 8 ) BCOEF(N,M).  For each J, B(*,J) is the 
  !    B-spline coefficient sequence, of length N.
  ! 
  !    Input, integer ( kind = 4 ) N, the length of BCOEF.
  ! 
  !    Input, integer ( kind = 4 ) K, the order of the spline.
  ! 
  !    Input, integer ( kind = 4 ) M, the number of data sets.
  ! 
  !    Work array, real ( kind = 8 ) SCRTCH(K,K,M).
  ! 
  !    Output, real ( kind = 8 ) BREAK(L+1), the breakpoint sequence
  !    containing the distinct points in the sequence T(K),...,T(N+1)
  ! 
  !    Output, real ( kind = 8 ) COEF(M,K,N), with COEF(MM,I,J) = the (I-1)st
  !    derivative of the MM-th spline at BREAK(J) from the right, MM=1, ..., M.
  ! 
  !    Output, integer ( kind = 4 ) L, the number of polynomial pieces which 
  !    make up the spline in the interval (T(K), T(N+1)).
  !
  implicit none

  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) bcoef(n,m)
  real ( kind = 8 ) biatx(k)
  real ( kind = 8 ) break(*)
  real ( kind = 8 ) coef(m,k,*)
  real ( kind = 8 ) diff
  real ( kind = 8 ) fkmj
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jp1
  integer ( kind = 4 ) kmj
  integer ( kind = 4 ) l
  integer ( kind = 4 ) left
  integer ( kind = 4 ) lsofar
  integer ( kind = 4 ) mm
  real ( kind = 8 ) scrtch(k,k,m)
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t(n+k)

  lsofar = 0
  break(1) = t(k)

  do left = k, n
     !
     !  Find the next nontrivial knot interval.
     !
     if ( t(left+1) == t(left) ) then
        cycle
     end if

     lsofar = lsofar+1
     break(lsofar+1) = t(left+1)

     if ( k <= 1 ) then

        coef(1:m,1,lsofar) = bcoef(left,1:m)

        cycle

     end if
     !
     !  Store the K B-spline coefficients relevant to current knot interval
     !  in SCRTCH(.,1).
     !
     do i = 1, k
        scrtch(i,1,1:m) = bcoef(left-k+i,1:m)
     end do
     !
     !  For J = 1,...,K-1, compute the ( K - J ) B-spline coefficients relevant to
     !  current knot interval for the J-th derivative by differencing
     !  those for the (J-1)st derivative, and store in SCRTCH(.,J+1).
     !
     do jp1 = 2, k

        j = jp1 - 1
        kmj = k - j
        fkmj = real ( k - j, kind = 8 )

        do i = 1, k - j

           diff = ( t(left+i) - t(left+i-kmj) ) / fkmj

           if ( 0.0D+00 < diff ) then

              scrtch(i,jp1,1:m) = ( scrtch(i+1,j,1:m) - scrtch(i,j,1:m) ) / diff

           end if

        end do

     end do
     !
     !  For  J = 0, ..., K-1, find the values at T(LEFT) of the J+1
     !  B-splines of order J+1 whose support contains the current
     !  knot interval from those of order J (in  BIATX ), then combine
     !  with the B-spline coefficients (in SCRTCH(.,K-J) ) found earlier
     !  to compute the (K-J-1)st derivative at T(LEFT) of the given spline.
     !
     call bsplvb ( t, 1, 1, t(left), left, biatx )

     coef(1:m,k,lsofar) = scrtch(1,k,1:m)

     do jp1 = 2, k

        call bsplvb ( t, jp1, 2, t(left), left, biatx )
        kmj = k + 1 - jp1

        do mm = 1, m

           sum1 = 0.0D+00
           do i = 1, jp1
              sum1 = sum1 + biatx(i) * scrtch(i,kmj,mm)
           end do

           coef(mm,kmj,lsofar) = sum1

        end do

     end do

  end do

  l = lsofar

  return
end subroutine bspp2d