dbspvn Subroutine

private subroutine dbspvn(t, jhigh, k, index, x, ileft, vnikx, work, iwork, iflag)

Calculates the value of all (possibly) nonzero basis functions at x of order max(jhigh,(j+1)*(index-1)), where t(k) <= x <= t(n+1) and j=iwork is set inside the routine on the first call when index=1. ileft is such that t(ileft) <= x < t(ileft+1). a call to dintrv(t,n+1,x,ilo,ileft,mflag) produces the proper ileft. dbspvn calculates using the basic algorithm needed in dbspvd. if only basis functions are desired, setting jhigh=k and index=1 can be faster than calling dbspvd, but extra coding is required for derivatives (index=2) and dbspvd is set up for this purpose.

left limiting values are set up as described in dbspvd.

Error Conditions

  • improper input

History

  • bsplvn written by carl de boor [5]
  • dbspvn 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) :: t(*)

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

order of b-spline, 1 <= jhigh <= k
integer, intent(in) :: k

highest possible order
integer, intent(in) :: index

index = 1 gives basis functions of order jhigh = 2 denotes previous entry with work, iwork values saved for subsequent calls to dbspvn.
real(kind=wp), intent(in) :: x

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

largest integer such that t(ileft) <= x < t(ileft+1)
real(kind=wp), intent(out) :: vnikx(k)

vector of length k for spline values.
real(kind=wp), intent(out) :: work(*)

a work vector of length 2*k
integer, intent(out) :: iwork

a work parameter. both work and iwork contain information necessary to continue for index = 2. when index = 1 exclusively, these are scratch variables and can be used for other purposes.
integer, intent(out) :: iflag

if 0: no errors if 201: k does not satisfy k>=1 if 202: jhigh does not satisfy 1<=jhigh<=k if 203: index is not 1 or 2 if 204: x does not satisfy t(ileft)<=x<=t(ileft+1)

Called by

proc~~dbspvn~~CalledByGraph proc~dbspvn dbspvn proc~dbintk dbintk proc~dbintk->proc~dbspvn proc~dbtpcf dbtpcf proc~dbtpcf->proc~dbintk proc~db1ink db1ink proc~db1ink->proc~dbtpcf proc~db2ink db2ink proc~db2ink->proc~dbtpcf proc~interp_surf interp_surf proc~interp_surf->proc~db2ink program~test_bspline test_bspline program~test_bspline->proc~interp_surf

Source Code

    subroutine dbspvn(t,jhigh,k,index,x,ileft,vnikx,work,iwork,iflag)

    implicit none

    real(wp),intent(in)  :: t(*)     !! knot vector of length n+k, where
                                     !! n = number of b-spline basis functions
                                     !! n = sum of knot multiplicities-k
                                     !! dimension t(ileft+jhigh)
    integer,intent(in)   :: jhigh    !! order of b-spline, 1 <= jhigh <= k
    integer,intent(in)   :: k        !! highest possible order
    integer,intent(in)   :: index    !! index = 1 gives basis functions of order jhigh
                                     !!       = 2 denotes previous entry with work, iwork
                                     !!         values saved for subsequent calls to
                                     !!         dbspvn.
    real(wp),intent(in)  :: x        !! argument of basis functions, t(k) <= x <= t(n+1)
    integer,intent(in)   :: ileft    !! largest integer such that t(ileft) <= x < t(ileft+1)
    real(wp),intent(out) :: vnikx(k) !! vector of length k for spline values.
    real(wp),intent(out) :: work(*)  !! a work vector of length 2*k
    integer,intent(out)  :: iwork    !! a work parameter.  both work and iwork contain
                                     !! information necessary to continue for index = 2.
                                     !! when index = 1 exclusively, these are scratch
                                     !! variables and can be used for other purposes.
    integer,intent(out) :: iflag     !! if   0: no errors
                                     !! if 201: k does not satisfy k>=1
                                     !! if 202: jhigh does not satisfy 1<=jhigh<=k
                                     !! if 203: index is not 1 or 2
                                     !! if 204: x does not satisfy t(ileft)<=x<=t(ileft+1)

    integer :: imjp1, ipj, jp1, jp1ml, l
    real(wp) :: vm, vmprev

    ! content of j, deltam, deltap is expected unchanged between calls.
    ! work(i) = deltap(i),
    ! work(k+i) = deltam(i), i = 1,k

    if (k<1) then
        write(error_unit,'(A)') 'dbspvn - k does not satisfy k>=1'
        iflag = 201
        return
    endif
    if (jhigh>k .or. jhigh<1) then
        write(error_unit,'(A)') 'dbspvn - jhigh does not satisfy 1<=jhigh<=k'
        iflag = 202
        return
    endif
    if (index<1 .or. index>2) then
        write(error_unit,'(A)') 'dbspvn - index is not 1 or 2'
        iflag = 203
        return
    endif
    if (x<t(ileft) .or. x>t(ileft+1)) then
        write(error_unit,'(A)') 'dbspvn - x does not satisfy t(ileft)<=x<=t(ileft+1)'
        iflag = 204
        return
    endif

    iflag = 0

    if (index==1) then
        iwork = 1
        vnikx(1) = 1.0_wp
        if (iwork>=jhigh) return
    endif

    do
        ipj = ileft + iwork
        work(iwork) = t(ipj) - x
        imjp1 = ileft - iwork + 1
        work(k+iwork) = x - t(imjp1)
        vmprev = 0.0_wp
        jp1 = iwork + 1
        do l=1,iwork
            jp1ml = jp1 - l
            vm = vnikx(l)/(work(l)+work(k+jp1ml))
            vnikx(l) = vm*work(l) + vmprev
            vmprev = vm*work(k+jp1ml)
        enddo
        vnikx(jp1) = vmprev
        iwork = jp1
        if (iwork>=jhigh) exit
    enddo

    endsubroutine dbspvn