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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
real function bvalue ( t, bcoef, n, k, x, jderiv )
c from * a practical guide to splines * by c. de boor
calls interv
c
calculates value at x of jderiv-th derivative of spline from b-repr.
c the spline is taken to be continuous from the right.
c
c****** i n p u t ******
c t, bcoef, n, k......forms the b-representation of the spline f to
c be evaluated. specifically,
c t.....knot sequence, of length n+k, assumed nondecreasing.
c bcoef.....b-coefficient sequence, of length n .
c n.....length of bcoef and dimension of spline(k,t),
c a s s u m e d positive .
c k.....order of the spline .
c
c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed
c arbitrarily by the dimension statement for aj, dl, dr below,
c but is n o w h e r e c h e c k e d for.
c
c x.....the point at which to evaluate .
c jderiv.....integer giving the order of the derivative to be evaluated
c a s s u m e d to be zero or positive.
c
c****** o u t p u t ******
c bvalue.....the value of the (jderiv)-th derivative of f at x .
c
c****** m e t h o d ******
c the nontrivial knot interval (t(i),t(i+1)) containing x is lo-
c cated with the aid of interv . the k b-coeffs of f relevant for
c this interval are then obtained from bcoef (or taken to be zero if
c not explicitly available) and are then differenced jderiv times to
c obtain the b-coeffs of (d**jderiv)f relevant for that interval.
c precisely, with j = jderiv, we have from x.(12) of the text that
c
c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
c
c where
c / bcoef(.), , j .eq. 0
c /
c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1)
c / ----------------------------- , j .gt. 0
c / (t(.+k-j) - t(.))/(k-j)
c
c then, we use repeatedly the fact that
c
c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) )
c with
c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
c a(.,x) = ---------------------------------------
c (x - t(.)) + (t(.+m-1) - x)
c
c to write (d**j)f(x) eventually as a linear combination of b-splines
c of order 1 , and the coefficient for b(i,1,t)(x) must then be the
c desired number (d**j)f(x). (see x.(17)-(19) of text).
c
c parameter kmax = 20
integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1,mflag
* ,nmi,jdrvp1
real bcoef(n),t(1),x, aj(20),dl(20),dr(20),fkmj
c real bcoef(n),t(1),x, aj(kmax),dl(kmax),dr(kmax),fkmj
c dimension t(n+k)
current fortran standard makes it impossible to specify the length of t
c precisely without the introduction of otherwise superfluous addition-
c al arguments.
bvalue = 0.
if (jderiv .ge. k) go to 99
c
c *** find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
c t(i) .le. x .lt. t(i+1) . if no such i can be found, x lies
c outside the support of the spline f and bvalue = 0.
c (the asymmetry in this choice of i makes f rightcontinuous)
call interv ( t, n+k, x, i, mflag )
if (mflag .ne. 0) go to 99
c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
km1 = k - 1
if (km1 .gt. 0) go to 1
bvalue = bcoef(i)
go to 99
c
c *** store the k b-spline coefficients relevant for the knot interval
c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
c from input to zero. set any t.s not obtainable equal to t(1) or
c to t(n+k) appropriately.
1 jcmin = 1
imk = i - k
if (imk .ge. 0) go to 8
jcmin = 1 - imk
do 5 j=1,i
5 dl(j) = x - t(i+1-j)
do 6 j=i,km1
aj(k-j) = 0.
6 dl(j) = dl(i)
go to 10
8 do 9 j=1,km1
9 dl(j) = x - t(i+1-j)
c
10 jcmax = k
nmi = n - i
if (nmi .ge. 0) go to 18
jcmax = k + nmi
do 15 j=1,jcmax
15 dr(j) = t(i+j) - x
do 16 j=jcmax,km1
aj(j+1) = 0.
16 dr(j) = dr(jcmax)
go to 20
18 do 19 j=1,km1
19 dr(j) = t(i+j) - x
c
20 do 21 jc=jcmin,jcmax
21 aj(jc) = bcoef(imk + jc)
c
c *** difference the coefficients jderiv times.
if (jderiv .eq. 0) go to 30
do 23 j=1,jderiv
kmj = k-j
fkmj = float(kmj)
ilo = kmj
do 23 jj=1,kmj
aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
23 ilo = ilo - 1
c
c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative,
c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
30 if (jderiv .eq. km1) go to 39
jdrvp1 = jderiv + 1
do 33 j=jdrvp1,km1
kmj = k-j
ilo = kmj
do 33 jj=1,kmj
aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
33 ilo = ilo - 1
39 bvalue = aj(1)
c
99 return
end
|
