1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
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
|
