dbvalu_weights Subroutine

private pure subroutine dbvalu_weights(t, n, k, ideriv, x, inbv, left, b, iflag, extrap)

Computes the interval index left and weights b such that the ideriv-th derivative of the b-spline (t,a,n,k) at x is

  sum(b(1:k)*a(left-k+1:left))

This gives the same result as dbvalu, but since the weights only depend on the knots, they can be reused for every coefficient vector with the same knots. This is used by the multidimensional db*val routines so that the work in each direction is only done once.

The interval selection (including left limiting values at t(n+1)), extrapolation behavior, and error codes are the same as dbvalu.

Arguments

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

knot vector of length n+k
integer(kind=ip), intent(in) :: n

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

order of the b-spline, k >= 1
integer(kind=ip), intent(in) :: ideriv

order of the derivative, 0 <= ideriv <= k-1.
real(kind=wp), intent(in) :: x

argument, t(k) <= x <= t(n+1) (unless extrapolation is allowed)
integer(kind=ip), intent(inout) :: inbv

an initialization parameter which must be set to 1 the first time this routine is called. (see dbvalu)
integer(kind=ip), intent(out) :: left

the interval index, k <= left <= n
real(kind=wp), intent(out), dimension(k) :: b

weights for coefficients left-k+1:left
integer(kind=ip), intent(out) :: iflag

status flag (same codes as dbvalu)
logical, intent(in), optional :: extrap

if extrapolation is allowed (if not present, default is False)

Calls

proc~~dbvalu_weights~~CallsGraph proc~dbvalu_weights dbvalu_weights proc~dintrv dintrv proc~dbvalu_weights->proc~dintrv proc~get_temp_x_for_extrap get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap

Called by

proc~~dbvalu_weights~~CalledByGraph proc~dbvalu_weights dbvalu_weights proc~db1val_multi db1val_multi proc~db1val_multi->proc~dbvalu_weights proc~db2val_default db2val_default proc~db2val_default->proc~dbvalu_weights proc~db2val_multi db2val_multi proc~db2val_multi->proc~dbvalu_weights proc~db3val_default db3val_default proc~db3val_default->proc~dbvalu_weights proc~db3val_multi db3val_multi proc~db3val_multi->proc~dbvalu_weights proc~db4val_default db4val_default proc~db4val_default->proc~dbvalu_weights proc~db4val_multi db4val_multi proc~db4val_multi->proc~dbvalu_weights proc~db5val_default db5val_default proc~db5val_default->proc~dbvalu_weights proc~db5val_multi db5val_multi proc~db5val_multi->proc~dbvalu_weights proc~db6val_default db6val_default proc~db6val_default->proc~dbvalu_weights proc~db6val_multi db6val_multi proc~db6val_multi->proc~dbvalu_weights interface~db1val db1val interface~db1val->proc~db1val_multi interface~db2val db2val interface~db2val->proc~db2val_default interface~db2val->proc~db2val_multi interface~db3val db3val interface~db3val->proc~db3val_default interface~db3val->proc~db3val_multi interface~db4val db4val interface~db4val->proc~db4val_default interface~db4val->proc~db4val_multi interface~db5val db5val interface~db5val->proc~db5val_default interface~db5val->proc~db5val_multi interface~db6val db6val interface~db6val->proc~db6val_default interface~db6val->proc~db6val_multi proc~evaluate_1d bspline_1d%evaluate_1d proc~evaluate_1d->interface~db1val proc~evaluate_1d_multi bspline_1d%evaluate_1d_multi proc~evaluate_1d_multi->interface~db1val proc~evaluate_2d bspline_2d%evaluate_2d proc~evaluate_2d->interface~db2val proc~evaluate_2d_multi bspline_2d%evaluate_2d_multi proc~evaluate_2d_multi->interface~db2val proc~evaluate_3d bspline_3d%evaluate_3d proc~evaluate_3d->interface~db3val proc~evaluate_3d_multi bspline_3d%evaluate_3d_multi proc~evaluate_3d_multi->interface~db3val proc~evaluate_4d bspline_4d%evaluate_4d proc~evaluate_4d->interface~db4val proc~evaluate_4d_multi bspline_4d%evaluate_4d_multi proc~evaluate_4d_multi->interface~db4val proc~evaluate_5d bspline_5d%evaluate_5d proc~evaluate_5d->interface~db5val proc~evaluate_5d_multi bspline_5d%evaluate_5d_multi proc~evaluate_5d_multi->interface~db5val proc~evaluate_6d bspline_6d%evaluate_6d proc~evaluate_6d->interface~db6val proc~evaluate_6d_multi bspline_6d%evaluate_6d_multi proc~evaluate_6d_multi->interface~db6val

Source Code

    pure subroutine dbvalu_weights(t,n,k,ideriv,x,inbv,left,b,iflag,extrap)

    implicit none

    real(wp),dimension(:),intent(in)  :: t       !! knot vector of length `n+k`
    integer(ip),intent(in)            :: n       !! number of b-spline coefficients.
    integer(ip),intent(in)            :: k       !! order of the b-spline, `k >= 1`
    integer(ip),intent(in)            :: ideriv  !! order of the derivative, `0 <= ideriv <= k-1`.
    real(wp),intent(in)               :: x       !! argument, `t(k) <= x <= t(n+1)`
                                                 !! (unless extrapolation is allowed)
    integer(ip),intent(inout)         :: inbv    !! an initialization parameter which must be set
                                                 !! to 1 the first time this routine is called.
                                                 !! (see [[dbvalu]])
    integer(ip),intent(out)           :: left    !! the interval index, `k <= left <= n`
    real(wp),dimension(k),intent(out) :: b       !! weights for coefficients `left-k+1:left`
    integer(ip),intent(out)           :: iflag   !! status flag (same codes as [[dbvalu]])
    logical,intent(in),optional       :: extrap  !! if extrapolation is allowed
                                                 !! (if not present, default is False)

    integer(ip) :: i,j,jj,kmider,kmj,mflag
    real(wp) :: xt,saved,term,fac
    real(wp),dimension(k) :: dr  !! `t(i+j) - x`
    real(wp),dimension(k) :: dl  !! `x - t(i+1-j)`
    logical :: extrapolation_allowed  !! if extrapolation is allowed

    b = 0.0_wp
    left = k

    if (k<1_ip) then
        iflag = 401_ip  ! k does not satisfy k>=1
        return
    end if

    if (n<k) then
        iflag = 402_ip  ! n does not satisfy n>=k
        return
    end if

    if (ideriv<0_ip .or. ideriv>=k) then
        iflag = 403_ip  ! ideriv does not satisfy 0<=ideriv<k
        return
    end if

    if (present(extrap)) then
        extrapolation_allowed = extrap
    else
        extrapolation_allowed = .false.
    end if

    ! temp copy of x for computing the interval
    ! (clamped if extrapolation is allowed)
    xt = x
    if (extrapolation_allowed) xt = min(max(x,t(k)),t(n+1_ip))

    ! find *i* in (k,n) such that t(i) <= x < t(i+1)
    ! (or, <= t(i+1) if t(i) < t(i+1) = t(n+1)).
    call dintrv(t, n+1_ip, xt, inbv, i, mflag)
    if (xt<t(k)) then
        iflag = 404_ip  ! x is not greater than or equal to t(k)
        return
    end if
    if (mflag/=0_ip) then
        if (xt>t(i)) then
            iflag = 405_ip  ! x is not less than or equal to t(n+1)
            return
        end if
        do
            if (i==k) then
                iflag = 406_ip  ! a left limiting value cannot be obtained at t(k)
                return
            end if
            i = i - 1_ip
            if (xt/=t(i)) exit
        end do
    end if

    ! values at x of the k-ideriv nonzero b-splines of order k-ideriv
    ! (the recurrence is a polynomial in x, so this also works when extrapolating)
    kmider = k - ideriv
    b(1_ip) = 1.0_wp
    do j=1_ip,kmider-1_ip
        dr(j) = t(i+j) - x
        dl(j) = x - t(i+1_ip-j)
        saved = 0.0_wp
        do jj=1_ip,j
            term = b(jj)/(dr(jj)+dl(j+1_ip-jj))
            b(jj) = saved + dr(jj)*term
            saved = dl(j+1_ip-jj)*term
        end do
        b(j+1_ip) = saved
    end do

    ! apply the transpose of the coefficient differencing in [[dbvalu]]
    ! so the weights apply to the original coefficients
    do j=ideriv,1_ip,-1_ip
        kmj = k - j
        b(kmj+1_ip) = 0.0_wp
        do jj=kmj,1_ip,-1_ip
            fac = b(jj)*real(kmj,wp)/(t(i+jj)-t(i+jj-kmj))
            b(jj+1_ip) = b(jj+1_ip) + fac
            b(jj) = -fac
        end do
    end do

    left = i
    iflag = 0_ip

    end subroutine dbvalu_weights