subroutine bsplvd8 ( t, k, x, left, dbiatx, nderiv )
!  REWRITTEN  in doulbe precesition and fortran 90 by Baturin, 2015
!
!  from  * a practical guide to splines *  by c. de Boor (7 may 92)    
!   calls bsplvb
!   calculates value and deriv.s of all b-splines which do not vanish at x
!
!******  i n p u t  ******
!  t     the knot array, of length left+k (at least)
!  k     the order of the b-splines to be evaluated (power + 1)
!  x     the point at which these values are sought
!  left  an integer indicating 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.
!  a s s u m p t i o n  - - -  it is assumed that
!               t(left) .lt. t(left+1)
!        division by zero will result otherwise (in  b s p l v b ).
!        also, the output is as advertised only if
!               t(left) .le. x .le. t(left+1) .
!  nderiv   an integer indicating that values of b-splines and their
!        derivatives up to but not including the  nderiv-th  are asked
!        for. ( nderiv  is replaced internally by the integer  m h i g h
!        in  (1,k)  closest to it.)
!
!******  w o r k   a r e a  ******  <just pre-allocated outside array, with dynamic dimensions>
!  a     an array of order (k,k), to contain b-coeff.s of the derivat-
!        ives of a certain order of the  k  b-splines of interest.
!
!******  o u t p u t  ******
!  dbiatx   an array of order (k,nderiv). its entry  (i,m)  contains
!        value of  (m-1)st  derivative of  (left-k+i)-th  b-spline of
!        order  k  for knot sequence  t , i=1,...,k, m=1,...,nderiv.
!
!******  m e t h o d  ******
!  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-coeffs of the required derivatives of the b-
!  splines of interest are generated by differencing, each from the pre-
!  ceding one of lower order, and combined with the values of b-splines
!  of corresponding order in  dbiatx  to produce the desired values .
!
      implicit none
      integer, intent (in) :: k,left,nderiv
      
      real*8, intent (out) :: dbiatx(k,nderiv)
      real*8, intent (in) :: t(:),x
!  locals:
      real*8 dbiat1(k)   ! one-dimensional array, use in bsplvb call
      real*8 factor,fkp1mm,sum
      integer ideriv, mhigh
      Integer i,il,j,jlow,jp1mid,kp1,kp1mm,ldummy,m
      integer ierr
      
      mhigh = max0(min0(nderiv,k),1)
!     mhigh is usually equal to nderiv.
      kp1 = k+1
      call bsplvb8(t,kp1-mhigh,1,x,left,dbiat1)
      dbiatx(mhigh:k,mhigh) = dbiat1(1:kp1-mhigh)
      if (mhigh .eq. 1) return
!     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
!            jp1mid = jp1mid + 1
!         enddo
         ideriv = ideriv - 1
         call bsplvb8(t,kp1-ideriv,2,x,left,dbiat1)
          dbiatx(ideriv:k,ideriv) = dbiat1(1:kp1-ideriv)
      enddo
!
!     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-
!     responding derivatives of b-splines in subsequent columns, gene-
!     rate their b-repr. by differencing, then evaluate at  x.
!
!  Esencialy, this block is initiation of a() to indentity matrix
      if (.not.allocated(a)) then
          allocate (a(k,k), stat=ierr)
          if (ierr.ne.0) stop 'Allocation of <a> error'
      elseif (size(a,dim=1).ne.k) then
          deallocate(a)
          allocate (a(k,k), stat=ierr)
          if (ierr.ne.0) stop 'Allocation of <a> error'
      end if

       a = 0.d0                           ! Initialize the array.
       forall(j = 1:k) a(j,j) = 1.d0      ! Set the diagonal.
!
!      jlow = 1
!      do i=1,k
!         do j=jlow,k
!         a(j,i) = 0.d0
!         enddo
!         jlow = i
!      a(i,i) = 1.
!      enddo
!
!     at this point, a(.,j) contains the b-coeffs for the j-th of the
!     k  b-splines of interest here.
!
      do m=2,mhigh
         kp1mm = kp1 - m
         fkp1mm = float(kp1mm)
         il = left
         i = k
!
!        for j=1,...,k, construct b-coeffs 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 .lt. j  is used.
         do ldummy=1,kp1mm
            factor = fkp1mm/(t(il+kp1mm) - t(il))
!           the assumption that t(left).lt.t(left+1) makes denominator
!           in  factor  nonzero.
            do j=1,i
               a(i,j) = (a(i,j) - a(i-1,j))*factor
            enddo
            il = il - 1
         i = i - 1
         enddo
!
!        for i=1,...,k, combine b-coeffs 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 .lt. i .
         do i=1,k
            sum = 0.d0
            jlow = max0(i,m)
            do j=jlow,k
              sum = a(j,i)*dbiatx(j,m) + sum
            enddo
         dbiatx(i,m) = sum
         enddo
         enddo
      return
      end subroutine