1 subroutine bsplvd ( t, lent, k, x, left, a, dbiatx, nderiv ) 2c -------- ------ 3c implicit none 4 5C calculates value and deriv.s of all b-splines which do not vanish at x 6C calls bsplvb 7c 8c****** i n p u t ****** 9c t the knot array, of length left+k (at least) 10c k the order of the b-splines to be evaluated 11c x the point at which these values are sought 12c left an integer indicating the left endpoint of the interval of 13c interest. the k b-splines whose support contains the interval 14c (t(left), t(left+1)) 15c are to be considered. 16c a s s u m p t i o n - - - it is assumed that 17c t(left) < t(left+1) 18c division by zero will result otherwise (in b s p l v b ). 19c also, the output is as advertised only if 20c t(left) <= x <= t(left+1) . 21c nderiv an integer indicating that values of b-splines and their 22c derivatives up to but not including the nderiv-th are asked 23c for. ( nderiv is replaced internally by the integer in (1,k) 24c closest to it.) 25c 26c****** w o r k a r e a ****** 27c a an array of order (k,k), to contain b-coeff.s of the derivat- 28c ives of a certain order of the k b-splines of interest. 29c 30c****** o u t p u t ****** 31c dbiatx an array of order (k,nderiv). its entry (i,m) contains 32c value of (m-1)st derivative of (left-k+i)-th b-spline of 33c order k for knot sequence t , i=m,...,k; m=1,...,nderiv. 34c 35c****** m e t h o d ****** 36c values at x of all the relevant b-splines of order k,k-1,..., 37c k+1-nderiv are generated via bsplvb and stored temporarily 38c in dbiatx . then, the b-coeffs of the required derivatives of the 39c b-splines of interest are generated by differencing, each from the 40c preceding one of lower order, and combined with the values of b- 41c splines of corresponding order in dbiatx to produce the desired 42c values. 43 44C Args 45 integer lent,k,left,nderiv 46 double precision t(lent),x, dbiatx(k,nderiv), a(k,k) 47C Locals 48 double precision factor,fkp1mm,sum 49 integer i,ideriv,il,j,jlow,jp1mid, kp1,kp1mm,ldummy,m,mhigh 50 51 mhigh = max0(min0(nderiv,k),1) 52c mhigh is usually equal to nderiv. 53 kp1 = k+1 54 call bsplvb(t,lent,kp1-mhigh,1,x,left,dbiatx) 55 if (mhigh .eq. 1) return 56c the first column of dbiatx always contains the b-spline values 57c for the current order. these are stored in column k+1-current 58c order before bsplvb is called to put values for the next 59c higher order on top of it. 60 ideriv = mhigh 61 do 15 m=2,mhigh 62 jp1mid = 1 63 do 11 j=ideriv,k 64 dbiatx(j,ideriv) = dbiatx(jp1mid,1) 65 jp1mid = jp1mid + 1 66 11 continue 67 ideriv = ideriv - 1 68 call bsplvb(t,lent,kp1-ideriv,2,x,left,dbiatx) 69 15 continue 70c 71c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for 72c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the 73c first column of dbiatx is already in final form. to obtain cor- 74c responding derivatives of b-splines in subsequent columns, gene- 75c rate their b-repr. by differencing, then evaluate at x. 76c 77 jlow = 1 78 do 20 i=1,k 79 do 19 j=jlow,k 80 a(j,i) = 0d0 81 19 continue 82 jlow = i 83 a(i,i) = 1d0 84 20 continue 85c at this point, a(.,j) contains the b-coeffs for the j-th of the 86c k b-splines of interest here. 87c 88 do 45 m=2,mhigh 89 kp1mm = kp1 - m 90 fkp1mm = dble(kp1mm) 91 il = left 92 i = k 93c 94c for j=1,...,k, construct b-coeffs of (m-1)st derivative of 95c b-splines from those for preceding derivative by differencing 96c and store again in a(.,j) . the fact that a(i,j) = 0 for 97c i < j is used.sed. 98 do 25 ldummy=1,kp1mm 99 factor = fkp1mm/(t(il+kp1mm) - t(il)) 100c the assumption that t(left) < t(left+1) makes denominator 101c in factor nonzero. 102 do 24 j=1,i 103 a(i,j) = (a(i,j) - a(i-1,j))*factor 104 24 continue 105 il = il - 1 106 i = i - 1 107 25 continue 108c 109c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values 110c stored in dbiatx(.,m) to get value of (m-1)st derivative of 111c i-th b-spline (of interest here) at x , and store in 112c dbiatx(i,m). storage of this value over the value of a b-spline 113c of order m there is safe since the remaining b-spline derivat- 114c ive of the same order do not use this value due to the fact 115c that a(j,i) = 0 for j < i . 116 do 40 i=1,k 117 sum = 0.d0 118 jlow = max0(i,m) 119 do 35 j=jlow,k 120 sum = a(j,i)*dbiatx(j,m) + sum 121 35 continue 122 dbiatx(i,m) = sum 123 40 continue 124 45 continue 125 end 126 127 subroutine bsplvb ( t, lent,jhigh, index, x, left, biatx ) 128c implicit none 129c ------------- 130 131calculates the value of all possibly nonzero b-splines at x of order 132c 133c jout = dmax( jhigh , (j+1)*(index-1) ) 134c 135c with knot sequence t . 136c 137c****** i n p u t ****** 138c t.....knot sequence, of length left + jout , assumed to be nonde- 139c creasing. 140c a s s u m p t i o n : t(left) < t(left + 1) 141c d i v i s i o n b y z e r o will result if t(left) = t(left+1) 142c 143c jhigh, 144c index.....integers which determine the order jout = max(jhigh, 145c (j+1)*(index-1)) of the b-splines whose values at x are to 146c be returned. index is used to avoid recalculations when seve- 147c ral columns of the triangular array of b-spline values are nee- 148c ded (e.g., in bvalue or in bsplvd ). precisely, 149c if index = 1 , 150c the calculation starts from scratch and the entire triangular 151c array of b-spline values of orders 1,2,...,jhigh is generated 152c order by order , i.e., column by column . 153c if index = 2 , 154c only the b-spline values of order j+1, j+2, ..., jout are ge- 155c nerated, the assumption being that biatx , j , deltal , deltar 156c are, on entry, as they were on exit at the previous call. 157c in particular, if jhigh = 0, then jout = j+1, i.e., just 158c the next column of b-spline values is generated. 159c 160c w a r n i n g . . . the restriction jout <= jmax (= 20) is 161c imposed arbitrarily by the dimension statement for deltal and 162c deltar below, but is n o w h e r e c h e c k e d for . 163c 164c x.....the point at which the b-splines are to be evaluated. 165c left.....an integer chosen (usually) so that 166c t(left) <= x <= t(left+1) . 167c 168c****** o u t p u t ****** 169c biatx.....array of length jout , with biatx(i) containing the val- 170c ue at x of the polynomial of order jout which agrees with 171c the b-spline b(left-jout+i,jout,t) on the interval (t(left), 172c t(left+1)) . 173c 174c****** m e t h o d ****** 175c the recurrence relation 176c 177c x - t(i) t(i+j+1) - x 178c b(i,j+1)(x) = ----------- b(i,j)(x) + --------------- b(i+1,j)(x) 179c t(i+j)-t(i) t(i+j+1)-t(i+1) 180c 181c is used (repeatedly) to generate the 182c (j+1)-vector b(left-j,j+1)(x),...,b(left,j+1)(x) 183c from the j-vector b(left-j+1,j)(x),...,b(left,j)(x), 184c storing the new values in biatx over the old. the facts that 185c b(i,1) = 1 if t(i) <= x < t(i+1) 186c and that 187c b(i,j)(x) = 0 unless t(i) <= x < t(i+j) 188c are used. the particular organization of the calculations follows 189c algorithm (8) in chapter x of the text. 190c 191 192C Arguments 193 integer lent, jhigh, index, left 194 double precision t(lent),x, biatx(jhigh) 195c dimension t(left+jout), biatx(jout) 196c ----------------------------------- 197c current fortran standard makes it impossible to specify the length of 198c t and of biatx precisely without the introduction of otherwise 199c superfluous additional arguments. 200 201C Local Variables 202 integer jmax 203 parameter(jmax = 20) 204 integer i,j,jp1 205 double precision deltal(jmax), deltar(jmax),saved,term 206 207 save j,deltal,deltar 208 data j/1/ 209c 210c go to (10,20), index 211 if (index .eq. 2) go to 20 212 j = 1 213 biatx(1) = 1d0 214 if (j .ge. jhigh) return 215c 216 20 jp1 = j + 1 217 deltar(j) = t(left+j) - x 218 deltal(j) = x - t(left+1-j) 219 saved = 0d0 220 do 26 i=1,j 221 term = biatx(i)/(deltar(i) + deltal(jp1-i)) 222 biatx(i) = saved + deltar(i)*term 223 saved = deltal(jp1-i)*term 224 26 continue 225 biatx(jp1) = saved 226 j = jp1 227 if (j .lt. jhigh) go to 20 228 end 229