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 t(n),c(n),zero(mest)
c ..local scalars..
integer i,j,j1,l,n4
real 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 y(3)
c ..
c set some constants
two = 0.2e+01
three = 0.3e+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.) 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.) z3 = .false.
nz3 = .not.z3
if(a0*b0.le.0.) go to 100
z0 = .true.
if(a0.lt.0.) z0 = .false.
nz0 = .not.z0
z2 = .true.
if(a2.lt.0.) z2 = .false.
nz2 = .not.z2
z4 = .true.
if(3.0*a3+a2.lt.0.) 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. .or. y(i).gt.1.0) 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
c GH211017
c if(zero(i).eq.zero(m)) go to 500
if(abs(zero(i)-zero(m))<epsilon(zero(i))) go to 500
m = m+1
zero(m) = zero(i)
500 continue
go to 800
700 ier = 1
800 return
end