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
|
subroutine bspp2d ( t, bcoef, n, k, m, scrtch, break, coef, l )
c from * a practical guide to splines * by c. de boor
calls bsplvb
c this is an extended version of bsplpp for use with tensor products
c
converts the b-representation t, bcoef(.,j), n, k of some spline into
c its pp-representation break, coef(j,.,.), l, k , j=1, ..., m .
c
c****** i n p u t ******
c t.....knot sequence, of length n+k
c bcoef(.,j) b-spline coefficient sequence, of length n ,j=1,...,m
c n.....length of bcoef and dimension of spline space spline(k,t)
c k.....order of the spline
c
c w a r n i n g . . . the restriction k .le. kmax (= 20) is impo-
c sed by the arbitrary dimension statement for biatx below, but
c is n o w h e r e c h e c k e d for.
c
c m number of data sets
c
c****** w o r k a r e a ******
c scrtch of size (k,k,m), needed to contain bcoeffs of a piece of
c the spline and its k-1 derivatives for each of the m sets
c
c****** o u t p u t ******
c break.....breakpoint sequence, of length l+1, contains (in increas-
c ing order) the distinct points in the sequence t(k),...,t(n+1)
c coef(mm,.,.) array of size (k,n), with coef(mm,i,j) = (i-1)st der-
c ivative of mm-th spline at break(j) from the right, mm=1,.,m
c l.....number of polynomial pieces which make up the spline in the in-
c terval (t(k), t(n+1))
c
c****** m e t h o d ******
c for each breakpoint interval, the k relevant b-coeffs of the
c spline are found and then differenced repeatedly to get the b-coeffs
c of all the derivatives of the spline on that interval. the spline and
c its first k-1 derivatives are then evaluated at the left end point
c of that interval, using bsplvb repeatedly to obtain the values of
c all b-splines of the appropriate order at that point.
c
c parameter kmax = 20
integer k,l,m,n, i,j,jp1,kmj,left,lsofar,mm
real bcoef(n,m),break(1),coef(m,k,1),t(1), scrtch(k,k,m)
* ,biatx(20),diff,fkmj,sum
c
c * ,biatx(kmax),diff,fkmj,sum
c dimension break(l+1),coef(k,l),t(n+k)
current fortran standard makes it impossible to specify the length of
c break , coef and t precisely without the introduction of otherwise
c superfluous additional arguments.
lsofar = 0
break(1) = t(k)
do 50 left=k,n
c find the next nontrivial knot interval.
if (t(left+1) .eq. t(left)) go to 50
lsofar = lsofar + 1
break(lsofar+1) = t(left+1)
if (k .gt. 1) go to 9
do 5 mm=1,m
5 coef(mm,1,lsofar) = bcoef(left,mm)
go to 50
c store the k b-spline coeff.s relevant to current knot interval
c in scrtch(.,1) .
9 do 10 i=1,k
do 10 mm=1,m
10 scrtch(i,1,mm) = bcoef(left-k+i,mm)
c
c for j=1,...,k-1, compute the k-j b-spline coeff.s relevant to
c current knot interval for the j-th derivative by differencing
c those for the (j-1)st derivative, and store in scrtch(.,j+1) .
do 20 jp1=2,k
j = jp1 - 1
kmj = k - j
fkmj = float(kmj)
do 20 i=1,kmj
diff = (t(left+i) - t(left+i - kmj))/fkmj
if (diff .le. 0.) go to 20
do 15 mm=1,m
15 scrtch(i,jp1,mm) =
* (scrtch(i+1,j,mm) - scrtch(i,j,mm))/diff
20 continue
c
c for j = 0, ..., k-1, find the values at t(left) of the j+1
c b-splines of order j+1 whose support contains the current
c knot interval from those of order j (in biatx ), then comb-
c ine with the b-spline coeff.s (in scrtch(.,k-j) ) found earlier
c to compute the (k-j-1)st derivative at t(left) of the given
c spline.
c note. if the repeated calls to bsplvb are thought to gene-
c rate too much overhead, then replace the first call by
c biatx(1) = 1.
c and the subsequent call by the statement
c j = jp1 - 1
c followed by a direct copy of the lines
c deltar(j) = t(left+j) - x
c ......
c biatx(j+1) = saved
c from bsplvb . deltal(kmax) and deltar(kmax) would have to
c appear in a dimension statement, of course.
c
call bsplvb ( t, 1, 1, t(left), left, biatx )
do 25 mm=1,m
25 coef(mm,k,lsofar) = scrtch(1,k,mm)
do 30 jp1=2,k
call bsplvb ( t, jp1, 2, t(left), left, biatx )
kmj = k+1 - jp1
do 30 mm=1,m
sum = 0.
do 28 i=1,jp1
28 sum = biatx(i)*scrtch(i,kmj,mm) + sum
30 coef(mm,kmj,lsofar) = sum
50 continue
l = lsofar
return
end
|
