| 1 | subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv ) |
|---|
| 2 | c from * a practical guide to splines * by c. de Boor (7 may 92) |
|---|
| 3 | calls bsplvb |
|---|
| 4 | calculates value and deriv.s of all b-splines which do not vanish at x |
|---|
| 5 | c |
|---|
| 6 | c****** i n p u t ****** |
|---|
| 7 | c t the knot array, of length left+k (at least) |
|---|
| 8 | c k the order of the b-splines to be evaluated |
|---|
| 9 | c x the point at which these values are sought |
|---|
| 10 | c left an integer indicating the left endpoint of the interval of |
|---|
| 11 | c interest. the k b-splines whose support contains the interval |
|---|
| 12 | c (t(left), t(left+1)) |
|---|
| 13 | c are to be considered. |
|---|
| 14 | c a s s u m p t i o n - - - it is assumed that |
|---|
| 15 | c t(left) .lt. t(left+1) |
|---|
| 16 | c division by zero will result otherwise (in b s p l v b ). |
|---|
| 17 | c also, the output is as advertised only if |
|---|
| 18 | c t(left) .le. x .le. t(left+1) . |
|---|
| 19 | c nderiv an integer indicating that values of b-splines and their |
|---|
| 20 | c derivatives up to but not including the nderiv-th are asked |
|---|
| 21 | c for. ( nderiv is replaced internally by the integer m h i g h |
|---|
| 22 | c in (1,k) closest to it.) |
|---|
| 23 | c |
|---|
| 24 | c****** w o r k a r e a ****** |
|---|
| 25 | c a an array of order (k,k), to contain b-coeff.s of the derivat- |
|---|
| 26 | c ives of a certain order of the k b-splines of interest. |
|---|
| 27 | c |
|---|
| 28 | c****** o u t p u t ****** |
|---|
| 29 | c dbiatx an array of order (k,nderiv). its entry (i,m) contains |
|---|
| 30 | c value of (m-1)st derivative of (left-k+i)-th b-spline of |
|---|
| 31 | c order k for knot sequence t , i=1,...,k, m=1,...,nderiv. |
|---|
| 32 | c |
|---|
| 33 | c****** m e t h o d ****** |
|---|
| 34 | c values at x of all the relevant b-splines of order k,k-1,..., |
|---|
| 35 | c k+1-nderiv are generated via bsplvb and stored temporarily in |
|---|
| 36 | c dbiatx . then, the b-coeffs of the required derivatives of the b- |
|---|
| 37 | c splines of interest are generated by differencing, each from the pre- |
|---|
| 38 | c ceding one of lower order, and combined with the values of b-splines |
|---|
| 39 | c of corresponding order in dbiatx to produce the desired values . |
|---|
| 40 | c |
|---|
| 41 | integer k,left,nderiv, i,ideriv,il,j,jlow,jp1mid,kp1,kp1mm |
|---|
| 42 | * ,ldummy,m,mhigh |
|---|
| 43 | real a(k,k),dbiatx(k,nderiv),t(1),x, factor,fkp1mm,sum |
|---|
| 44 | mhigh = max0(min0(nderiv,k),1) |
|---|
| 45 | c mhigh is usually equal to nderiv. |
|---|
| 46 | kp1 = k+1 |
|---|
| 47 | call bsplvb(t,kp1-mhigh,1,x,left,dbiatx) |
|---|
| 48 | if (mhigh .eq. 1) go to 99 |
|---|
| 49 | c the first column of dbiatx always contains the b-spline values |
|---|
| 50 | c for the current order. these are stored in column k+1-current |
|---|
| 51 | c order before bsplvb is called to put values for the next |
|---|
| 52 | c higher order on top of it. |
|---|
| 53 | ideriv = mhigh |
|---|
| 54 | do 15 m=2,mhigh |
|---|
| 55 | jp1mid = 1 |
|---|
| 56 | do 11 j=ideriv,k |
|---|
| 57 | dbiatx(j,ideriv) = dbiatx(jp1mid,1) |
|---|
| 58 | 11 jp1mid = jp1mid + 1 |
|---|
| 59 | ideriv = ideriv - 1 |
|---|
| 60 | call bsplvb(t,kp1-ideriv,2,x,left,dbiatx) |
|---|
| 61 | 15 continue |
|---|
| 62 | c |
|---|
| 63 | c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for |
|---|
| 64 | c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the |
|---|
| 65 | c first column of dbiatx is already in final form. to obtain cor- |
|---|
| 66 | c responding derivatives of b-splines in subsequent columns, gene- |
|---|
| 67 | c rate their b-repr. by differencing, then evaluate at x. |
|---|
| 68 | c |
|---|
| 69 | jlow = 1 |
|---|
| 70 | do 20 i=1,k |
|---|
| 71 | do 19 j=jlow,k |
|---|
| 72 | 19 a(j,i) = 0. |
|---|
| 73 | jlow = i |
|---|
| 74 | 20 a(i,i) = 1. |
|---|
| 75 | c at this point, a(.,j) contains the b-coeffs for the j-th of the |
|---|
| 76 | c k b-splines of interest here. |
|---|
| 77 | c |
|---|
| 78 | do 40 m=2,mhigh |
|---|
| 79 | kp1mm = kp1 - m |
|---|
| 80 | fkp1mm = float(kp1mm) |
|---|
| 81 | il = left |
|---|
| 82 | i = k |
|---|
| 83 | c |
|---|
| 84 | c for j=1,...,k, construct b-coeffs of (m-1)st derivative of |
|---|
| 85 | c b-splines from those for preceding derivative by differencing |
|---|
| 86 | c and store again in a(.,j) . the fact that a(i,j) = 0 for |
|---|
| 87 | c i .lt. j is used. |
|---|
| 88 | do 25 ldummy=1,kp1mm |
|---|
| 89 | factor = fkp1mm/(t(il+kp1mm) - t(il)) |
|---|
| 90 | c the assumption that t(left).lt.t(left+1) makes denominator |
|---|
| 91 | c in factor nonzero. |
|---|
| 92 | do 24 j=1,i |
|---|
| 93 | 24 a(i,j) = (a(i,j) - a(i-1,j))*factor |
|---|
| 94 | il = il - 1 |
|---|
| 95 | 25 i = i - 1 |
|---|
| 96 | c |
|---|
| 97 | c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values |
|---|
| 98 | c stored in dbiatx(.,m) to get value of (m-1)st derivative of |
|---|
| 99 | c i-th b-spline (of interest here) at x , and store in |
|---|
| 100 | c dbiatx(i,m). storage of this value over the value of a b-spline |
|---|
| 101 | c of order m there is safe since the remaining b-spline derivat- |
|---|
| 102 | c ives of the same order do not use this value due to the fact |
|---|
| 103 | c that a(j,i) = 0 for j .lt. i . |
|---|
| 104 | do 40 i=1,k |
|---|
| 105 | sum = 0. |
|---|
| 106 | jlow = max0(i,m) |
|---|
| 107 | do 35 j=jlow,k |
|---|
| 108 | 35 sum = a(j,i)*dbiatx(j,m) + sum |
|---|
| 109 | 40 dbiatx(i,m) = sum |
|---|
| 110 | 99 return |
|---|
| 111 | end |
|---|
