bvalue.f

   1       real function bvalue ( t, bcoef, n, k, x, jderiv )
   2 c  from  * a practical guide to splines *  by c. de boor
   3 calls  interv
   4 c
   5 calculates value at  x  of  jderiv-th derivative of spline from b-repr.
   6 c  the spline is taken to be continuous from the right, EXCEPT at the
   7 c  rightmost knot, where it is taken to be continuous from the left.
   8 c
   9 c******  i n p u t ******
  10 c  t, bcoef, n, k......forms the b-representation of the spline  f  to
  11 c        be evaluated. specifically,
  12 c  t.....knot sequence, of length  n+k, assumed nondecreasing.
  13 c  bcoef.....b-coefficient sequence, of length  n .
  14 c  n.....length of  bcoef  and dimension of spline(k,t),
  15 c        a s s u m e d  positive .
  16 c  k.....order of the spline .
  17 c
  18 c  w a r n i n g . . .   the restriction  k .le. kmax (=20)  is imposed
  19 c        arbitrarily by the dimension statement for  aj, dl, dr  below,
  20 c        but is  n o w h e r e  c h e c k e d  for.
  21 c
  22 c  x.....the point at which to evaluate .
  23 c  jderiv.....integer giving the order of the derivative to be evaluated
  24 c        a s s u m e d  to be zero or positive.
  25 c
  26 c******  o u t p u t  ******
  27 c  bvalue.....the value of the (jderiv)-th derivative of  f  at  x .
  28 c
  29 c******  m e t h o d  ******
  30 c     The nontrivial knot interval  (t(i),t(i+1))  containing  x  is lo-
  31 c  cated with the aid of  interv . The  k  b-coeffs of  f  relevant for
  32 c  this interval are then obtained from  bcoef (or taken to be zero if
  33 c  not explicitly available) and are then differenced  jderiv  times to
  34 c  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval.
  35 c  Precisely, with  j = jderiv, we have from x.(12) of the text that
  36 c
  37 c     (d**j)f  =  sum ( bcoef(.,j)*b(.,k-j,t) )
  38 c
  39 c  where
  40 c                   / bcoef(.),                     ,  j .eq. 0
  41 c                   /
  42 c    bcoef(.,j)  =  / bcoef(.,j-1) - bcoef(.-1,j-1)
  43 c                   / ----------------------------- ,  j .gt. 0
  44 c                   /    (t(.+k-j) - t(.))/(k-j)
  45 c
  46 c     Then, we use repeatedly the fact that
  47 c
  48 c    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) )
  49 c  with
  50 c                 (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
  51 c    a(.,x)  =    ---------------------------------------
  52 c                 (x - t(.))      + (t(.+m-1) - x)
  53 c
  54 c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
  55 c  of order  1 , and the coefficient for  b(i,1,t)(x)  must then be the
  56 c  desired number  (d**j)f(x). (see x.(17)-(19) of text).
  57 c
  58       integer jderiv,k,n,   i,ilo,imk,j,jc,jcmin,jcmax,jj,kmax,kmj,km1
  59      *                     ,mflag,nmi,jdrvp1
  60       parameter (kmax = 20)
  61       real bcoef(n),t(n+k),x,   aj(kmax),dl(kmax),dr(kmax),fkmj
  62       bvalue = 0.
  63       if (jderiv .ge. k)                go to 99
  64 c
  65 c  *** Find  i   s.t.   1 .le. i .lt. n+k   and   t(i) .lt. t(i+1)   and
  66 c      t(i) .le. x .lt. t(i+1) . If no such i can be found,  x  lies
  67 c      outside the support of  the spline  f , hence  bvalue = 0.
  68 c      (The asymmetry in this choice of  i  makes  f  rightcontinuous, except
  69 c      at  t(n+k) where it is leftcontinuous.)
  70       call interv ( t, n+k, x, i, mflag )
  71       if (mflag .ne. 0)                 go to 99
  72 c  *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
  73       km1 = k - 1
  74       if (km1 .gt. 0)                   go to 1
  75       bvalue = bcoef(i)
  76                                         go to 99
  77 c
  78 c  *** store the k b-spline coefficients relevant for the knot interval
  79 c     (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
  80 c     dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
  81 c     from input to zero. set any t.s not obtainable equal to t(1) or
  82 c     to t(n+k) appropriately.
  83     1 jcmin = 1
  84       imk = i - k
  85       if (imk .ge. 0)                   go to 8
  86       jcmin = 1 - imk
  87       do 5 j=1,i
  88     5    dl(j) = x - t(i+1-j)
  89       do 6 j=i,km1
  90          aj(k-j) = 0.
  91     6    dl(j) = dl(i)
  92                                         go to 10
  93     8 do 9 j=1,km1
  94     9    dl(j) = x - t(i+1-j)
  95 c
  96    10 jcmax = k
  97       nmi = n - i
  98       if (nmi .ge. 0)                   go to 18
  99       jcmax = k + nmi
 100       do 15 j=1,jcmax
 101    15    dr(j) = t(i+j) - x
 102       do 16 j=jcmax,km1
 103          aj(j+1) = 0.
 104    16    dr(j) = dr(jcmax)
 105                                         go to 20
 106    18 do 19 j=1,km1
 107    19    dr(j) = t(i+j) - x
 108 c
 109    20 do 21 jc=jcmin,jcmax
 110    21    aj(jc) = bcoef(imk + jc)
 111 c
 112 c               *** difference the coefficients  jderiv  times.
 113       if (jderiv .eq. 0)                go to 30
 114       do 23 j=1,jderiv
 115          kmj = k-j
 116          fkmj = float(kmj)
 117          ilo = kmj
 118          do 23 jj=1,kmj
 119             aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
 120    23       ilo = ilo - 1
 121 c
 122 c  *** compute value at  x  in (t(i),t(i+1)) of jderiv-th derivative,
 123 c     given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
 124    30 if (jderiv .eq. km1)              go to 39
 125       jdrvp1 = jderiv + 1
 126       do 33 j=jdrvp1,km1
 127          kmj = k-j
 128          ilo = kmj
 129          do 33 jj=1,kmj
 130             aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
 131    33       ilo = ilo - 1
 132    39 bvalue = aj(1)
 133 c
 134    99                                   return
 135       end