!>
!> @file bsplvd.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 bsplvd ( t, k, x, left, a, dbiatx, nderiv )

!*************************************************************************
!
!! BSPLVD calculates the nonvanishing B-splines and derivatives at X.
!
!  Discussion:
!
!    Values at X of all the relevant B-splines of order K, K-1,..., K+1-NDERIV
!    are generated via BSPLVB and stored temporarily in DBIATX.
!
!    Then, the B-spline coefficients of the required derivatives
!    of the B-splines of interest are generated by differencing,
!    each from the preceding one of lower order, and combined with
!    the values of B-splines of corresponding order in DBIATX
!    to produce the desired values.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) T(LEFT+K), the knot sequence.  It is assumed that
!    T(LEFT) < T(LEFT+1).  Also, the output is correct only if
!    T(LEFT) <= X <= T(LEFT+1) .
!
!    Input, integer K, the order of the B-splines to be evaluated.
!
!    Input, real ( kind = 8 ) X, the point at which these values are sought.
!
!    Input, integer LEFT, indicates the left endpoint of the interval of
!    interest.  The K B-splines whose support contains the interval
!    (T(LEFT), T(LEFT+1)) are to be considered.
!
!    Workspace, real ( kind = 8 ) A(K,K).
!
!    Output, real ( kind = 8 ) DBIATX(K,NDERIV).  DBIATX(I,M) contains
!    the value of the (M-1)st derivative of the (LEFT-K+I)-th B-spline
!    of order K for knot sequence T, I=M,...,K, M=1,...,NDERIV.
!
!    Input, integer NDERIV, indicates that values of
!    B-splines and their derivatives up to but not
!    including the NDERIV-th are asked for.
!
  implicit none

  integer k
  integer left
  integer nderiv

  real ( kind = 8 ) a(k,k)
  real ( kind = 8 ) dbiatx(k,nderiv)
  real ( kind = 8 ) factor
  real ( kind = 8 ) fkp1mm
  integer i
  integer ideriv
  integer il
  integer j
  integer jlow
  integer jp1mid
  integer ldummy
  integer m
  integer mhigh
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t(left+k)
  real ( kind = 8 ) x

  mhigh = max ( min ( nderiv, k ), 1 )
!
!  MHIGH is usually equal to nderiv.
!
  call bsplvb ( t, k+1-mhigh, 1, x, left, dbiatx )

  if ( mhigh == 1 ) then
    return
  end if
!
!  The first column of DBIATX always contains the B-spline values
!  for the current order.  These are stored in column K+1-current
!  order  before BSPLVB is called to put values for the next
!  higher order on top of it.
!
  ideriv = mhigh
  do m = 2, mhigh
    jp1mid = 1
    do j = ideriv, k
      dbiatx(j,ideriv) = dbiatx(jp1mid,1)
      jp1mid = jp1mid+1
    end do
    ideriv = ideriv-1
    call bsplvb(t,k+1-ideriv,2,x,left,dbiatx)
  end do
!
!  At this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j) for
!  i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
!  first column of  dbiatx  is already in final form. to obtain cor-
!  ???  LOST A LINE ???
!  rate their b-repr. by differencing, then evaluate at  x.
!
  jlow = 1
  do i = 1, k
    do j = jlow,k
      a(j,i) = 0.0D+00
    end do
    jlow = i
    a(i,i) = 1.0D+00
  end do
!
!  At this point, a(.,j) contains the b-coefficients for the J-th of the
!  k  b-splines of interest here.
!
  do m = 2, mhigh

    fkp1mm = real ( k + 1 - m, kind = 8 )
    il = left
    i = k
!
!  For j=1,...,k, construct b-coefficients of  (m-1)st  derivative of
!  b-splines from those for preceding derivative by differencing
!  and store again in  a(.,j) .  The fact that  a(i,j)=0  for
!  i < j  is used.
!
    do ldummy = 1, k+1-m

      factor = fkp1mm/(t(il+k+1-m)-t(il))
!
!  The assumption that t(left) < t(left+1) makes denominator
!  in  factor  nonzero.
!
      do j = 1, i
        a(i,j) = (a(i,j)-a(i-1,j))*factor
      end do

      il = il-1
      i = i-1

    end do
!
!  For i=1,...,k, combine b-coefficients a(.,i) with B-spline values
!  stored in dbiatx(.,m) to get value of  (m-1)st  derivative of
!  i-th b-spline (of interest here) at  x , and store in
!  dbiatx(i,m). storage of this value over the value of a b-spline
!  of order m there is safe since the remaining b-spline derivat-
!  ives of the same order do not use this value due to the fact
!  that  a(j,i)=0  for j < i.
!
    do i = 1, k

      sum1 = 0.0D+00
      jlow = max(i,m)
      do j = jlow,k
        sum1 = sum1 + a(j,i) * dbiatx(j,m)
      end do

      dbiatx(i,m) = sum1
    end do

  end do

  return
end