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.
| Type | Intent | Optional | 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, | intent(in) | :: | n |
number of b-spline coefficients. (sum of knot multiplicities-k) |
||
| integer, | intent(in) | :: | k |
order of the b-spline, k >= 1 |
||
| integer, | 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, | intent(inout) | :: | inbv |
an initialization parameter which must be set to 1 the first time dbvalu is called. inbv contains information for efficient process- ing after the initial call and inbv must not be changed by the user. distinct splines require distinct inbv parameters. |
||
| real(kind=wp), | dimension(:) | :: | work |
work vector of length 3*k |
||
| integer, | intent(out) | :: | iflag |
if 0: no errors if 401: k does not satisfy k>=1 if 402: n does not satisfy n>=k if 403: ideriv does not satisfy 0<=ideriv<k if 404: x is not greater than or equal to t(k) if 405: x is not less than or equal to t(n+1) if 406: a left limiting value cannot be obtained at t(k) |
real(wp) function dbvalu(t,a,n,k,ideriv,x,inbv,work,iflag) implicit none integer,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,intent(in) :: k !! order of the b-spline, k >= 1 integer,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,intent(inout) :: inbv !! an initialization parameter which must be set !! to 1 the first time dbvalu is called. !! inbv contains information for efficient process- !! ing after the initial call and inbv must not !! be changed by the user. distinct splines require !! distinct inbv parameters. real(wp),dimension(:) :: work !! work vector of length 3*k integer,intent(out) :: iflag !! if 0: no errors !! if 401: k does not satisfy k>=1 !! if 402: n does not satisfy n>=k !! if 403: ideriv does not satisfy 0<=ideriv<k !! if 404: x is not greater than or equal to t(k) !! if 405: x is not less than or equal to t(n+1) !! if 406: a left limiting value cannot be obtained at t(k) integer :: i,iderp1,ihi,ihmkmj,ilo,imk,imkpj,ipj,& ip1,ip1mj,j,jj,j1,j2,kmider,kmj,km1,kpk,mflag real(wp) :: fkmj dbvalu = 0.0_wp if (k<1) then write(error_unit,'(A)') 'dbvalu - k does not satisfy k>=1' iflag = 401 return endif if (n<k) then write(error_unit,'(A)') 'dbvalu - n does not satisfy n>=k' iflag = 402 return endif if (ideriv<0 .or. ideriv>=k) then write(error_unit,'(A)') 'dbvalu - ideriv does not satisfy 0<=ideriv<k' iflag = 403 return endif 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 call dintrv(t, n+1, x, inbv, i, mflag) if (x<t(k)) then write(error_unit,'(A)') 'dbvalu - x is not greater than or equal to t(k)' iflag = 404 return endif if (mflag/=0) then if (x>t(i)) then write(error_unit,'(A)') 'dbvalu - x is not less than or equal to t(n+1)' iflag = 405 return endif do if (i==k) then write(error_unit,'(A)') 'dbvalu - a left limiting value cannot be obtained at t(k)' iflag = 406 return endif i = i - 1 if (x/=t(i)) exit enddo endif ! 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,k imkpj = imk + j work(j) = a(imkpj) enddo if (ideriv/=0) then do j=1,ideriv kmj = k - j fkmj = real(kmj,wp) do jj=1,kmj ihi = i + jj ihmkmj = ihi - kmj work(jj) = (work(jj+1)-work(jj))/(t(ihi)-t(ihmkmj))*fkmj enddo enddo endif ! 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 kpk = k + k j1 = k + 1 j2 = kpk + 1 do j=1,kmider ipj = i + j work(j1) = t(ipj) - x ip1mj = ip1 - j work(j2) = x - t(ip1mj) j1 = j1 + 1 j2 = j2 + 1 enddo iderp1 = ideriv + 1 do j=iderp1,km1 kmj = k - j ilo = kmj do jj=1,kmj work(jj) = (work(jj+1)*work(kpk+ilo)+work(jj)*& work(k+jj))/(work(kpk+ilo)+work(k+jj)) ilo = ilo - 1 enddo enddo endif iflag = 0 dbvalu = work(1) endfunction dbvalu