dbvalu Subroutine

private pure subroutine dbvalu(t, a, n, k, ideriv, x, inbv, work, iflag, val, extrap)

Evaluates the b-representation (t,a,n,k) of a b-spline at x for the function value on ideriv=0 or any of its derivatives on ideriv=1,2,...,k-1. right limiting values (right derivatives) are returned except at the right end point x=t(n+1) where left limiting values are computed. the spline is defined on t(k) x t(n+1). dbvalu returns a fatal error message when x is outside of this interval.

To compute left derivatives or left limiting values at a knot t(i), replace n by i-1 and set x=t(i), i=k+1,n+1.

Error Conditions

  • improper input

History

  • bvalue written by carl de boor [5]
  • dbvalu author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • 000330 modified array declarations. (jec)
  • Jacob Williams, 2/24/2015 : extensive refactoring of CMLIB routine.

Arguments

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

knot vector of length n+k
real(kind=wp), intent(in), dimension(n) :: a

b-spline coefficient vector of length n
integer(kind=ip), intent(in) :: n

number of b-spline coefficients. (sum of knot multiplicities-k)
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. ideriv = 0 returns the b-spline value
real(kind=wp), intent(in) :: x

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

an initialization parameter which must be set to 1 the first time dbvalu is called. inbv contains information for efficient processing after the initial call and inbv must not be changed by the user. distinct splines require distinct inbv parameters.
real(kind=wp), intent(inout), dimension(:) :: work

work vector of length at least 3*k
integer(kind=ip), intent(out) :: iflag

status flag:

  • 0: no errors
  • 401: k does not satisfy k 1
  • 402: n does not satisfy n k
  • 403: ideriv does not satisfy 0 ideriv k
  • 404: x is not greater than or equal to t(k)
  • 405: x is not less than or equal to t(n+1)
  • 406: a left limiting value cannot be obtained at t(k)
real(kind=wp), intent(out) :: val

the interpolated value
logical, intent(in), optional :: extrap

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

Calls

proc~~dbvalu~~CallsGraph proc~dbvalu bspline_sub_module::dbvalu proc~dintrv bspline_sub_module::dintrv proc~dbvalu->proc~dintrv proc~get_temp_x_for_extrap bspline_sub_module::get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap

Called by

proc~~dbvalu~~CalledByGraph proc~dbvalu bspline_sub_module::dbvalu proc~db1val_alt bspline_sub_module::db1val_alt proc~db1val_alt->proc~dbvalu proc~db1val_default bspline_sub_module::db1val_default proc~db1val_default->proc~dbvalu proc~db2val bspline_sub_module::db2val proc~db2val->proc~dbvalu proc~db3val bspline_sub_module::db3val proc~db3val->proc~dbvalu proc~db4val bspline_sub_module::db4val proc~db4val->proc~dbvalu proc~db5val bspline_sub_module::db5val proc~db5val->proc~dbvalu proc~db6val bspline_sub_module::db6val proc~db6val->proc~dbvalu proc~dbsgq8 bspline_sub_module::dbsgq8 proc~dbsgq8->proc~dbvalu proc~dbsqad bspline_sub_module::dbsqad proc~dbsqad->proc~dbvalu interface~db1val bspline_sub_module::db1val interface~db1val->proc~db1val_alt interface~db1val->proc~db1val_default proc~db1sqad bspline_sub_module::db1sqad proc~db1sqad->proc~dbsqad proc~dbfqad bspline_sub_module::dbfqad proc~dbfqad->proc~dbsgq8 proc~evaluate_2d bspline_oo_module::bspline_2d%evaluate_2d proc~evaluate_2d->proc~db2val proc~evaluate_3d bspline_oo_module::bspline_3d%evaluate_3d proc~evaluate_3d->proc~db3val proc~evaluate_4d bspline_oo_module::bspline_4d%evaluate_4d proc~evaluate_4d->proc~db4val proc~evaluate_5d bspline_oo_module::bspline_5d%evaluate_5d proc~evaluate_5d->proc~db5val proc~evaluate_6d bspline_oo_module::bspline_6d%evaluate_6d proc~evaluate_6d->proc~db6val proc~db1fqad bspline_sub_module::db1fqad proc~db1fqad->proc~dbfqad proc~evaluate_1d bspline_oo_module::bspline_1d%evaluate_1d proc~evaluate_1d->interface~db1val proc~integral_1d bspline_oo_module::bspline_1d%integral_1d proc~integral_1d->proc~db1sqad proc~fintegral_1d bspline_oo_module::bspline_1d%fintegral_1d proc~fintegral_1d->proc~db1fqad

Source Code

    pure subroutine dbvalu(t,a,n,k,ideriv,x,inbv,work,iflag,val,extrap)

    implicit none

    real(wp),intent(out)             :: val     !! the interpolated value
    integer(ip),intent(in)           :: n       !! number of b-spline coefficients.
                                                !! (sum of knot multiplicities-`k`)
    real(wp),dimension(:),intent(in) :: t       !! knot vector of length `n+k`
    real(wp),dimension(n),intent(in) :: a       !! b-spline coefficient vector of length `n`
    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`.
                                                !! `ideriv = 0` returns the b-spline value
    real(wp),intent(in)              :: x       !! argument, `t(k) <= x <= t(n+1)`
    integer(ip),intent(inout)        :: inbv    !! an initialization parameter which must be set
                                                !! to 1 the first time [[dbvalu]] is called.
                                                !! `inbv` contains information for efficient processing
                                                !! after the initial call and `inbv` must not
                                                !! be changed by the user.  distinct splines require
                                                !! distinct `inbv` parameters.
    real(wp),dimension(:),intent(inout) :: work !! work vector of length at least `3*k`
    integer(ip),intent(out)          :: iflag   !! status flag:
                                                !!
                                                !! * 0: no errors
                                                !! * 401: `k` does not satisfy `k` \( \ge \) 1
                                                !! * 402: `n` does not satisfy `n` \( \ge \) `k`
                                                !! * 403: `ideriv` does not satisfy 0 \( \le \) `ideriv` \(<\) `k`
                                                !! * 404: `x` is not greater than or equal to `t(k)`
                                                !! * 405: `x` is not less than or equal to `t(n+1)`
                                                !! * 406: a left limiting value cannot be obtained at `t(k)`
    logical,intent(in),optional :: extrap   !! if extrapolation is allowed
                                            !! (if not present, default is False)

    integer(ip) :: i,iderp1,ihi,ihmkmj,ilo,imk,imkpj,ipj,&
               ip1,ip1mj,j,jj,j1,j2,kmider,kmj,km1,kpk,mflag
    real(wp) :: fkmj
    real(wp) :: xt
    logical :: extrapolation_allowed  !! if extrapolation is allowed

    val = 0.0_wp

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

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

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

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

    ! make a temp copy of x (for computing the
    ! interval) in case extrapolation is allowed
    if (extrapolation_allowed) then
        if (x<t(k)) then
            xt = t(k)
        else if (x>t(n+1_ip)) then
            xt = t(n+1_ip)
        else
            xt = x
        end if
    else
        xt = x
    end if

    kmider = k - ideriv

    ! 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)).

    km1 = k - 1_ip
    call dintrv(t, n+1, xt, inbv, i, mflag)
    if (xt<t(k)) then
        iflag = 404_ip  ! dbvalu - 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  ! dbvalu - x is not less than or equal to t(n+1)
            return
        end if

        do
            if (i==k) then
                iflag = 406_ip  ! dbvalu - 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

    ! difference the coefficients *ideriv* times
    ! work(i) = aj(i), work(k+i) = dp(i), work(k+k+i) = dm(i), i=1.k

    imk = i - k
    do j=1_ip,k
        imkpj = imk + j
        work(j) = a(imkpj)
    end do

    if (ideriv/=0_ip) then
        do j=1_ip,ideriv
            kmj = k - j
            fkmj = real(kmj,wp)
            do jj=1_ip,kmj
                ihi = i + jj
                ihmkmj = ihi - kmj
                work(jj) = (work(jj+1_ip)-work(jj))/(t(ihi)-t(ihmkmj))*fkmj
            end do
        end do
    end if

    ! compute value at *x* in (t(i),(t(i+1)) of ideriv-th derivative,
    ! given its relevant b-spline coeff. in aj(1),...,aj(k-ideriv).

    if (ideriv/=km1) then
        ip1 = i + 1_ip
        kpk = k + k
        j1 = k + 1_ip
        j2 = kpk + 1_ip
        do j=1_ip,kmider
            ipj = i + j
            work(j1) = t(ipj) - x
            ip1mj = ip1 - j
            work(j2) = x - t(ip1mj)
            j1 = j1 + 1_ip
            j2 = j2 + 1_ip
        end do
        iderp1 = ideriv + 1_ip
        do j=iderp1,km1
            kmj = k - j
            ilo = kmj
            do jj=1_ip,kmj
                work(jj) = (work(jj+1_ip)*work(kpk+ilo)+work(jj)*&
                            work(k+jj))/(work(kpk+ilo)+work(k+jj))
                ilo = ilo - 1
            end do
        end do
    end if

    iflag = 0_ip
    val = work(1_ip)

    end subroutine dbvalu