! Copyright (c) INRIA
! License: CECILL-B
!
!> @ingroup splines
!> @brief   Interpolator for 1D splines of arbitrary degree,
!>          on uniform and non-uniform grids
!> @author  Yaman Güçlü  - IPP Garching
!> @author  Edoardo Zoni - IPP Garching

module sll_m_spline_interpolator_1d
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_errors.h"

  use sll_m_working_precision, only: f64

  use sll_m_boundary_condition_descriptors, only: &
    sll_p_periodic, &
    sll_p_hermite , &
    sll_p_greville

  use sll_m_bsplines_base, only: sll_c_bsplines
  use sll_m_spline_1d, only: sll_t_spline_1d

  use sll_m_spline_matrix, only: &
    sll_c_spline_matrix, &
    sll_s_spline_matrix_new

  implicit none

  public :: &
    sll_t_spline_interpolator_1d, &
    sll_s_spline_1d_compute_num_cells

  private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Working precision
  integer, parameter :: wp = f64

  !> Allowed boundary conditions
  integer, parameter :: &
    allowed_bcs(1:3) = [sll_p_periodic, sll_p_hermite, sll_p_greville]

  !-----------------------------------------------------------------------------
  !> 1D spline interpolator
  !-----------------------------------------------------------------------------
  type :: sll_t_spline_interpolator_1d

    ! Private attributes
    class(sll_c_bsplines),          pointer, private :: bspl => null()
    integer                                , private ::  bc_xmin
    integer                                , private ::  bc_xmax
    integer                                , private :: nbc_xmin
    integer                                , private :: nbc_xmax
    integer                                , private :: odd
    real(wp)                               , private :: dx
    real(wp),                   allocatable, private :: tau(:)
    class(sll_c_spline_matrix), allocatable, private :: matrix

  contains

    procedure :: init                => s_spline_interpolator_1d__init
    procedure :: free                => s_spline_interpolator_1d__free
    procedure :: get_interp_points   => s_spline_interpolator_1d__get_interp_points
    procedure :: compute_interpolant => s_spline_interpolator_1d__compute_interpolant

  end type sll_t_spline_interpolator_1d

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !-----------------------------------------------------------------------------
  !> @brief      Calculate number of cells from number of interpolation points
  !> @details    Important for parallelization: for given spline degree and BCs,
  !>             calculate the number of grid cells that yields the desired
  !>             number of interpolation points
  !>
  !> @param[in]  degree   spline degree
  !> @param[in]  bc_xmin  boundary condition type at left  boundary (x=xmin)
  !> @param[in]  bc_xmax  boundary condition type at right boundary (x=xmax)
  !> @param[in]  nipts    desired number of interpolation points
  !> @param[out] ncells   calculated number of cells in domain
  !-----------------------------------------------------------------------------
  subroutine sll_s_spline_1d_compute_num_cells( &
      degree , &
      bc_xmin, &
      bc_xmax, &
      nipts  , &
      ncells )

    integer, intent(in   ) :: degree
    integer, intent(in   ) :: bc_xmin
    integer, intent(in   ) :: bc_xmax
    integer, intent(in   ) :: nipts
    integer, intent(  out) :: ncells

    ! Sanity checks
    SLL_ASSERT( degree > 0  )
    SLL_ASSERT( any( bc_xmin == allowed_bcs ) )
    SLL_ASSERT( any( bc_xmax == allowed_bcs ) )

    if (any( [bc_xmin,bc_xmax]==sll_p_periodic ) .and. bc_xmin /= bc_xmax) then
      SLL_ERROR( "sll_s_spline_1d_compute_num_cells", "Incompatible BCs" )
    end if

    if (bc_xmin == sll_p_periodic) then
      ncells = nipts
    else
      associate( nbc_xmin => merge( degree/2, 0, bc_xmin == sll_p_hermite ), &
                 nbc_xmax => merge( degree/2, 0, bc_xmax == sll_p_hermite ) )
        ncells = nipts + nbc_xmin + nbc_xmax - degree
      end associate
    end if

  end subroutine sll_s_spline_1d_compute_num_cells

  !-----------------------------------------------------------------------------
  !> @brief      Initialize a 1D spline interpolator object
  !> @param[out] self     1D spline interpolator
  !> @param[in]  bspl     B-splines (basis)
  !> @param[in]  bc_xmin  boundary condition at xmin
  !> @param[in]  bc_xmax  boundary condition at xmax
  !-----------------------------------------------------------------------------
  subroutine s_spline_interpolator_1d__init( self, bspl, bc_xmin, bc_xmax )

    class(sll_t_spline_interpolator_1d), intent(  out) :: self
    class(sll_c_bsplines),       target, intent(in   ) :: bspl
    integer                            , intent(in   ) :: bc_xmin
    integer                            , intent(in   ) :: bc_xmax

    integer :: kl, ku
    character(len=32) :: matrix_type

    ! Sanity checks
    SLL_ASSERT( any( bc_xmin == allowed_bcs ) )
    SLL_ASSERT( any( bc_xmax == allowed_bcs ) )
    if (bspl % periodic) then
      SLL_ASSERT( bc_xmin == sll_p_periodic )
      SLL_ASSERT( bc_xmax == sll_p_periodic )
    end if

    ! Save pointer to B-splines
    ! (later needed to verify 1D spline input to 'compute_interpolant')
    self % bspl => bspl

    ! Save other data
    self % bc_xmin  = bc_xmin
    self % bc_xmax  = bc_xmax
    self % nbc_xmin = merge ( bspl%degree/2, 0, bc_xmin == sll_p_hermite )
    self % nbc_xmax = merge ( bspl%degree/2, 0, bc_xmax == sll_p_hermite )
    self % odd      = modulo( bspl%degree  , 2 )

    ! Save average cell size for normalization of derivatives
    self%dx = (bspl%xmax - bspl%xmin) / bspl%ncells

    ! Compute interpolation points and number of diagonals in linear system
    if (bspl % uniform) then
      call s_compute_interpolation_points_uniform( self, self%tau )
      call s_compute_num_diags_uniform( self, kl, ku )
    else
      call s_compute_interpolation_points_non_uniform( self, self%tau )
      call s_compute_num_diags_non_uniform( self, kl, ku )
    end if

    ! Special case: linear spline
    ! No need for matrix assembly
    if (self % bspl % degree == 1) return

    ! Determine matrix storage type and allocate matrix
    associate( nbasis => self % bspl % nbasis )

      if (bc_xmin == sll_p_periodic) then
        if (kl+1+ku >= nbasis) then
          matrix_type = "dense"
        else
          matrix_type = "periodic_banded"
        end if
      else
        matrix_type = "banded"
      end if

      call sll_s_spline_matrix_new( self%matrix, matrix_type, nbasis, kl, ku )

    end associate ! nbasis

    ! Fill in entries A_ij of linear system A*x = b for interpolation
    call s_build_system( self, self%matrix )

    ! Factorize matrix A to speed-up solution for any given b
    call self % matrix % factorize()

  end subroutine s_spline_interpolator_1d__init

  !-----------------------------------------------------------------------------
  !> @brief        Destroy local objects and free allocated memory
  !> @param[inout] self  1D spline interpolator
  !-----------------------------------------------------------------------------
  subroutine s_spline_interpolator_1d__free( self )

    class(sll_t_spline_interpolator_1d), intent(inout) :: self

    nullify   ( self % bspl )
    deallocate( self % tau  )

    if (allocated( self % matrix )) then
      call self % matrix % free()
      deallocate( self % matrix )
    end if

  end subroutine s_spline_interpolator_1d__free

  !-----------------------------------------------------------------------------
  !> @brief      Get coordinates of interpolation points (1D grid)
  !> @param[in]  self  1D spline interpolator
  !> @param[out] tau   x coordinates of interpolation points
  !-----------------------------------------------------------------------------
  subroutine s_spline_interpolator_1d__get_interp_points( self, tau )

    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    real(wp),               allocatable, intent(  out) :: tau(:)

    SLL_ASSERT( allocated( self%tau ) )
    allocate( tau(size(self%tau)), source=self%tau )

  end subroutine s_spline_interpolator_1d__get_interp_points

  !-----------------------------------------------------------------------------
  !> @brief        Compute interpolating 1D spline
  !> @details      Compute coefficients of 1D spline that interpolates function
  !>               values on grid. If Hermite BCs are used, function derivatives
  !>               at appropriate boundary are also needed.
  !>
  !> @param[in]    self        1D spline interpolator
  !> @param[inout] spline      1D spline
  !> @param[in]    gtau        function values of interpolation points
  !> @param[in]    derivs_xmin (optional) array with boundary conditions at xmin
  !> @param[in]    derivs_xmax (optional) array with boundary conditions at xmax
  !-----------------------------------------------------------------------------
  subroutine s_spline_interpolator_1d__compute_interpolant( self, &
      bcoef, gtau, derivs_xmin, derivs_xmax )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
  !  type (sll_t_spline_1d)             , intent(inout) :: spline
    real(wp)                           , intent(  out) :: &  ! not using spline to separate data and basis
    bcoef(1-MERGE(1+self%bspl%degree/2,0,self%bspl%periodic):self%bspl%nbasis+MERGE(1+self%bspl%degree/2,0,self%bspl%periodic))
    real(wp)                           , intent(in   ) :: gtau(:)
    real(wp),                  optional, intent(in   ) :: derivs_xmin(:)
    real(wp),                  optional, intent(in   ) :: derivs_xmax(:)

    character(len=*), parameter :: &
      this_sub_name = "sll_t_spline_interpolator_1d % compute_interpolant"

    integer :: i

    SLL_ASSERT( size(gtau) == self%bspl%nbasis - self%nbc_xmin - self%nbc_xmax )
!    SLL_ASSERT( spline % belongs_to_space( self % bspl ) )

    associate ( nbasis   =>   self % bspl % nbasis, &
                degree   =>   self % bspl % degree, &
                nbc_xmin =>   self % nbc_xmin     , &
                nbc_xmax =>   self % nbc_xmax     )!, &
!                bcoef    => spline % bcoef )

      ! Special case: linear spline
      if (degree == 1) then
        bcoef(1:nbasis) = gtau(1:nbasis)
        ! Periodic only: "wrap around" coefficients onto extended array
        if (self%bspl%periodic) then
          bcoef(0)        = bcoef(nbasis)
          bcoef(nbasis+1) = bcoef(1)
        end if
        return
      end if

      ! Hermite boundary conditions at xmin, if any
      ! NOTE: For consistency with the linear system, the i-th derivative
      !       provided by the user must be multiplied by dx^i
      if ( self%bc_xmin == sll_p_hermite ) then
        if ( present( derivs_xmin ) ) then
          bcoef(1:nbc_xmin) = &
                     [(derivs_xmin(i)*self%dx**(i+self%odd-1), i=nbc_xmin,1,-1)]
        else
          SLL_ERROR(this_sub_name,"Hermite BC at xmin requires derivatives")
        end if
      end if

      ! Interpolation points
      bcoef(nbc_xmin+1:nbasis-nbc_xmax) = gtau(:)

      ! Hermite boundary conditions at xmax, if any
      ! NOTE: For consistency with the linear system, the i-th derivative
      !       provided by the user must be multiplied by dx^i
      if ( self%bc_xmax == sll_p_hermite ) then
        if ( present( derivs_xmax ) ) then
          bcoef(nbasis-nbc_xmax+1:nbasis) = &
                        [(derivs_xmax(i)*self%dx**(i+self%odd-1), i=1,nbc_xmax)]
        else
          SLL_ERROR(this_sub_name,"Hermite BC at xmax requires derivatives")
        end if
      end if

      ! Solve linear system and compute coefficients
      call self % matrix % solve_inplace(1, bcoef(1:nbasis) )

      ! Periodic only: "wrap around" coefficients onto extended array
      if (self%bc_xmin == sll_p_periodic) then
        associate( g => 1+degree/2 )
          bcoef(1-g:0)             = bcoef(nbasis-g+1:nbasis)
          bcoef(nbasis+1:nbasis+g) = bcoef(1:g)
        end associate
      end if

    end associate ! nbasis, degree, bcoef

  end subroutine s_spline_interpolator_1d__compute_interpolant

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!                            PRIVATE SUBROUTINES
  !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! TODO: move subroutine to init()
  !-----------------------------------------------------------------------------
  !> @brief        Private subroutine for assembling and factorizing linear
  !>               system for any kind of boundary conditions at xmin and xmax
  !> @param[in]    self   1D spline interpolator object
  !> @param[inout] matrix generic 'spline' matrix (dense/banded/periodic-banded)
  !-----------------------------------------------------------------------------
  subroutine s_build_system( self, matrix )

    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    class(sll_c_spline_matrix)         , intent(inout) :: matrix

    integer  :: i,j,d,s
    integer  :: j0,d0
    integer  :: order
    real(wp) :: x
    real(wp) :: values(self%bspl%degree+1)
    real(wp), allocatable :: derivs(:,:)

    ! NEW
    integer :: jmin

    associate( nbasis   => self % bspl % nbasis, &
               degree   => self % bspl % degree, &
               nbc_xmin => self % nbc_xmin     , &
               nbc_xmax => self % nbc_xmax )

      if ( any( [self%bc_xmin,self%bc_xmax] == sll_p_hermite ) ) then
        allocate ( derivs (0:degree/2, 1:degree+1) )
      end if

      ! Hermite boundary conditions at xmin, if any
      if ( self%bc_xmin == sll_p_hermite ) then
        x = self%bspl%xmin
        call self % bspl % eval_basis_and_n_derivs( x, nbc_xmin, derivs, jmin )

        ! In order to improve the condition number of the matrix, we normalize all
        ! derivatives by multiplying the i-th derivative by dx^i
        associate( h => [(self%dx**i, i=1, ubound(derivs,1))] )
          do j = lbound(derivs,2), ubound(derivs,2)
            derivs(1:,j) = derivs(1:,j) * h(1:)
          end do
        end associate

        do i = 1, nbc_xmin
          ! iterate only to deg as last bspline is 0
          order = nbc_xmin-i+self%odd
          do j = 1, degree
            call matrix % set_element( i, j, derivs(order,j) )
          end do
        end do

      end if

      ! Interpolation points
      do i = nbc_xmin+1, nbasis-nbc_xmax
        x = self%tau(i-nbc_xmin)
        call self % bspl % eval_basis( x, values, jmin )
        do s = 1, degree+1
          j = modulo( jmin+s-2, nbasis ) + 1
          call matrix % set_element( i, j, values(s) )
        end do
      end do

      ! Hermite boundary conditions at xmax, if any
      if ( self%bc_xmax == sll_p_hermite ) then
        x = self%bspl%xmax
        call self % bspl % eval_basis_and_n_derivs( x, nbc_xmax, derivs, jmin )

        ! In order to improve the condition number of the matrix, we normalize all
        ! derivatives by multiplying the i-th derivative by dx^i
        associate( h => [(self%dx**i, i=1, ubound(derivs,1))] )
          do j = lbound(derivs,2), ubound(derivs,2)
            derivs(1:,j) = derivs(1:,j) * h(1:)
          end do
        end associate

        do i = nbasis-nbc_xmax+1, nbasis
          order = i-(nbasis-nbc_xmax+1)+self%odd
          j0 = nbasis-degree
          d0 = 1
          do s = 1, degree
            j = j0 + s
            d = d0 + s
            call matrix % set_element( i, j, derivs(order,d) )
          end do
        end do

      end if

      if ( allocated( derivs ) ) deallocate( derivs )

    end associate ! nbasis, degree

  end subroutine s_build_system

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!                       PRIVATE SUBROUTINES, UNIFORM
  !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !-----------------------------------------------------------------------------
  subroutine s_compute_interpolation_points_uniform( self, tau )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    real(wp),               allocatable, intent(  out) :: tau(:)

    integer :: i, ntau, isum
    integer, allocatable :: iknots(:)

    associate( nbasis   => self % bspl % nbasis, &
               ncells   => self % bspl % ncells, &
               degree   => self % bspl % degree, &
               xmin     => self % bspl % xmin  , &
               xmax     => self % bspl % xmax  , &
               dx       => self % dx           , &
               bc_xmin  => self %  bc_xmin     , &
               bc_xmax  => self %  bc_xmax     , &
               nbc_xmin => self % nbc_xmin     , &
               nbc_xmax => self % nbc_xmax )

      ! Determine size of tau and allocate tau
      ntau = nbasis - nbc_xmin - nbc_xmax
      allocate( tau(1:ntau) )

      if ( bc_xmin == sll_p_periodic ) then

      ! Periodic case:
      !   . for odd  degree, interpolation points are breakpoints (last excluded)
      !   . for even degree, interpolation points are cell centers
      associate( shift => merge( 0.5_wp, 0.0_wp, self%odd == 0 ) )
        tau = [(xmin + (real(i-1,wp)+shift)*dx, i=1,ntau)]
      end associate

      else

      ! Non-periodic case: create array of temporary knots (integer shifts only)
      ! in order to compute interpolation points using Greville-style averaging:
      ! tau(i) = xmin + average(knots(i+1-degree:i)) * dx
      allocate( iknots (2-degree:ntau) )

      ! Additional knots near x=xmin
      associate( r => 2-degree, s => -nbc_xmin )
        select case (bc_xmin)
          case (sll_p_greville); iknots(r:s) = 0
          case (sll_p_hermite ); iknots(r:s) = [(i,i=r-s-1,-1)]
        end select
      end associate

      ! Knots inside the domain
      associate( r => -nbc_xmin+1, s => -nbc_xmin+1+ncells )
        iknots(r:s) = [(i,i=0,ncells)]
      end associate

      ! Additional knots near x=xmax
      associate( r => -nbc_xmin+1+ncells+1, s => ntau )
        select case (bc_xmax)
          case (sll_p_greville); iknots(r:s) = ncells
          case (sll_p_hermite ); iknots(r:s) = [(i,i=ncells+1,ncells+1+s-r)]
        end select
      end associate

      ! Compute interpolation points using Greville-style averaging
      associate( inv_deg => 1.0_wp / real( degree, wp ) )
        do i = 1, ntau
          isum = sum( iknots(i+1-degree:i) )
          if (modulo( isum, degree ) == 0) then
            tau(i) = xmin + real(isum/degree,wp) * dx
          else
            tau(i) = xmin + real(isum,wp) * inv_deg * dx
          end if
        end do
      end associate

      ! Non-periodic case, odd degree: fix round-off issues
      if ( self%odd == 1 ) then
        tau(1)    = xmin
        tau(ntau) = xmax
      end if

      deallocate( iknots )

      end if  ! (bc_xmin == sll_p_periodic)

    end associate

  end subroutine s_compute_interpolation_points_uniform

  !-----------------------------------------------------------------------------
  subroutine s_compute_num_diags_uniform( self, kl, ku )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    integer                            , intent(  out) :: kl
    integer                            , intent(  out) :: ku

    associate( degree => self % bspl % degree )

      ! FIXME: In Hermite case ku and kl computed in general case when derivatives
      !        of B-splines do not vanish at boundary
      ! TODO: reduce kl and ku to take advantage of uniform grid
      select case( self % bc_xmin )
        case ( sll_p_periodic ); ku = (degree+1)/2
        case ( sll_p_hermite  ); ku = max( (degree+1)/2, degree-1 )
        case ( sll_p_greville ); ku = degree
      end select

      select case( self % bc_xmax )
        case ( sll_p_periodic ); kl = (degree+1)/2
        case ( sll_p_hermite  ); kl = max( (degree+1)/2, degree-1 )
        case ( sll_p_greville ); kl = degree
      end select

    end associate

  end subroutine s_compute_num_diags_uniform

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!                       PRIVATE SUBROUTINES, NON-UNIFORM
  !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !-----------------------------------------------------------------------------
  subroutine s_compute_interpolation_points_non_uniform( self, tau )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    real(wp),               allocatable, intent(  out) :: tau(:)

   integer :: i, ntau
   real(wp), allocatable :: temp_knots(:)

   associate( nbasis   => self % bspl % nbasis, &
              ncells   => self % bspl % ncells, &
              degree   => self % bspl % degree, &
              xmin     => self % bspl % xmin  , &
              xmax     => self % bspl % xmax  , &
              dx       => self % dx           , &
              bc_xmin  => self %  bc_xmin     , &
              bc_xmax  => self %  bc_xmax     , &
              nbc_xmin => self % nbc_xmin     , &
              nbc_xmax => self % nbc_xmax )

   associate( breaks   => self % bspl % knots(1:ncells+1) )

      ! Determine size of tau and allocate tau
      ntau = nbasis - nbc_xmin - nbc_xmax
      allocate( tau(1:ntau) )

      ! Array of temporary knots needed to compute interpolation points
      ! using Greville-style averaging: tau(i) = average(temp_knots(i+1-degree:i))
      allocate( temp_knots( 2-degree : ntau ) )

      if ( bc_xmin == sll_p_periodic ) then

        associate( k => degree/2 )
          temp_knots(:) = self%bspl%knots(2-degree+k:ntau+k)
        end associate

      else

        associate( r => 2-degree, s => -nbc_xmin )
          select case (bc_xmin)
          case (sll_p_greville); temp_knots(r:s) = breaks(1)
          case (sll_p_hermite ); temp_knots(r:s) = 2.0_wp*breaks(1) - breaks(2+s-r:2:-1)
          end select
        end associate

        associate( r => -nbc_xmin+1, s => -nbc_xmin+1+ncells )
          temp_knots(r:s) = breaks(:)
        end associate

        associate( r => -nbc_xmin+1+ncells+1, s => ntau )
          select case (bc_xmax)
          case (sll_p_greville)
            temp_knots(r:s) = breaks(ncells+1)
          case (sll_p_hermite )
            temp_knots(r:s) = 2.0_wp*breaks(ncells+1) - breaks(ncells:ncells+r-s:-1)
          end select
        end associate

      end if

      ! Compute interpolation points using Greville-style averaging
      associate( inv_deg => 1.0_wp / real( degree, wp ) )
        do i = 1, ntau
          tau(i) = sum( temp_knots(i+1-degree:i) ) * inv_deg
        end do
      end associate

      ! Periodic case: apply periodic BCs to interpolation points
      if ( bc_xmin == sll_p_periodic ) then
        tau(:) = modulo( tau(:)-xmin, xmax-xmin ) + xmin

      ! Non-periodic case, odd degree: fix round-off issues
      else if ( self%odd == 1 ) then
        tau(1)    = xmin
        tau(ntau) = xmax
      end if

    end associate ! breaks
    end associate ! all other variables

  end subroutine s_compute_interpolation_points_non_uniform

  !-----------------------------------------------------------------------------
  subroutine s_compute_num_diags_non_uniform( self, kl, ku )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    integer                            , intent(  out) :: kl
    integer                            , intent(  out) :: ku

    integer  :: i,j,s,d,icell
    real(wp) :: x

    ku = 0
    kl = 0

    associate( nbasis   => self % bspl % nbasis, &
               degree   => self % bspl % degree, &
               nbc_xmin => self % nbc_xmin     , &
               nbc_xmax => self % nbc_xmax )

      if (self%bc_xmin == sll_p_periodic) then

        do i = 1, nbasis
          x = self % tau(i)
          icell = self % bspl % find_cell( x )
          do s = 1, degree+1
            j = modulo(icell-self%bspl%offset-2+s,nbasis)+1
            d = j-i
            if (d >  nbasis/2) then; d = d-nbasis; else &
            if (d < -nbasis/2) then; d = d+nbasis; end if
            ku = max( ku,  d )
            kl = max( kl, -d )
          end do
        end do

      else

        do i = nbc_xmin+1, nbasis-nbc_xmax
          x = self % tau(i-nbc_xmin)
          icell = self % bspl % find_cell( x )
          do s = 1, degree+1
            j = icell-1+s
            d = j-i
            ku = max( ku,  d )
            kl = max( kl, -d )
          end do
        end do
        ! FIXME: In Hermite case ku and kl computed in general case where
        !        derivatives of B-splines do not vanish at boundary
        ku = ku + nbc_xmin
        kl = kl + nbc_xmax

      end if

    end associate

  end subroutine s_compute_num_diags_non_uniform

end module sll_m_spline_interpolator_1d