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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
subroutine sproot(t,n,c,zero,mest,m,ier)
c subroutine sproot finds the zeros of a cubic spline s(x),which is
c given in its normalized b-spline representation.
c
c calling sequence:
c call sproot(t,n,c,zero,mest,m,ier)
c
c input parameters:
c t : real array,length n, containing the knots of s(x).
c n : integer, containing the number of knots. n>=8
c c : real array,length n, containing the b-spline coefficients.
c mest : integer, specifying the dimension of array zero.
c
c output parameters:
c zero : real array,lenth mest, containing the zeros of s(x).
c m : integer,giving the number of zeros.
c ier : error flag:
c ier = 0: normal return.
c ier = 1: the number of zeros exceeds mest.
c ier =10: invalid input data (see restrictions).
c
c other subroutines required: fpcuro
c
c restrictions:
c 1) n>= 8.
c 2) t(4) < t(5) < ... < t(n-4) < t(n-3).
c t(1) <= t(2) <= t(3) <= t(4)
c t(n-3) <= t(n-2) <= t(n-1) <= t(n)
c
c author :
c p.dierckx
c dept. computer science, k.u.leuven
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c latest update : march 1987
c
c ..
c ..scalar arguments..
integer n,mest,m,ier
c ..array arguments..
real*8 t(n),c(n),zero(mest)
c ..local scalars..
integer i,j,j1,l,n4
real*8 ah,a0,a1,a2,a3,bh,b0,b1,c1,c2,c3,c4,c5,d4,d5,h1,h2,
* three,two,t1,t2,t3,t4,t5,zz
logical z0,z1,z2,z3,z4,nz0,nz1,nz2,nz3,nz4
c ..local array..
real*8 y(3)
c ..
c set some constants
two = 0.2d+01
three = 0.3d+01
c before starting computations a data check is made. if the input data
c are invalid, control is immediately repassed to the calling program.
n4 = n-4
ier = 10
if(n.lt.8) go to 800
j = n
do 10 i=1,3
if(t(i).gt.t(i+1)) go to 800
if(t(j).lt.t(j-1)) go to 800
j = j-1
10 continue
do 20 i=4,n4
if(t(i).ge.t(i+1)) go to 800
20 continue
c the problem considered reduces to finding the zeros of the cubic
c polynomials pl(x) which define the cubic spline in each knot
c interval t(l)<=x<=t(l+1). a zero of pl(x) is also a zero of s(x) on
c the condition that it belongs to the knot interval.
c the cubic polynomial pl(x) is determined by computing s(t(l)),
c s'(t(l)),s(t(l+1)) and s'(t(l+1)). in fact we only have to compute
c s(t(l+1)) and s'(t(l+1)); because of the continuity conditions of
c splines and their derivatives, the value of s(t(l)) and s'(t(l))
c is already known from the foregoing knot interval.
ier = 0
c evaluate some constants for the first knot interval
h1 = t(4)-t(3)
h2 = t(5)-t(4)
t1 = t(4)-t(2)
t2 = t(5)-t(3)
t3 = t(6)-t(4)
t4 = t(5)-t(2)
t5 = t(6)-t(3)
c calculate a0 = s(t(4)) and ah = s'(t(4)).
c1 = c(1)
c2 = c(2)
c3 = c(3)
c4 = (c2-c1)/t4
c5 = (c3-c2)/t5
d4 = (h2*c1+t1*c2)/t4
d5 = (t3*c2+h1*c3)/t5
a0 = (h2*d4+h1*d5)/t2
ah = three*(h2*c4+h1*c5)/t2
z1 = .true.
if(ah.lt.0.0d0) z1 = .false.
nz1 = .not.z1
m = 0
c main loop for the different knot intervals.
do 300 l=4,n4
c evaluate some constants for the knot interval t(l) <= x <= t(l+1).
h1 = h2
h2 = t(l+2)-t(l+1)
t1 = t2
t2 = t3
t3 = t(l+3)-t(l+1)
t4 = t5
t5 = t(l+3)-t(l)
c find a0 = s(t(l)), ah = s'(t(l)), b0 = s(t(l+1)) and bh = s'(t(l+1)).
c1 = c2
c2 = c3
c3 = c(l)
c4 = c5
c5 = (c3-c2)/t5
d4 = (h2*c1+t1*c2)/t4
d5 = (h1*c3+t3*c2)/t5
b0 = (h2*d4+h1*d5)/t2
bh = three*(h2*c4+h1*c5)/t2
c calculate the coefficients a0,a1,a2 and a3 of the cubic polynomial
c pl(x) = ql(y) = a0+a1*y+a2*y**2+a3*y**3 ; y = (x-t(l))/(t(l+1)-t(l)).
a1 = ah*h1
b1 = bh*h1
a2 = three*(b0-a0)-b1-two*a1
a3 = two*(a0-b0)+b1+a1
c test whether or not pl(x) could have a zero in the range
c t(l) <= x <= t(l+1).
z3 = .true.
if(b1.lt.0.0d0) z3 = .false.
nz3 = .not.z3
if(a0*b0.le.0.0d0) go to 100
z0 = .true.
if(a0.lt.0.0d0) z0 = .false.
nz0 = .not.z0
z2 = .true.
if(a2.lt.0.) z2 = .false.
nz2 = .not.z2
z4 = .true.
if(3.0d0*a3+a2.lt.0.0d0) z4 = .false.
nz4 = .not.z4
if(.not.((z0.and.(nz1.and.(z3.or.z2.and.nz4).or.nz2.and.
* z3.and.z4).or.nz0.and.(z1.and.(nz3.or.nz2.and.z4).or.z2.and.
* nz3.and.nz4))))go to 200
c find the zeros of ql(y).
100 call fpcuro(a3,a2,a1,a0,y,j)
if(j.eq.0) go to 200
c find which zeros of pl(x) are zeros of s(x).
do 150 i=1,j
if(y(i).lt.0.0d0 .or. y(i).gt.1.0d0) go to 150
c test whether the number of zeros of s(x) exceeds mest.
if(m.ge.mest) go to 700
m = m+1
zero(m) = t(l)+h1*y(i)
150 continue
200 a0 = b0
ah = bh
z1 = z3
nz1 = nz3
300 continue
c the zeros of s(x) are arranged in increasing order.
if(m.lt.2) go to 800
do 400 i=2,m
j = i
350 j1 = j-1
if(j1.eq.0) go to 400
if(zero(j).ge.zero(j1)) go to 400
zz = zero(j)
zero(j) = zero(j1)
zero(j1) = zz
j = j1
go to 350
400 continue
j = m
m = 1
do 500 i=2,j
if(zero(i).eq.zero(m)) go to 500
m = m+1
zero(m) = zero(i)
500 continue
go to 800
700 ier = 1
800 return
end
|
