1 subroutine interv ( xt, lxt, x, left, mflag )
2 c from * a practical guide to splines * by C. de Boor
3 computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) .
5 c****** i n p u t ******
6 c xt.....a real sequence, of length lxt , assumed to be nondecreasing
7 c lxt.....number of terms in the sequence xt .
8 c x.....the point whose location with respect to the sequence xt is
11 c****** o u t p u t ******
12 c left, mflag.....both integers, whose value is
14 c 1 -1 if x .lt. xt(1)
15 c i 0 if xt(i) .le. x .lt. xt(i+1)
16 c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt)
17 c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x
19 c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0
20 c indicates that x lies outside the CLOSED interval
21 c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the
22 c intervals is due to the decision to make all pp functions cont-
23 c inuous from the right, but, by returning mflag = 0 even if
24 C x = xt(lxt), there is the option of having the computed pp function
25 c continuous from the left at xt(lxt) .
27 c****** m e t h o d ******
28 c The program is designed to be efficient in the common situation that
29 c it is called repeatedly, with x taken from an increasing or decrea-
30 c sing sequence. This will happen, e.g., when a pp function is to be
31 c graphed. The first guess for left is therefore taken to be the val-
32 c ue returned at the previous call and stored in the l o c a l varia-
33 c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec-
34 c essary since the present call may have nothing to do with the previ-
35 c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left =
36 c ilo and are done after just three comparisons.
37 c Otherwise, we repeatedly double the difference istep = ihi - ilo
38 c while also moving ilo and ihi in the direction of x , until
39 c xt(ilo) .le. x .lt. xt(ihi) ,
40 c after which we use bisection to get, in addition, ilo+1 = ihi .
41 c left = ilo is then returned.
43 integer left,lxt,mflag, ihi,ilo,istep,middle
44 double precision x,xt(lxt)
48 if (ihi .lt. lxt) go to 20
49 if (x .ge. xt(lxt)) go to 110
50 if (lxt .le. 1) go to 90
54 20 if (x .ge. xt(ihi)) go to 40
55 if (x .ge. xt(ilo)) go to 100
57 c **** now x .lt. xt(ilo) . decrease ilo to capture x .
61 if (ilo .le. 1) go to 35
62 if (x .ge. xt(ilo)) go to 50
66 if (x .lt. xt(1)) go to 90
68 c **** now x .ge. xt(ihi) . increase ihi to capture x .
72 if (ihi .ge. lxt) go to 45
73 if (x .lt. xt(ihi)) go to 50
76 45 if (x .ge. xt(lxt)) go to 110
79 c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
80 50 middle = (ilo + ihi)/2
81 if (middle .eq. ilo) go to 100
82 c note. it is assumed that middle = ilo in case ihi = ilo+1 .
83 if (x .lt. xt(middle)) go to 53
88 c**** set output and return.
96 if (x .eq. xt(lxt)) mflag = 0
98 111 if (left .eq. 1) return
100 if (xt(left) .lt. xt(lxt)) return