interv.f

   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       double precision 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