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