| 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 ) . |
|---|
| 4 | c |
|---|
| 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 |
|---|
| 9 | c to be determined. |
|---|
| 10 | c |
|---|
| 11 | c****** o u t p u t ****** |
|---|
| 12 | c left, mflag.....both integers, whose value is |
|---|
| 13 | c |
|---|
| 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 |
|---|
| 18 | c |
|---|
| 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) . |
|---|
| 26 | c |
|---|
| 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. |
|---|
| 42 | c |
|---|
| 43 | integer left,lxt,mflag, ihi,ilo,istep,middle |
|---|
| 44 | real x,xt(lxt) |
|---|
| 45 | data ilo /1/ |
|---|
| 46 | save ilo |
|---|
| 47 | ihi = ilo + 1 |
|---|
| 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 |
|---|
| 51 | ilo = lxt - 1 |
|---|
| 52 | ihi = lxt |
|---|
| 53 | c |
|---|
| 54 | 20 if (x .ge. xt(ihi)) go to 40 |
|---|
| 55 | if (x .ge. xt(ilo)) go to 100 |
|---|
| 56 | c |
|---|
| 57 | c **** now x .lt. xt(ilo) . decrease ilo to capture x . |
|---|
| 58 | istep = 1 |
|---|
| 59 | 31 ihi = ilo |
|---|
| 60 | ilo = ihi - istep |
|---|
| 61 | if (ilo .le. 1) go to 35 |
|---|
| 62 | if (x .ge. xt(ilo)) go to 50 |
|---|
| 63 | istep = istep*2 |
|---|
| 64 | go to 31 |
|---|
| 65 | 35 ilo = 1 |
|---|
| 66 | if (x .lt. xt(1)) go to 90 |
|---|
| 67 | go to 50 |
|---|
| 68 | c **** now x .ge. xt(ihi) . increase ihi to capture x . |
|---|
| 69 | 40 istep = 1 |
|---|
| 70 | 41 ilo = ihi |
|---|
| 71 | ihi = ilo + istep |
|---|
| 72 | if (ihi .ge. lxt) go to 45 |
|---|
| 73 | if (x .lt. xt(ihi)) go to 50 |
|---|
| 74 | istep = istep*2 |
|---|
| 75 | go to 41 |
|---|
| 76 | 45 if (x .ge. xt(lxt)) go to 110 |
|---|
| 77 | ihi = lxt |
|---|
| 78 | c |
|---|
| 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 |
|---|
| 84 | ilo = middle |
|---|
| 85 | go to 50 |
|---|
| 86 | 53 ihi = middle |
|---|
| 87 | go to 50 |
|---|
| 88 | c**** set output and return. |
|---|
| 89 | 90 mflag = -1 |
|---|
| 90 | left = 1 |
|---|
| 91 | return |
|---|
| 92 | 100 mflag = 0 |
|---|
| 93 | left = ilo |
|---|
| 94 | return |
|---|
| 95 | 110 mflag = 1 |
|---|
| 96 | if (x .eq. xt(lxt)) mflag = 0 |
|---|
| 97 | left = lxt |
|---|
| 98 | 111 if (left .eq. 1) return |
|---|
| 99 | left = left - 1 |
|---|
| 100 | if (xt(left) .lt. xt(lxt)) return |
|---|
| 101 | go to 111 |
|---|
| 102 | end |
|---|
