!>
!> @file bspp2d.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine bspp2d ( t, bcoef, n, k, m, scrtch, break, coef, l )

!*************************************************************************
!
!! 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.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  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 N, the length of BCOEF.
!
!    Input, integer K, the order of the spline.
!
!    Input, integer 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 L, the number of polynomial pieces which make up the
!    spline in the interval (T(K), T(N+1)).
!
  implicit none

  integer k
  integer m
  integer 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 i
  integer j
  integer jp1
  integer kmj
  integer l
  integer left
  integer lsofar
  integer 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

      do mm = 1, m
        coef(mm,1,lsofar) = bcoef(left,mm)
      end do

      cycle

    end if
!
!  Store the K b-spline coefficients relevant to current knot interval
!  in  scrtch(.,1) .
!
    do i = 1, k
      do mm = 1, m
        scrtch(i,1,mm) = bcoef(left-k+i,mm)
      end do
    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

          do mm = 1, m
            scrtch(i,jp1,mm)=(scrtch(i+1,j,mm)-scrtch(i,j,mm))/diff
          end do

        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 comb-
!  ine 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 )

    do mm = 1, m
      coef(mm,k,lsofar) = scrtch(1,k,mm)
    end do

    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