diff options
Diffstat (limited to 'math/deboor/bsplvd.f')
| -rw-r--r-- | math/deboor/bsplvd.f | 111 |
1 files changed, 111 insertions, 0 deletions
diff --git a/math/deboor/bsplvd.f b/math/deboor/bsplvd.f new file mode 100644 index 00000000..d48d0adc --- /dev/null +++ b/math/deboor/bsplvd.f @@ -0,0 +1,111 @@ + subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv ) +c from * a practical guide to splines * by c. de boor +calls bsplvb +calculates value and deriv.s of all b-splines which do not vanish at x +c +c****** i n p u t ****** +c t the knot array, of length left+k (at least) +c k the order of the b-splines to be evaluated +c x the point at which these values are sought +c left an integer indicating the left endpoint of the interval of +c interest. the k b-splines whose support contains the interval +c (t(left), t(left+1)) +c are to be considered. +c a s s u m p t i o n - - - it is assumed that +c t(left) .lt. t(left+1) +c division by zero will result otherwise (in b s p l v b ). +c also, the output is as advertised only if +c t(left) .le. x .le. t(left+1) . +c nderiv an integer indicating that values of b-splines and their +c derivatives up to but not including the nderiv-th are asked +c for. ( nderiv is replaced internally by the integer m h i g h +c in (1,k) closest to it.) +c +c****** w o r k a r e a ****** +c a an array of order (k,k), to contain b-coeff.s of the derivat- +c ives of a certain order of the k b-splines of interest. +c +c****** o u t p u t ****** +c dbiatx an array of order (k,nderiv). its entry (i,m) contains +c value of (m-1)st derivative of (left-k+i)-th b-spline of +c order k for knot sequence t , i=m,...,k, m=1,...,nderiv. +c +c****** m e t h o d ****** +c values at x of all the relevant b-splines of order k,k-1,..., +c k+1-nderiv are generated via bsplvb and stored temporarily in +c dbiatx . then, the b-coeffs of the required derivatives of the b- +c splines of interest are generated by differencing, each from the pre- +c ceding one of lower order, and combined with the values of b-splines +c of corresponding order in dbiatx to produce the desired values . +c + integer k,left,nderiv, i,ideriv,il,j,jlow,jp1mid,kp1,kp1mm + * ,ldummy,m,mhigh + real a(k,k),dbiatx(k,nderiv),t(1),x, factor,fkp1mm,sum + mhigh = max0(min0(nderiv,k),1) +c mhigh is usually equal to nderiv. + kp1 = k+1 + call bsplvb(t,kp1-mhigh,1,x,left,dbiatx) + if (mhigh .eq. 1) go to 99 +c the first column of dbiatx always contains the b-spline values +c for the current order. these are stored in column k+1-current +c order before bsplvb is called to put values for the next +c higher order on top of it. + ideriv = mhigh + do 15 m=2,mhigh + jp1mid = 1 + do 11 j=ideriv,k + dbiatx(j,ideriv) = dbiatx(jp1mid,1) + 11 jp1mid = jp1mid + 1 + ideriv = ideriv - 1 + call bsplvb(t,kp1-ideriv,2,x,left,dbiatx) + 15 continue +c +c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for +c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the +c first column of dbiatx is already in final form. to obtain cor- +c responding derivatives of b-splines in subsequent columns, gene- +c rate their b-repr. by differencing, then evaluate at x. +c + jlow = 1 + do 20 i=1,k + do 19 j=jlow,k + 19 a(j,i) = 0. + jlow = i + 20 a(i,i) = 1. +c at this point, a(.,j) contains the b-coeffs for the j-th of the +c k b-splines of interest here. +c + do 40 m=2,mhigh + kp1mm = kp1 - m + fkp1mm = float(kp1mm) + il = left + i = k +c +c for j=1,...,k, construct b-coeffs of (m-1)st derivative of +c b-splines from those for preceding derivative by differencing +c and store again in a(.,j) . the fact that a(i,j) = 0 for +c i .lt. j is used. + do 25 ldummy=1,kp1mm + factor = fkp1mm/(t(il+kp1mm) - t(il)) +c the assumption that t(left).lt.t(left+1) makes denominator +c in factor nonzero. + do 24 j=1,i + 24 a(i,j) = (a(i,j) - a(i-1,j))*factor + il = il - 1 + 25 i = i - 1 +c +c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values +c stored in dbiatx(.,m) to get value of (m-1)st derivative of +c i-th b-spline (of interest here) at x , and store in +c dbiatx(i,m). storage of this value over the value of a b-spline +c of order m there is safe since the remaining b-spline derivat- +c ives of the same order do not use this value due to the fact +c that a(j,i) = 0 for j .lt. i . + 30 do 40 i=1,k + sum = 0. + jlow = max0(i,m) + do 35 j=jlow,k + 35 sum = a(j,i)*dbiatx(j,m) + sum + 40 dbiatx(i,m) = sum + 99 return + end |