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 nderiv1
right derivatives at x
, i=1,k, j=1,nderiv
.
These basis functions have indices ileftk+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.
Note
DBSPVD
is the BSPLVD
routine of the reference.
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(in),  dimension(:)  ::  t 
knot vector of length 

integer(kind=ip),  intent(in)  ::  k 
order of the bspline, 

integer(kind=ip),  intent(in)  ::  nderiv 
number of derivatives = 

real(kind=wp),  intent(in)  ::  x 
argument of basis functions,


integer(kind=ip),  intent(in)  ::  ileft 
largest integer such that


integer(kind=ip),  intent(in)  ::  ldvnik 
leading dimension of matrix 

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

real(kind=wp),  intent(out),  dimension(*)  ::  work 
a work vector of length 

integer(kind=ip),  intent(out)  ::  iflag 
status flag:

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 bspline basis functions !! `n` = sum of knot multiplicitiesk integer(ip),intent(in) :: k !! order of the bspline, `k >= 1` integer(ip),intent(in) :: nderiv !! number of derivatives = `nderiv1`, !! `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,k1 ! a(i,k) = work(i+k*(k1)/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