interv.f
102 lines
| 4.1 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine interv ( xt, lxt, x, left, mflag ) | |||
| c from * a practical guide to splines * by C. de Boor | ||||
| computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) . | ||||
| c | ||||
| c****** i n p u t ****** | ||||
| c xt.....a real sequence, of length lxt , assumed to be nondecreasing | ||||
| c lxt.....number of terms in the sequence xt . | ||||
| c x.....the point whose location with respect to the sequence xt is | ||||
| c to be determined. | ||||
| c | ||||
| c****** o u t p u t ****** | ||||
| c left, mflag.....both integers, whose value is | ||||
| c | ||||
| c 1 -1 if x .lt. xt(1) | ||||
| c i 0 if xt(i) .le. x .lt. xt(i+1) | ||||
| c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt) | ||||
| c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x | ||||
| c | ||||
| c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0 | ||||
| c indicates that x lies outside the CLOSED interval | ||||
| c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the | ||||
| c intervals is due to the decision to make all pp functions cont- | ||||
| c inuous from the right, but, by returning mflag = 0 even if | ||||
| C x = xt(lxt), there is the option of having the computed pp function | ||||
| c continuous from the left at xt(lxt) . | ||||
| c | ||||
| c****** m e t h o d ****** | ||||
| c The program is designed to be efficient in the common situation that | ||||
| c it is called repeatedly, with x taken from an increasing or decrea- | ||||
| c sing sequence. This will happen, e.g., when a pp function is to be | ||||
| c graphed. The first guess for left is therefore taken to be the val- | ||||
| c ue returned at the previous call and stored in the l o c a l varia- | ||||
| c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec- | ||||
| c essary since the present call may have nothing to do with the previ- | ||||
| c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left = | ||||
| c ilo and are done after just three comparisons. | ||||
| c Otherwise, we repeatedly double the difference istep = ihi - ilo | ||||
| c while also moving ilo and ihi in the direction of x , until | ||||
| c xt(ilo) .le. x .lt. xt(ihi) , | ||||
| c after which we use bisection to get, in addition, ilo+1 = ihi . | ||||
| c left = ilo is then returned. | ||||
| c | ||||
| integer left,lxt,mflag, ihi,ilo,istep,middle | ||||
| real x,xt(lxt) | ||||
| data ilo /1/ | ||||
| save ilo | ||||
| ihi = ilo + 1 | ||||
| if (ihi .lt. lxt) go to 20 | ||||
| if (x .ge. xt(lxt)) go to 110 | ||||
| if (lxt .le. 1) go to 90 | ||||
| ilo = lxt - 1 | ||||
| ihi = lxt | ||||
| c | ||||
| 20 if (x .ge. xt(ihi)) go to 40 | ||||
| if (x .ge. xt(ilo)) go to 100 | ||||
| c | ||||
| c **** now x .lt. xt(ilo) . decrease ilo to capture x . | ||||
| istep = 1 | ||||
| 31 ihi = ilo | ||||
| ilo = ihi - istep | ||||
| if (ilo .le. 1) go to 35 | ||||
| if (x .ge. xt(ilo)) go to 50 | ||||
| istep = istep*2 | ||||
| go to 31 | ||||
| 35 ilo = 1 | ||||
| if (x .lt. xt(1)) go to 90 | ||||
| go to 50 | ||||
| c **** now x .ge. xt(ihi) . increase ihi to capture x . | ||||
| 40 istep = 1 | ||||
| 41 ilo = ihi | ||||
| ihi = ilo + istep | ||||
| if (ihi .ge. lxt) go to 45 | ||||
| if (x .lt. xt(ihi)) go to 50 | ||||
| istep = istep*2 | ||||
| go to 41 | ||||
| 45 if (x .ge. xt(lxt)) go to 110 | ||||
| ihi = lxt | ||||
| c | ||||
| c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval. | ||||
| 50 middle = (ilo + ihi)/2 | ||||
| if (middle .eq. ilo) go to 100 | ||||
| c note. it is assumed that middle = ilo in case ihi = ilo+1 . | ||||
| if (x .lt. xt(middle)) go to 53 | ||||
| ilo = middle | ||||
| go to 50 | ||||
| 53 ihi = middle | ||||
| go to 50 | ||||
| c**** set output and return. | ||||
| 90 mflag = -1 | ||||
| left = 1 | ||||
| return | ||||
| 100 mflag = 0 | ||||
| left = ilo | ||||
| return | ||||
| 110 mflag = 1 | ||||
| if (x .eq. xt(lxt)) mflag = 0 | ||||
| left = lxt | ||||
| 111 if (left .eq. 1) return | ||||
| left = left - 1 | ||||
| if (xt(left) .lt. xt(lxt)) return | ||||
| go to 111 | ||||
| end | ||||
