s_build_system Subroutine

  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