bsplvd.f

   1       subroutine bsplvd ( t, lent, k, x, left, a, dbiatx, nderiv )
   2 c     --------   ------
   3 c      implicit none
   4
   5 C calculates value and deriv.s of all b-splines which do not vanish at x
   6 C calls bsplvb
   7 c
   8 c******  i n p u t  ******
   9 c  t     the knot array, of length left+k (at least)
  10 c  k     the order of the b-splines to be evaluated
  11 c  x     the point at which these values are sought
  12 c  left  an integer indicating the left endpoint of the interval of
  13 c        interest. the  k  b-splines whose support contains the interval
  14 c               (t(left), t(left+1))
  15 c        are to be considered.
  16 c  a s s u m p t i o n  - - -  it is assumed that
  17 c               t(left) < t(left+1)
  18 c        division by zero will result otherwise (in  b s p l v b ).
  19 c        also, the output is as advertised only if
  20 c               t(left) <= x <= t(left+1) .
  21 c  nderiv   an integer indicating that values of b-splines and their
  22 c        derivatives up to but not including the  nderiv-th  are asked
  23 c        for. ( nderiv  is replaced internally by the integer in (1,k)
  24 c        closest to it.)
  25 c
  26 c******  w o r k   a r e a  ******
  27 c  a     an array of order (k,k), to contain b-coeff.s of the derivat-
  28 c        ives of a certain order of the  k  b-splines of interest.
  29 c
  30 c******  o u t p u t  ******
  31 c  dbiatx   an array of order (k,nderiv). its entry  (i,m)  contains
  32 c        value of  (m-1)st  derivative of  (left-k+i)-th  b-spline of
  33 c        order  k  for knot sequence  t , i=m,...,k; m=1,...,nderiv.
  34 c
  35 c******  m e t h o d  ******
  36 c  values at  x  of all the relevant b-splines of order k,k-1,...,
  37 c  k+1-nderiv  are generated via  bsplvb  and stored temporarily
  38 c  in  dbiatx .  then, the b-coeffs of the required derivatives of the
  39 c  b-splines of interest are generated by differencing, each from the
  40 c  preceding one of lower order, and combined with the values of b-
  41 c  splines of corresponding order in  dbiatx  to produce the desired
  42 c  values.
  43
  44 C Args
  45       integer lent,k,left,nderiv
  46       double precision t(lent),x, dbiatx(k,nderiv), a(k,k)
  47 C Locals
  48       double precision factor,fkp1mm,sum
  49       integer i,ideriv,il,j,jlow,jp1mid, kp1,kp1mm,ldummy,m,mhigh
  50
  51       mhigh = max0(min0(nderiv,k),1)
  52 c     mhigh is usually equal to nderiv.
  53       kp1 = k+1
  54       call bsplvb(t,lent,kp1-mhigh,1,x,left,dbiatx)
  55       if (mhigh .eq. 1)                 go to 99
  56 c     the first column of  dbiatx  always contains the b-spline values
  57 c     for the current order. these are stored in column k+1-current
  58 c     order  before  bsplvb  is called to put values for the next
  59 c     higher order on top of it.
  60       ideriv = mhigh
  61       do 15 m=2,mhigh
  62          jp1mid = 1
  63          do 11 j=ideriv,k
  64             dbiatx(j,ideriv) = dbiatx(jp1mid,1)
  65    11       jp1mid = jp1mid + 1
  66          ideriv = ideriv - 1
  67          call bsplvb(t,lent,kp1-ideriv,2,x,left,dbiatx)
  68    15    continue
  69 c
  70 c     at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j) for
  71 c     i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
  72 c     first column of  dbiatx  is already in final form. to obtain cor-
  73 c     responding derivatives of b-splines in subsequent columns, gene-
  74 c     rate their b-repr. by differencing, then evaluate at  x.
  75 c
  76       jlow = 1
  77       do 20 i=1,k
  78          do 19 j=jlow,k
  79    19       a(j,i) = 0d0
  80          jlow = i
  81    20    a(i,i) = 1d0
  82 c     at this point, a(.,j) contains the b-coeffs for the j-th of the
  83 c     k  b-splines of interest here.
  84 c
  85       do 40 m=2,mhigh
  86          kp1mm = kp1 - m
  87          fkp1mm = dble(kp1mm)
  88          il = left
  89          i = k
  90 c
  91 c        for j=1,...,k, construct b-coeffs of  (m-1)st  derivative of
  92 c        b-splines from those for preceding derivative by differencing
  93 c        and store again in  a(.,j) . the fact that  a(i,j) = 0  for
  94 c        i < j  is used.sed.
  95          do 25 ldummy=1,kp1mm
  96             factor = fkp1mm/(t(il+kp1mm) - t(il))
  97 c           the assumption that t(left) < t(left+1) makes denominator
  98 c           in  factor  nonzero.
  99             do 24 j=1,i
 100    24          a(i,j) = (a(i,j) - a(i-1,j))*factor
 101             il = il - 1
 102    25       i = i - 1
 103 c
 104 c        for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
 105 c        stored in dbiatx(.,m) to get value of  (m-1)st  derivative of
 106 c        i-th b-spline (of interest here) at  x , and store in
 107 c        dbiatx(i,m). storage of this value over the value of a b-spline
 108 c        of order m there is safe since the remaining b-spline derivat-
 109 c        ive of the same order do not use this value due to the fact
 110 c        that  a(j,i) = 0  for j < i .
 111          do 40 i=1,k
 112             sum = 0.d0
 113             jlow = max0(i,m)
 114             do 35 j=jlow,k
 115    35          sum = a(j,i)*dbiatx(j,m) + sum
 116    40       dbiatx(i,m) = sum
 117    99 return
 118       end