s_bsplines_non_uniform__eval_basis_and_n_derivs Subroutine

  SLL_PURE subroutine s_bsplines_non_uniform__eval_basis_and_n_derivs( self, x, n, derivs, jmin )
    class(sll_t_bsplines_non_uniform), intent(in   ) :: self
    real(wp)                         , intent(in   ) :: x
    integer                          , intent(in   ) :: n
    real(wp)                         , intent(  out) :: derivs(0:,0:)
    integer                          , intent(  out) :: jmin

    integer  :: icell

    integer  :: j, r
    real(wp) :: temp
    real(wp) :: saved

    integer  :: k, s1, s2, rk, pk, j1, j2
    real(wp) :: d

    ! GFortran: to allocate on stack use -fstack-arrays
    real(wp) :: left (1:self%degree)
    real(wp) :: right(1:self%degree)
    real(wp) :: ndu  (0:self%degree,0:self%degree)
    real(wp) :: a    (0:1          ,0:self%degree)

    ! Check on inputs
    SLL_ASSERT( x > self%xmin - 1.0d-14 )
    SLL_ASSERT( x < self%xmax + 1.0d-14 )
    SLL_ASSERT( n >= 0 )
    SLL_ASSERT( n <= self%degree )
    SLL_ASSERT( size(derivs,1) == 1+n           )
    SLL_ASSERT( size(derivs,2) == 1+self%degree )

    ! 1. Compute cell index 'icell' and x_offset
!    icell = f_bsplines_non_uniform__find_cell( self, x )
    icell = self % find_cell( x )
    SLL_ASSERT( icell >= 1               )
    SLL_ASSERT( icell <= self%ncells     )
    SLL_ASSERT( self%knots(icell) <= x   )
    SLL_ASSERT( x <= self%knots(icell+1) )

    ! 2. Compute index range of B-splines with support over cell 'icell'
    jmin = icell - self%offset

    ! 3. Compute nonzero basis functions and knot differences for splines
    !    up to degree (degree-1) which are needed to compute derivative
    !    Algorithm  A2.3 of NURBS book
    !
    !    21.08.2017: save inverse of knot differences to avoid unnecessary divisions
    !                [Yaman Güçlü, Edoardo Zoni]
    associate( degree => self%degree )

      ndu(0,0) = 1.0_wp
      do j = 1, degree
         left(j)  = x - self%knots(icell+1-j)
         right(j) = self%knots(icell+j) - x
         saved    = 0.0_wp
         do r = 0, j-1
            ! compute inverse of knot differences and save them into lower triangular part of ndu
            ndu(j,r) = 1.0_wp / (right(r+1) + left(j-r))
            ! compute basis functions and save them into upper triangular part of ndu
            temp     = ndu(r,j-1) * ndu(j,r)
            ndu(r,j) = saved + right(r+1) * temp
            saved    = left(j-r) * temp
         end do
         ndu(j,j) = saved
      end do
      derivs(0,:) = ndu(:,degree)

      do r = 0, degree
         s1 = 0
         s2 = 1
         a(0,0) = 1.0_wp
         do k = 1, n
            d  = 0.0_wp
            rk = r-k
            pk = degree-k
            if (r >= k) then
               a(s2,0) = a(s1,0) * ndu(pk+1,rk)
               d = a(s2,0) * ndu(rk,pk)
            end if
            if (rk > -1) then
               j1 = 1
            else
               j1 = -rk
            end if
            if (r-1 <= pk) then
               j2 = k-1
            else
               j2 = degree-r
            end if
            do j = j1, j2
               a(s2,j) = (a(s1,j) - a(s1,j-1)) * ndu(pk+1,rk+j)
               d = d + a(s2,j) * ndu(rk+j,pk)
            end do
            if (r <= pk) then
               a(s2,k) = - a(s1,k-1) * ndu(pk+1,r)
               d = d + a(s2,k) * ndu(r,pk)
            end if
            derivs(k,r) = d
            j  = s1
            s1 = s2
            s2 = j
         end do
      end do
      r = degree
      do k = 1, n
         derivs(k,:) = derivs(k,:) * r
         r = r * (degree-k)
      end do

    end associate

  end subroutine s_bsplines_non_uniform__eval_basis_and_n_derivs