bsplvb.f

   1       subroutine bsplvb ( t, lent,jhigh, index, x, left, biatx )
   2 c      implicit none
   3 c     -------------
   4
   5 calculates the value of all possibly nonzero b-splines at  x  of order
   6 c
   7 c               jout  =  dmax( jhigh , (j+1)*(index-1) )
   8 c
   9 c  with knot sequence  t .
  10 c
  11 c******  i n p u t  ******
  12 c  t.....knot sequence, of length  left + jout  , assumed to be nonde-
  13 c        creasing.
  14 c    a s s u m p t i o n  :  t(left)  <  t(left + 1)
  15 c    d i v i s i o n  b y  z e r o  will result if  t(left) = t(left+1)
  16 c
  17 c  jhigh,
  18 c  index.....integers which determine the order  jout = max(jhigh,
  19 c        (j+1)*(index-1))  of the b-splines whose values at  x  are to
  20 c        be returned.  index  is used to avoid recalculations when seve-
  21 c        ral columns of the triangular array of b-spline values are nee-
  22 c        ded (e.g., in  bvalue  or in  bsplvd ). precisely,
  23 c                     if  index = 1 ,
  24 c        the calculation starts from scratch and the entire triangular
  25 c        array of b-spline values of orders 1,2,...,jhigh  is generated
  26 c        order by order , i.e., column by column .
  27 c                     if  index = 2 ,
  28 c        only the b-spline values of order  j+1, j+2, ..., jout  are ge-
  29 c        nerated, the assumption being that  biatx , j , deltal , deltar
  30 c        are, on entry, as they were on exit at the previous call.
  31 c           in particular, if  jhigh = 0, then  jout = j+1, i.e., just
  32 c        the next column of b-spline values is generated.
  33 c
  34 c  w a r n i n g . . .  the restriction   jout <= jmax (= 20)  is
  35 c        imposed arbitrarily by the dimension statement for  deltal and
  36 c        deltar  below, but is  n o w h e r e  c h e c k e d  for .
  37 c
  38 c  x.....the point at which the b-splines are to be evaluated.
  39 c  left.....an integer chosen (usually) so that
  40 c                  t(left) <= x <= t(left+1)  .
  41 c
  42 c******  o u t p u t  ******
  43 c  biatx.....array of length  jout , with  biatx(i)  containing the val-
  44 c        ue at  x  of the polynomial of order  jout  which agrees with
  45 c        the b-spline  b(left-jout+i,jout,t)  on the interval (t(left),
  46 c        t(left+1)) .
  47 c
  48 c******  m e t h o d  ******
  49 c  the recurrence relation
  50 c
  51 c                       x - t(i)               t(i+j+1) - x
  52 c     b(i,j+1)(x)  =  ----------- b(i,j)(x) + --------------- b(i+1,j)(x)
  53 c                     t(i+j)-t(i)             t(i+j+1)-t(i+1)
  54 c
  55 c  is used (repeatedly) to generate the
  56 c  (j+1)-vector  b(left-j,j+1)(x),...,b(left,j+1)(x)
  57 c  from the j-vector  b(left-j+1,j)(x),...,b(left,j)(x),
  58 c  storing the new values in  biatx  over the old.  the facts that
  59 c            b(i,1) = 1         if  t(i) <= x < t(i+1)
  60 c  and that
  61 c            b(i,j)(x) = 0  unless  t(i) <= x < t(i+j)
  62 c  are used. the particular organization of the calculations follows
  63 c  algorithm (8)  in chapter x of the text.
  64 c
  65
  66 C Arguments
  67       integer lent, jhigh, index, left
  68       double precision t(lent),x, biatx(jhigh)
  69 c     dimension     t(left+jout), biatx(jout)
  70 c     -----------------------------------
  71 c current fortran standard makes it impossible to specify the length of
  72 c  t  and of  biatx  precisely without the introduction of otherwise
  73 c  superfluous additional arguments.
  74
  75 C Local Variables
  76       integer jmax
  77       parameter(jmax = 20)
  78       integer i,j,jp1
  79       double precision deltal(jmax), deltar(jmax),saved,term
  80
  81       save j,deltal,deltar
  82       data j/1/
  83 c
  84                                         go to (10,20), index
  85    10 j = 1
  86       biatx(1) = 1d0
  87       if (j .ge. jhigh)                 go to 99
  88 c
  89    20    jp1 = j + 1
  90          deltar(j) = t(left+j) - x
  91          deltal(j) = x - t(left+1-j)
  92          saved = 0d0
  93          do 26 i=1,j
  94             term = biatx(i)/(deltar(i) + deltal(jp1-i))
  95             biatx(i) = saved + deltar(i)*term
  96    26       saved = deltal(jp1-i)*term
  97          biatx(jp1) = saved
  98          j = jp1
  99          if (j .lt. jhigh)              go to 20
 100 c
 101    99                                   return
 102       end