dbspvd Subroutine

private pure subroutine dbspvd(t, k, nderiv, x, ileft, ldvnik, vnikx, work, iflag)

DBSPVD calculates the value and all derivatives of order less than nderiv of all basis functions which do not (possibly) vanish at x. ileft is input such that t(ileft) <= x < t(ileft+1). A call to dintrv(t,n+1,x, ilo,ileft,mflag) will produce the proper ileft. The output of dbspvd is a matrix vnikx(i,j) of dimension at least (k,nderiv) whose columns contain the k nonzero basis functions and their nderiv-1 right derivatives at x, i=1,k, j=1,nderiv. These basis functions have indices ileft-k+i, i=1,k, k <= ileft <= n. The nonzero part of the i-th basis function lies in (t(i),t(i+k)), i=1,n).

If x=t(ileft+1) then vnikx contains left limiting values (left derivatives) at t(ileft+1). In particular, ileft = n produces left limiting values at the right end point x=t(n+1). To obtain left limiting values at t(i), i=k+1,n+1, set x = next lower distinct knot, call dintrv to get ileft, set x=t(i), and then call dbspvd.

History

  • Written by Carl de Boor and modified by D. E. Amos
  • date written 800901
  • revision date 820801
  • 000330 Modified array declarations. (JEC)
  • Jacob Williams, 8/30/2018 : refactored to modern Fortran.

Note

DBSPVD is the BSPLVD routine of the reference.

Arguments

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

knot vector of length n+k, where n = number of b-spline basis functions n = sum of knot multiplicities-k
integer(kind=ip), intent(in) :: k

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

number of derivatives = nderiv-1, 1 <= nderiv <= k
real(kind=wp), intent(in) :: x

argument of basis functions, t(k) <= x <= t(n+1)
integer(kind=ip), intent(in) :: ileft

largest integer such that t(ileft) <= x < t(ileft+1)
integer(kind=ip), intent(in) :: ldvnik

leading dimension of matrix vnikx
real(kind=wp), intent(out), dimension(ldvnik,nderiv) :: vnikx

matrix of dimension at least (k,nderiv) containing the nonzero basis functions at x and their derivatives columnwise.
real(kind=wp), intent(out), dimension(*) :: work

a work vector of length (k+1)*(k+2)/2
integer(kind=ip), intent(out) :: iflag

status flag:

  • 0: no errors
  • 3001: k does not satisfy k>=1
  • 3002: nderiv does not satisfy 1<=nderiv<=k
  • 3003: ldvnik does not satisfy ldvnik>=k

Calls

proc~~dbspvd~~CallsGraph proc~dbspvd bspline_sub_module::dbspvd proc~dbspvn bspline_sub_module::dbspvn proc~dbspvd->proc~dbspvn

Called by

proc~~dbspvd~~CalledByGraph proc~dbspvd bspline_sub_module::dbspvd proc~dbint4 bspline_sub_module::dbint4 proc~dbint4->proc~dbspvd proc~db1ink_alt bspline_sub_module::db1ink_alt proc~db1ink_alt->proc~dbint4 proc~db1ink_alt_2 bspline_sub_module::db1ink_alt_2 proc~db1ink_alt_2->proc~dbint4 interface~db1ink bspline_sub_module::db1ink interface~db1ink->proc~db1ink_alt interface~db1ink->proc~db1ink_alt_2 proc~initialize_1d_auto_knots bspline_oo_module::bspline_1d%initialize_1d_auto_knots proc~initialize_1d_auto_knots->interface~db1ink proc~initialize_1d_specify_knots bspline_oo_module::bspline_1d%initialize_1d_specify_knots proc~initialize_1d_specify_knots->interface~db1ink proc~bspline_1d_constructor_auto_knots bspline_oo_module::bspline_1d_constructor_auto_knots proc~bspline_1d_constructor_auto_knots->proc~initialize_1d_auto_knots proc~bspline_1d_constructor_specify_knots bspline_oo_module::bspline_1d_constructor_specify_knots proc~bspline_1d_constructor_specify_knots->proc~initialize_1d_specify_knots interface~bspline_1d bspline_oo_module::bspline_1d interface~bspline_1d->proc~bspline_1d_constructor_auto_knots interface~bspline_1d->proc~bspline_1d_constructor_specify_knots

Source Code

    pure subroutine dbspvd(t,k,nderiv,x,ileft,ldvnik,vnikx,work,iflag)

    implicit none

    real(wp),dimension(:),intent(in)              :: t       !! knot vector of length `n+k`, where
                                                             !! `n` = number of b-spline basis functions
                                                             !! `n` = sum of knot multiplicities-k
    integer(ip),intent(in)                        :: k       !! order of the b-spline, `k >= 1`
    integer(ip),intent(in)                        :: nderiv  !! number of derivatives = `nderiv-1`,
                                                             !! `1 <= nderiv <= k`
    real(wp),intent(in)                           :: x       !! argument of basis functions,
                                                             !! `t(k) <= x <= t(n+1)`
    integer(ip),intent(in)                        :: ileft   !! largest integer such that
                                                             !! `t(ileft) <= x < t(ileft+1)`
    integer(ip),intent(in)                        :: ldvnik  !! leading dimension of matrix `vnikx`
    real(wp),dimension(ldvnik,nderiv),intent(out) :: vnikx   !! matrix of dimension at least `(k,nderiv)`
                                                             !! containing the nonzero basis functions
                                                             !! at `x` and their derivatives columnwise.
    real(wp),dimension(*),intent(out)             :: work    !! a work vector of length `(k+1)*(k+2)/2`
    integer(ip),intent(out)                       :: iflag   !! status flag:
                                                             !!
                                                             !! * 0: no errors
                                                             !! * 3001: `k` does not satisfy `k>=1`
                                                             !! * 3002: `nderiv` does not satisfy `1<=nderiv<=k`
                                                             !! * 3003: `ldvnik` does not satisfy `ldvnik>=k`

    integer(ip) :: i,ideriv,ipkmd,j,jj,jlow,jm,jp1mid,kmd,kp1,l,ldummy,m,mhigh,iwork
    real(wp) :: factor, fkmd, v

    ! dimension t(ileft+k), work((k+1)*(k+2)/2)
    ! a(i,j) = work(i+j*(j+1)/2),  i=1,j+1  j=1,k-1
    ! a(i,k) = work(i+k*(k-1)/2)  i=1.k
    ! work(1) and work((k+1)*(k+2)/2) are not used.

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

    if (nderiv<1 .or. nderiv>k) then
        iflag = 3002_ip ! nderiv does not satisfy 1<=nderiv<=k
        return
    end if

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

    iflag = 0_ip

    ideriv = nderiv
    kp1 = k + 1
    jj = kp1 - ideriv
    call dbspvn(t, jj, k, 1_ip, x, ileft, vnikx, work, iwork, iflag)
    if (iflag/=0 .or. ideriv==1) return
    mhigh = ideriv
    do m=2,mhigh
        jp1mid = 1
        do j=ideriv,k
            vnikx(j,ideriv) = vnikx(jp1mid,1)
            jp1mid = jp1mid + 1
        end do
        ideriv = ideriv - 1
        jj = kp1 - ideriv
        call dbspvn(t, jj, k, 2_ip, x, ileft, vnikx, work, iwork, iflag)
        if (iflag/=0) return
    end do

    jm = kp1*(kp1+1)/2
    do l = 1,jm
        work(l) = 0.0_wp
    end do
    ! a(i,i) = work(i*(i+3)/2) = 1.0       i = 1,k
    l = 2
    j = 0
    do i = 1,k
        j = j + l
        work(j) = 1.0_wp
        l = l + 1
    end do
    kmd = k
    do m=2,mhigh
        kmd = kmd - 1
        fkmd = real(kmd,wp)
        i = ileft
        j = k
        jj = j*(j+1)/2
        jm = jj - j
        do ldummy=1,kmd
            ipkmd = i + kmd
            factor = fkmd/(t(ipkmd)-t(i))
            do l=1,j
                work(l+jj) = (work(l+jj)-work(l+jm))*factor
            end do
            i = i - 1
            j = j - 1
            jj = jm
            jm = jm - j
        end do

        do i=1,k
            v = 0.0_wp
            jlow = max(i,m)
            jj = jlow*(jlow+1)/2
            do j=jlow,k
                v = work(i+jj)*vnikx(j,m) + v
                jj = jj + j + 1
            end do
            vnikx(i,m) = v
        end do
    end do

    end subroutine dbspvd