spline_derivative Subroutine

private pure subroutine spline_derivative(t, c, n, k, td, d)

The derivative of a b-spline, as a b-spline of order k-1 with n-1 coefficients and knots t(2:n+k-1), times a positive constant (so it has the same roots).

The constant is chosen (as a power of 2, so the scaling is exact) to avoid overflow for large coefficients or closely spaced knots. It is only used for the roots of .

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(n+k) :: t

knot array
real(kind=wp), intent(in), dimension(n) :: c

b-spline coefficient array
integer(kind=ip), intent(in) :: n

number of b-spline coefficients
integer(kind=ip), intent(in) :: k

order of the b-spline, k >= 2
real(kind=wp), intent(out), dimension(:), allocatable :: td

knots of the derivative (size n+k-2)
real(kind=wp), intent(out), dimension(:), allocatable :: d

coefficients of the derivative (size n-1)

Called by

proc~~spline_derivative~~CalledByGraph proc~spline_derivative spline_derivative proc~db1extremum db1extremum proc~db1extremum->proc~spline_derivative proc~db1root db1root proc~db1root->proc~spline_derivative proc~db1max db1max proc~db1max->proc~db1extremum proc~db1min db1min proc~db1min->proc~db1extremum proc~roots_1d bspline_1d%roots_1d proc~roots_1d->proc~db1root proc~inverse_1d bspline_1d%inverse_1d proc~inverse_1d->proc~roots_1d proc~maximum_1d bspline_1d%maximum_1d proc~maximum_1d->proc~db1max proc~minimum_1d bspline_1d%minimum_1d proc~minimum_1d->proc~db1min

Source Code

    pure subroutine spline_derivative(t,c,n,k,td,d)

    implicit none

    integer(ip),intent(in)                        :: n   !! number of b-spline coefficients
    integer(ip),intent(in)                        :: k   !! order of the b-spline, `k >= 2`
    real(wp),dimension(n+k),intent(in)            :: t   !! knot array
    real(wp),dimension(n),intent(in)              :: c   !! b-spline coefficient array
    real(wp),dimension(:),allocatable,intent(out) :: td  !! knots of the derivative (size `n+k-2`)
    real(wp),dimension(:),allocatable,intent(out) :: d   !! coefficients of the derivative (size `n-1`)

    integer(ip) :: i   !! counter
    real(wp) :: dt     !! support width of a b-spline basis function: `t(i+k-1) - t(i)` (scaled)
    real(wp) :: cs     !! scale factor for the coefficients (a power of 2)
    real(wp) :: ts     !! scale factor for the knots (a power of 2)

    cs = pow2(maxval(abs(c)))
    ts = pow2(max(abs(t(1)),abs(t(n+k))))
    td = t(2_ip:n+k-1_ip)
    allocate(d(n-1_ip))
    do i = 2_ip, n
        dt = t(i+k-1_ip)/ts - t(i)/ts
        if (dt>0.0_wp) then
            d(i-1_ip) = real(k-1_ip,wp) * (c(i)/cs - c(i-1_ip)/cs) / dt
        else
            d(i-1_ip) = 0.0_wp  ! the basis function is identically zero
        end if
    end do

    contains

        pure function pow2(x) result(p)
        !! The largest power of 2 not greater than `x` (1 if `x` is 0).
        !! (so it is representable for any finite `x`)
        real(wp),intent(in) :: x  !! the value (`x >= 0`)
        real(wp) :: p             !! the power of 2
        if (x>0.0_wp) then
            p = scale(1.0_wp, exponent(x)-1)
        else
            p = 1.0_wp
        end if
        end function pow2

    end subroutine spline_derivative