commit 73292653e364a5af76ccb2ddd57facc492fb4773
Author: Rhys Ulerich <rhys.ulerich@gmail.com>
Date:   Wed Nov 19 15:15:14 2008 -0600

    Adding derivative capabilities to gsl bsplines

diff --git a/bspline/TODO b/bspline/TODO
index 67f8aa8..41e9bf0 100644
--- a/bspline/TODO
+++ b/bspline/TODO
@@ -1,7 +1,5 @@
 Add functions:
 
-gsl_bspline_eval_deriv (see pppack bsplvd) to compute B-spline derivatives
-
 gsl_bspline_smooth to fit smoothing splines to data more efficiently
 than the standard least squares inversion (see pppack l2appr and
 smooth.spline() from GNU R)
diff --git a/bspline/bspline.c b/bspline/bspline.c
index d38e855..f7bd978 100644
--- a/bspline/bspline.c
+++ b/bspline/bspline.c
@@ -1,17 +1,18 @@
 /* bspline/bspline.c
- * 
+ *
  * Copyright (C) 2006, 2007, 2008 Patrick Alken
- * 
+ * Copyright (C) 2008 Rhys Ulerich
+ *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 3 of the License, or (at
  * your option) any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
@@ -27,18 +28,47 @@
  *
  * [1] Carl de Boor, "A Practical Guide to Splines", Springer
  *     Verlag, 1978.
+ *
+ * The bspline_pppack_* internal routines contain code adapted from
+ *
+ * [2] "PPPACK - Piecewise Polynomial Package",
+ *     http://www.netlib.org/pppack/
+ *
  */
 
-static int bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
-                            gsl_bspline_workspace *w);
 
 static inline size_t bspline_find_interval(const double x, int *flag,
-                                           gsl_bspline_workspace *w);
+    gsl_bspline_workspace *w);
+
+static inline int bspline_process_interval_for_eval(
+  const double x, size_t *i, int flag, gsl_bspline_workspace *w);
+
+static void
+bspline_pppack_bsplvb(const gsl_vector *t,
+                      const size_t jhigh,
+                      const size_t index,
+                      const double x,
+                      const size_t left,
+                      size_t *j,
+                      gsl_vector *deltal,
+                      gsl_vector *deltar,
+                      gsl_vector *biatx);
+
+static void
+bspline_pppack_bsplvd(const gsl_vector *t,
+                      const size_t k,
+                      const double x,
+                      const size_t left,
+                      gsl_vector *deltal,
+                      gsl_vector *deltar,
+                      gsl_matrix *a,
+                      gsl_matrix *dbiatx,
+                      const size_t nderiv);
 
 /*
 gsl_bspline_alloc()
   Allocate space for a bspline workspace. The size of the
-workspace is O(5k + nbreak)
+workspace is O(5k + nbreak + 2k^2)
 
 Inputs: k      - spline order (cubic = 4)
         nbreak - number of breakpoints
@@ -89,6 +119,13 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak)
           GSL_ERROR_NULL("failed to allocate space for delta vectors", GSL_ENOMEM);
         }
 
+      w->A = gsl_matrix_alloc(k, k);
+      if (w->A == 0)
+        {
+          gsl_bspline_free(w);
+          GSL_ERROR_NULL("failed to allocate space for derivative work matrix", GSL_ENOMEM);
+        }
+
       w->B = gsl_vector_alloc(k);
       if (w->B == 0)
         {
@@ -96,6 +133,14 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak)
           GSL_ERROR_NULL("failed to allocate space for temporary spline vector", GSL_ENOMEM);
         }
 
+      w->dB = gsl_matrix_alloc(k, k + 1);
+      if (w->dB == 0)
+        {
+          gsl_bspline_free(w);
+          GSL_ERROR_NULL("failed to allocate space for temporary derivative matrix", GSL_ENOMEM);
+        }
+
+
       return (w);
     }
 } /* gsl_bspline_alloc() */
@@ -149,9 +194,15 @@ gsl_bspline_free(gsl_bspline_workspace *w)
   if (w->deltar)
     gsl_vector_free(w->deltar);
 
+  if (w->A)
+    gsl_matrix_free(w->A);
+
   if (w->B)
     gsl_vector_free(w->B);
 
+  if (w->dB)
+    gsl_matrix_free(w->dB);
+
   free(w);
 } /* gsl_bspline_free() */
 
@@ -256,8 +307,8 @@ gsl_bspline_knots_uniform(const double a, const double b,
 /*
 gsl_bspline_eval()
   Evaluate the basis functions B_i(x) for all i. This is
-basically a wrapper function for bspline_eval_all()
-which formats the output in a nice way.
+a wrapper function for gsl_bspline_eval_nonzero() which
+formats the output in a nice way.
 
 Inputs: x - point for evaluation
         B - (output) where to store B_i(x) values
@@ -277,169 +328,266 @@ gsl_bspline_eval(const double x, gsl_vector *B,
 {
   if (B->size != w->n)
     {
-      GSL_ERROR("B vector length does not match n", GSL_EBADLEN);
+      GSL_ERROR("vector B not of length n", GSL_EBADLEN);
     }
   else
     {
-      size_t i;     /* looping */
-      size_t idx = 0;
-      size_t start; /* first non-zero spline */
+      size_t i;      /* looping */
+      size_t istart; /* first non-zero spline for x */
+      size_t iend;   /* last non-zero spline for x, knot for x*/
+      int error;     /* error handling */
 
       /* find all non-zero B_i(x) values */
-      bspline_eval_all(x, w->B, &idx, w);
+      error = gsl_bspline_eval_nonzero(x, w->B, &istart, &iend, w);
+      if (error)
+        {
+          return error;
+        }
 
       /* store values in appropriate part of given vector */
-
-      start = idx - w->k + 1;
-      for (i = 0; i < start; ++i)
+      for (i = 0; i < istart; ++i)
         gsl_vector_set(B, i, 0.0);
 
-      for (i = start; i <= idx; ++i)
-        gsl_vector_set(B, i, gsl_vector_get(w->B, i - start));
+      for (i = istart; i <= iend; ++i)
+        gsl_vector_set(B, i, gsl_vector_get(w->B, i - istart));
 
-      for (i = idx + 1; i < w->n; ++i)
+      for (i = iend + 1; i < w->n; ++i)
         gsl_vector_set(B, i, 0.0);
 
       return GSL_SUCCESS;
     }
 } /* gsl_bspline_eval() */
 
-/****************************************
- *          INTERNAL ROUTINES           *
- ****************************************/
 
 /*
-bspline_eval_all()
-  Evaluate all non-zero B-splines B_i(x) using algorithm (8)
-of chapter X of [1]
-
-The idea is something like this. Given x \in [t_i, t_{i+1}]
-with t_i < t_{i+1} and the t_i are knots, the values of the
-B-splines not automatically zero fit into a triangular
-array as follows:
-                         0
-            0
-0                        B_{i-2,3}
-            B_{i-1,2}
-B_{i,1}                  B_{i-1,3}       ...
-            B_{i,2}
-0                        B_{i,3}
-            0
-                         0
-
-where B_{i,k} is the ith B-spline of order k. The boundary 0s
-indicate that those B-splines not in the table vanish at x.
-
-To compute the non-zero B-splines of a given order k, we use
-Eqs. (4) and (5) of chapter X of [1]:
-
-(4) B_{i,1}(x) = { 1, t_i <= x < t_{i+1}
-                   0, else }
-
-(5) B_{i,k}(x) =     (x - t_i)
-                 ----------------- B_{i,k-1}(x)
-                 (t_{i+k-1} - t_i)
-
-                      t_{i+k} - x
-                 + ----------------- B_{i+1,k-1}(x)
-                   t_{i+k} - t_{i+1}
-
-So (4) gives us the first column of the table and we can use
-the recurrence relation (5) to get the rest of the columns.
-
-Inputs: x   - point at which to evaluate splines
-        B   - (output) where to store B-spline values (length k)
-        idx - (output) B-spline function index of last output
-              value (B_{idx}(x) is stored in the last slot of 'B')
-        w   - bspline workspace
+gsl_bspline_eval_nonzero()
+  Evaluate all non-zero B-spline functions at point x.
+These are the B_i(x) for i in [istart, iend].
+Always B_i(x) = 0 for i < istart and for i > iend.
+
+Inputs: x      - point at which to evaluate splines
+        Bk     - (output) where to store B-spline values (length k)
+        istart - (output) B-spline function index of
+                    first non-zero basis for given x
+        iend   - (output) B-spline function index of
+                 last non-zero basis for given x.
+                 This is also the knot index corresponding to x.
+        w      - bspline workspace
 
 Return: success or error
 
 Notes: 1) the w->knots vector must be initialized before calling
           this function
 
-       2) On output, B contains:
-
-       B = [B_{i-k+1,k}, B_{i-k+2,k}, ..., B_{i-1,k}, B_{i,k}]
+       2) On output, B contains
 
-          where 'i' is stored in 'idx' on output
+              [B_{istart,k}, B_{istart+1,k},
+               ..., B_{iend-1,k}, B_{iend,k}]
 
-       3) based on PPPACK bsplvb
+          evaluated at the given x.
 */
 
-static int
-bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
-                 gsl_bspline_workspace *w)
+int
+gsl_bspline_eval_nonzero(const double x,
+                         gsl_vector *Bk,
+                         size_t *istart,
+                         size_t *iend,
+                         gsl_bspline_workspace *w)
 {
-  if (B->size != w->k)
+  if (Bk->size != w->k)
     {
-      GSL_ERROR("B vector not of length k", GSL_EBADLEN);
+      GSL_ERROR("Bk vector length does not match order k", GSL_EBADLEN);
     }
   else
     {
       size_t i;  /* spline index */
       size_t j;  /* looping */
-      size_t ii; /* looping */
-      int flag = 0;  /* error flag */
-      double saved;
-      double term;
+      int flag  = 0; /* interval search flag */
+      int error = 0; /* error flag */
 
       i = bspline_find_interval(x, &flag, w);
-
-      if (flag == -1)
+      error = bspline_process_interval_for_eval(x, &i, flag, w);
+      if (error)
         {
-          GSL_ERROR("x outside of knot interval", GSL_EINVAL);
+          return error;
         }
-      else if (flag == 1)
+
+      *istart = i - w->k + 1;
+      *iend = i;
+
+      bspline_pppack_bsplvb(w->knots, w->k, 1, x, *iend,
+                            &j, w->deltal, w->deltar, Bk);
+
+      return GSL_SUCCESS;
+    }
+} /* gsl_bspline_eval_nonzero() */
+
+/*
+gsl_bspline_eval_deriv()
+  Evaluate d^j/dx^j B_i(x) for all i, 0 <= j <= nderiv.
+This is a wrapper function for gsl_bspline_eval_deriv_nonzero()
+which formats the output in a nice way.
+
+Inputs: x      - point for evaluation
+        nderiv - number of derivatives to compute,
+                 inclusive.
+        dB     - (output) where to store d^j/dx^j B_i(x)
+                 values. the size of this matrix is
+                 (n = nbreak + k - 2 = l + k - 1 = w->n)
+                 by (nderiv + 1)
+        w      - bspline workspace
+
+Return: success or error
+
+Notes: 1) The w->knots vector must be initialized prior to calling
+          this function (see gsl_bspline_knots())
+
+       2) based on PPPACK's bsplvd
+*/
+
+int
+gsl_bspline_eval_deriv(const double x,
+                       const size_t nderiv,
+                       gsl_matrix *dB,
+                       gsl_bspline_workspace *w)
+{
+  if (dB->size1 != w->n)
+    {
+      GSL_ERROR("dB matrix first dimension not of length n", GSL_EBADLEN);
+    }
+  else if (dB->size2 < nderiv+1)
+    {
+      GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN);
+    }
+  else
+    {
+      size_t i;         /* looping */
+      size_t j;         /* looping */
+      size_t istart;    /* first non-zero spline for x */
+      size_t iend;      /* last non-zero spline for x, knot for x*/
+      int error;        /* error handling */
+
+      /* find all non-zero d^j/dx^j B_i(x) values */
+      error = gsl_bspline_eval_deriv_nonzero(
+                x, nderiv, w->dB, &istart, &iend, w);
+      if (error)
         {
-          if (x <= gsl_vector_get(w->knots, i) + GSL_DBL_EPSILON)
-            {
-              --i;
-            }
-          else
-            {
-              GSL_ERROR("x outside of knot interval", GSL_EINVAL);
-            }
+          return error;
         }
 
-      if (gsl_vector_get(w->knots, i) == gsl_vector_get(w->knots, i + 1))
+      /* store values in appropriate part of given matrix */
+      for (j = 0; j <= nderiv; ++j)
         {
-          GSL_ERROR("knot(i) = knot(i+1) will result in division by zero",
-                    GSL_EINVAL);
+          for (i = 0; i < istart; ++i)
+            gsl_matrix_set(dB, i, j, 0.0);
+
+          for (i = istart; i <= iend; ++i)
+            gsl_matrix_set(dB, i, j, gsl_matrix_get(w->dB, i - istart, j));
+
+          for (i = iend + 1; i < w->n; ++i)
+            gsl_matrix_set(dB, i, j, 0.0);
         }
 
-      *idx = i;
+      return GSL_SUCCESS;
+    }
+} /* gsl_bspline_eval_deriv() */
 
-      gsl_vector_set(B, 0, 1.0);
+/*
+gsl_bspline_eval_deriv_nonzero()
+  At point x evaluate all requested, non-zero B-spline function
+derivatives and store them in dB.  These are the
+d^j/dx^j B_i(x) with i in [istart, iend] and j in [0, nderiv].
+Always d^j/dx^j B_i(x) = 0 for i < istart and for i > iend.
+
+Inputs: x      - point at which to evaluate splines
+        nderiv - number of derivatives to request, inclusive
+        dB     - (output) where to store dB-spline derivatives
+                 (size k by nderiv + 1)
+        istart - (output) B-spline function index of
+                 first non-zero basis for given x
+        iend   - (output) B-spline function index of
+                 last non-zero basis for given x.
+                 This is also the knot index corresponding to x.
+        w      - bspline workspace
 
-      for (j = 0; j < w->k - 1; ++j)
-        {
-          gsl_vector_set(w->deltar, j,
-                         gsl_vector_get(w->knots, i + j + 1) - x);
-          gsl_vector_set(w->deltal, j,
-                         x - gsl_vector_get(w->knots, i - j));
+Return: success or error
 
-          saved = 0.0;
+Notes: 1) the w->knots vector must be initialized before calling
+          this function
 
-          for (ii = 0; ii <= j; ++ii)
-            {
-              term = gsl_vector_get(B, ii) /
-                     (gsl_vector_get(w->deltar, ii) +
-                      gsl_vector_get(w->deltal, j - ii));
+       2) On output, dB contains
 
-              gsl_vector_set(B, ii,
-                             saved +
-                             gsl_vector_get(w->deltar, ii) * term);
+       [[B_{istart,  k}, ..., d^nderiv/dx^nderiv B_{istart  ,k}],
+        [B_{istart+1,k}, ..., d^nderiv/dx^nderiv B_{istart+1,k}],
+        ...
+        [B_{iend-1,  k}, ..., d^nderiv/dx^nderiv B_{iend-1,  k}],
+        [B_{iend,    k}, ..., d^nderiv/dx^nderiv B_{iend,    k}]]
 
-              saved = gsl_vector_get(w->deltal, j - ii) * term;
-            }
+        evaluated at x.  B_{istart, k} is stored in dB(0,0).
+        Each additional column contains an additional derivative.
+
+       3) Note that the zero-th column of the result contains the
+          0th derivative, which is simply a function evaluation.
+
+       4) based on PPPACK's bsplvd
+*/
+
+int
+gsl_bspline_eval_deriv_nonzero(const double x,
+                               const size_t nderiv,
+                               gsl_matrix *dB,
+                               size_t *istart,
+                               size_t *iend,
+                               gsl_bspline_workspace *w)
+{
+  if (dB->size1 != w->k)
+    {
+      GSL_ERROR("dB matrix first dimension not of length k", GSL_EBADLEN);
+    }
+  else if (dB->size2 < nderiv+1)
+    {
+      GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN);
+    }
+  else
+    {
+      size_t i;      /* spline index */
+      size_t j;      /* looping */
+      int flag  = 0; /* interval search flag */
+      int error = 0; /* error flag */
+      size_t min_nderivk;
 
-          gsl_vector_set(B, j + 1, saved);
+      i = bspline_find_interval(x, &flag, w);
+      error = bspline_process_interval_for_eval(x, &i, flag, w);
+      if (error)
+        {
+          return error;
+        }
+
+      *istart = i - w->k + 1;
+      *iend = i;
+
+      bspline_pppack_bsplvd(w->knots, w->k, x, *iend,
+                            w->deltal, w->deltar, w->A, dB, nderiv);
+
+      /* An order k b-spline has at most k-1 nonzero derivatives
+         so we need to zero all requested higher order derivatives */
+      min_nderivk = GSL_MIN_INT(nderiv, w->k - 1);
+      for (j = min_nderivk + 1; j <= nderiv; ++j)
+        {
+          for (i = 0; i < w->k; ++i)
+            {
+              gsl_matrix_set(dB, i, j, 0.0);
+            }
         }
 
       return GSL_SUCCESS;
     }
-} /* bspline_eval_all() */
+} /* gsl_bspline_eval_deriv_nonzero() */
+
+
+/****************************************
+ *          INTERNAL ROUTINES           *
+ ****************************************/
 
 /*
 bspline_find_interval()
@@ -503,3 +651,319 @@ bspline_find_interval(const double x, int *flag,
 
   return i;
 } /* bspline_find_interval() */
+
+
+/*
+bspline_process_interval_for_eval()
+    Consumes an x location, left knot from bspline_find_interval, flag
+from bspline_find_interval, and a workspace.  Checks that x lies within
+the splines' knots, enforces some endpoint continuity requirements, and
+avoids divide by zero errors in the underlying bspline_pppack_* functions.
+ */
+static inline int bspline_process_interval_for_eval(
+  const double x, size_t *i, const int flag, gsl_bspline_workspace
+  *w)
+{
+  if (flag == -1)
+    {
+      GSL_ERROR("x outside of knot interval", GSL_EINVAL);
+    }
+  else if (flag == 1)
+    {
+      if (x <= gsl_vector_get(w->knots, *i) + GSL_DBL_EPSILON)
+        {
+          *i -= 1;
+        }
+      else
+        {
+          GSL_ERROR("x outside of knot interval", GSL_EINVAL);
+        }
+    }
+
+  if (gsl_vector_get(w->knots, *i) == gsl_vector_get(w->knots, *i + 1))
+    {
+      GSL_ERROR("knot(i) = knot(i+1) will result in division by zero",
+                GSL_EINVAL);
+    }
+
+  return GSL_SUCCESS;
+} /* bspline_process_interval_for_eval */
+
+
+/********************************************************************
+ * PPPACK ROUTINES
+ *
+ * The routines herein deliberately avoid using the bspline workspace,
+ * choosing instead to pass all work areas explicitly.  This allows
+ * others to more easily adapt these routines to low memory or
+ * parallel scenarios.
+ ********************************************************************/
+
+/*
+bspline_pppack_bsplvb()
+  calculates the value of all possibly nonzero b-splines at x of order
+jout = max( jhigh , (j+1)*(index-1) ) with knot sequence t.
+
+Parameters:
+        t      - knot sequence, of length left + jout , assumed to be
+                 nondecreasing.  assumption t(left).lt.t(left + 1).
+                 division by zero  will result if t(left) = t(left+1)
+        jhigh  -
+        index  - integers which determine the order jout = max(jhigh,
+                 (j+1)*(index-1))  of the b-splines whose values at x
+                 are to be returned.  index  is used to avoid
+                 recalculations when several columns of the triangular
+                 array of b-spline values are needed (e.g., in  bsplpp
+                 or in  bsplvd ).
+                 precisely, if  index = 1 ,
+                   the calculation starts from scratch and the entire
+                   triangular array of b-spline values of orders
+                   1,2,...,jhigh  is generated order by order , i.e.,
+                   column by column .
+                 if  index = 2 ,
+                   only the b-spline values of order j+1, j+2, ..., jout
+                   are generated, the assumption being that biatx, j,
+                   deltal, deltar are, on entry, as they were on exit
+                   at the previous call.
+                   in particular, if jhigh = 0, then jout = j+1, i.e.,
+                   just the next column of b-spline values is generated.
+        x      - the point at which the b-splines are to be evaluated.
+        left   - an integer chosen (usually) so that
+                 t(left) .le. x .le. t(left+1).
+        j      - (output) a working scalar for indexing
+        deltal - (output) a working area which must be of length at
+                 least jout
+        deltar - (output) a working area which must be of length at
+                 least jout
+        biatx  - (output) array of length jout, with  biatx(i)
+                 containing the value at  x  of the polynomial of order
+                 jout which agrees with the b-spline b(left-jout+i,jout,t)
+                 on the interval (t(left), t(left+1)) .
+
+Method:
+      the recurrence relation
+
+                        x - t(i)              t(i+j+1) - x
+        b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x)
+                      t(i+j)-t(i)            t(i+j+1)-t(i+1)
+
+      is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x),
+      ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),...,
+      b(left,j)(x), storing the new values in  biatx  over the old. the
+      facts that
+                b(i,1) = 1  if  t(i) .le. x .lt. t(i+1)
+      and that
+                b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j)
+      are used. the particular organization of the calculations follows
+      algorithm (8) in chapter x of [1].
+
+Notes:
+
+    (1) This is a direct translation of PPPACK's bsplvb routine with
+        j, deltal, deltar rewritten as input parameters and
+        utilizing zero-based indexing.
+
+    (2) This routine contains no error checking.  Please use routines
+        like gsl_bspline_eval().
+*/
+
+static void
+bspline_pppack_bsplvb(const gsl_vector *t,
+                      const size_t jhigh,
+                      const size_t index,
+                      const double x,
+                      const size_t left,
+                      size_t *j,
+                      gsl_vector *deltal,
+                      gsl_vector *deltar,
+                      gsl_vector *biatx)
+{
+  size_t i;   /* looping */
+  double saved;
+  double term;
+
+  if (index == 1)
+    {
+      *j = 0;
+      gsl_vector_set(biatx, 0, 1.0);
+    }
+
+  for (/* NOP */; *j < jhigh - 1; *j += 1)
+    {
+      gsl_vector_set(deltar, *j, gsl_vector_get(t, left + *j + 1) - x);
+      gsl_vector_set(deltal, *j, x - gsl_vector_get(t, left - *j));
+
+      saved = 0.0;
+
+      for (i = 0; i <= *j; ++i)
+        {
+          term = gsl_vector_get(biatx, i) /
+                 (gsl_vector_get(deltar, i) +
+                  gsl_vector_get(deltal, *j - i));
+
+          gsl_vector_set(biatx, i,
+                         saved + gsl_vector_get(deltar, i) * term);
+
+          saved = gsl_vector_get(deltal, *j - i) * term;
+        }
+
+      gsl_vector_set(biatx, *j + 1, saved);
+    }
+
+  return;
+} /* gsl_bspline_pppack_bsplvb */
+
+/*
+bspline_pppack_bsplvd()
+calculates value and derivs of all b-splines which do not vanish at x
+
+Parameters:
+        t      - the knot array, of length left+k (at least)
+        k      - the order of the b-splines to be evaluated
+        x      - the point at which these values are sought
+        left   - an integer indicating the left endpoint of the interval
+                 of interest. the k b-splines whose support contains the
+                 interval (t(left), t(left+1)) are to be considered.
+                 it is assumed that t(left) .lt. t(left+1)
+                 division by zero will result otherwise (in  bsplvb).
+                 also, the output is as advertised only if
+                 t(left) .le. x .le. t(left+1) .
+        deltal - a working area which must be of length at least k
+        deltar - a working area which must be of length at least k
+        a      - an array of order (k,k), to contain b-coeffs of the
+                 derivatives of a certain order of the k b-splines
+                 of interest.
+        dbiatx - an array of order (k,nderiv). its entry (i,m) contains
+                 value of (m)th derivative of (left-k+i)-th b-spline
+                 of order k for knot sequence  t, i=1,...,k,
+                 m=0,...,nderiv.
+        nderiv - an integer indicating that values of b-splines and
+                 their derivatives up to AND INCLUDING the nderiv-th
+                 are asked for. (nderiv is replaced internally by the
+                 integer mhigh in (1,k) closest to it.)
+
+Method:
+  values at x of all the relevant b-splines of order k,k-1,..., k+1-nderiv
+  are generated via bsplvb and stored temporarily in dbiatx.  then, the
+  b-coeffs of the required derivatives of the b-splines of interest are
+  generated by differencing, each from the preceeding one of lower order,
+  and combined with the values of b-splines of corresponding order in
+  dbiatx  to produce the desired values .
+
+Notes:
+
+    (1) This is a direct translation of PPPACK's bsplvd routine with
+        deltal, deltar rewritten as input parameters (to later feed them
+        to bspline_pppack_bsplvb) and utilizing zero-based indexing.
+
+    (2) This routine contains no error checking.
+*/
+
+static void
+bspline_pppack_bsplvd(const gsl_vector *t,
+                      const size_t k,
+                      const double x,
+                      const size_t left,
+                      gsl_vector *deltal,
+                      gsl_vector *deltar,
+                      gsl_matrix *a,
+                      gsl_matrix *dbiatx,
+                      const size_t nderiv)
+{
+  int i, ideriv, il, j, jlow, jp1mid, kmm, ldummy, m, mhigh;
+  double factor, fkmm, sum;
+
+  size_t bsplvb_j;
+  gsl_vector_view dbcol = gsl_matrix_column(dbiatx, 0);
+
+  mhigh = GSL_MIN_INT(nderiv, k - 1);
+  bspline_pppack_bsplvb(t, k-mhigh, 1, x, left,
+                        &bsplvb_j, deltal, deltar, &dbcol.vector);
+  if (mhigh > 0)
+    {
+      /* the first column of dbiatx always contains the b-spline
+         values for the current order. these are stored in column
+         k-current order before bsplvb is called to put values
+         for the next higher order on top of it.  */
+      ideriv = mhigh;
+      for (m = 1; m <= mhigh; ++m)
+        {
+          for (j = ideriv, jp1mid = 0; j < k; ++j, ++jp1mid)
+            {
+              gsl_matrix_set(dbiatx, j, ideriv,
+                             gsl_matrix_get(dbiatx, jp1mid, 0));
+            }
+          --ideriv;
+          bspline_pppack_bsplvb(t, k-ideriv, 2, x, left,
+                                &bsplvb_j, deltal, deltar, &dbcol.vector);
+        }
+
+      /* at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j)
+         for i=j,...,k-1 and j=0,...,mhigh. in particular, the
+         first column of dbiatx is already in final form. to obtain
+         corresponding derivatives of b-splines in subsequent columns,
+         generate their b-repr. by differencing, then evaluate at x. */
+      jlow = 0;
+      for (i = 0; i < k; ++i)
+        {
+          for (j = jlow; j < k; ++j)
+            {
+              gsl_matrix_set(a, j, i, 0.0);
+            }
+          jlow = i;
+          gsl_matrix_set(a, i, i, 1.0);
+        }
+
+      /* at this point, a(.,j) contains the b-coeffs for the j-th of the
+         k b-splines of interest here. */
+      for (m = 1; m <= mhigh; ++m)
+        {
+          kmm = k - m;
+          fkmm = (float) kmm;
+          il = left;
+          i = k - 1;
+
+          /* for j=1,...,k, construct b-coeffs of (m)th  derivative
+             of b-splines from those for preceding derivative by
+             differencing and store again in  a(.,j) . the fact that
+             a(i,j) = 0  for i .lt. j  is used. */
+          for (ldummy = 0; ldummy < kmm; ++ldummy)
+            {
+              factor = fkmm/(
+                         gsl_vector_get(t, il + kmm) - gsl_vector_get(t, il)
+                       );
+              /* the assumption that t(left).lt.t(left+1) makes
+                 denominator in factor nonzero. */
+              for (j = 0; j <=i; ++j)
+                {
+                  gsl_matrix_set(a, i, j, factor*(
+                                   gsl_matrix_get(a, i, j) - gsl_matrix_get(a, i-1, j)));
+                }
+              --il;
+              --i;
+            }
+
+          /* for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
+             stored in dbiatx(.,m) to get value of (m)th  derivative
+             of i-th b-spline (of interest here) at x, and store in
+             dbiatx(i,m). storage of this value over the value of a
+             b-spline of order m there is safe since the remaining
+             b-spline derivatives of the same order do not use this
+             value due to the fact that a(j,i) = 0 for j .lt. i . */
+          for (i = 0; i < k; ++i)
+            {
+              sum = 0;
+              jlow = GSL_MAX_INT(i, m);
+              for (j = jlow; j < k; ++j)
+                {
+                  sum += gsl_matrix_get(a, j, i)
+                         * gsl_matrix_get(dbiatx, j, m);
+                }
+              gsl_matrix_set(dbiatx, i, m, sum);
+            }
+        }
+    }
+
+  return;
+
+} /* bspline_pppack_bsplvd */
diff --git a/bspline/gsl_bspline.h b/bspline/gsl_bspline.h
index 0a1b248..8a2ffc5 100644
--- a/bspline/gsl_bspline.h
+++ b/bspline/gsl_bspline.h
@@ -1,17 +1,18 @@
 /* bspline/gsl_bspline.h
- * 
+ *
  * Copyright (C) 2006 Patrick Alken
- * 
+ * Copyright (C) 2008 Rhys Ulerich
+ *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 3 of the License, or (at
  * your option) any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
@@ -23,6 +24,7 @@
 #include <stdlib.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
 
 #undef __BEGIN_DECLS
 #undef __END_DECLS
@@ -37,18 +39,20 @@
 __BEGIN_DECLS
 
 typedef struct
-{
-  size_t k;      /* spline order */
-  size_t km1;    /* k - 1 (polynomial order) */
-  size_t l;      /* number of polynomial pieces on interval */
-  size_t nbreak; /* number of breakpoints (l + 1) */
-  size_t n;      /* number of bspline basis functions (l + k - 1) */
-
-  gsl_vector *knots;  /* knots vector */
-  gsl_vector *deltal; /* left delta */
-  gsl_vector *deltar; /* right delta */
-  gsl_vector *B;      /* temporary spline results */
-} gsl_bspline_workspace;
+  {
+    size_t k;      /* spline order */
+    size_t km1;    /* k - 1 (polynomial order) */
+    size_t l;      /* number of polynomial pieces on interval */
+    size_t nbreak; /* number of breakpoints (l + 1) */
+    size_t n;      /* number of bspline basis functions (l + k - 1) */
+
+    gsl_vector *knots;  /* knots vector */
+    gsl_vector *deltal; /* left delta */
+    gsl_vector *deltar; /* right delta */
+    gsl_matrix *A;      /* work matrix */
+    gsl_vector *B;      /* temporary spline results */
+    gsl_matrix *dB;     /* temporary derivative results */
+  } gsl_bspline_workspace;
 
 gsl_bspline_workspace *
 gsl_bspline_alloc(const size_t k, const size_t nbreak);
@@ -67,9 +71,31 @@ int gsl_bspline_knots_uniform(const double a, const double b,
                               gsl_bspline_workspace *w);
 
 int
-gsl_bspline_eval(const double x, gsl_vector *B,
+gsl_bspline_eval(const double x,
+                 gsl_vector *B,
                  gsl_bspline_workspace *w);
 
+int
+gsl_bspline_eval_nonzero(const double x,
+                         gsl_vector *Bk,
+                         size_t *istart,
+                         size_t *iend,
+                         gsl_bspline_workspace *w);
+
+int
+gsl_bspline_eval_deriv(const double x,
+                       const size_t nderiv,
+                       gsl_matrix *dB,
+                       gsl_bspline_workspace *w);
+
+int
+gsl_bspline_eval_deriv_nonzero(const double x,
+                               const size_t nderiv,
+                               gsl_matrix *dB,
+                               size_t *istart,
+                               size_t *iend,
+                               gsl_bspline_workspace *w);
+
 __END_DECLS
 
 #endif /* __GSL_BSPLINE_H__ */
diff --git a/bspline/test.c b/bspline/test.c
index bdc9188..393fbf3 100644
--- a/bspline/test.c
+++ b/bspline/test.c
@@ -1,17 +1,18 @@
 /* bspline/test.c
- * 
+ *
  * Copyright (C) 2006, 2007 Brian Gough
- * 
+ * Copyright (C) 2008 Rhys Ulerich
+ *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 3 of the License, or (at
  * your option) any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
@@ -22,53 +23,73 @@
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_bspline.h>
 #include <gsl/gsl_ieee_utils.h>
+#include <gsl/gsl_nan.h>
 
 void
-test_bspline (gsl_bspline_workspace * bw)
+test_bspline(gsl_bspline_workspace * bw)
 {
   gsl_vector *B;
+  gsl_matrix *dB;
   size_t i, j;
   size_t n = 100;
-  size_t ncoeffs = gsl_bspline_ncoeffs (bw);
-  size_t order = gsl_bspline_order (bw);
-  size_t nbreak = gsl_bspline_nbreak (bw);
-  double a = gsl_bspline_breakpoint (0, bw);
-  double b = gsl_bspline_breakpoint (nbreak - 1, bw);
+  size_t ncoeffs = gsl_bspline_ncoeffs(bw);
+  size_t order = gsl_bspline_order(bw);
+  size_t nbreak = gsl_bspline_nbreak(bw);
+  double a = gsl_bspline_breakpoint(0, bw);
+  double b = gsl_bspline_breakpoint(nbreak - 1, bw);
 
-  B = gsl_vector_alloc (ncoeffs);
+  B  = gsl_vector_alloc(ncoeffs);
+  dB = gsl_matrix_alloc(ncoeffs, 1);
 
+  /* Ensure B-splines form a partition of unity */
   for (i = 0; i < n; i++)
     {
       double xi = a + (b - a) * (i / (n - 1.0));
       double sum = 0;
-      gsl_bspline_eval (xi, B, bw);
+      gsl_bspline_eval(xi, B, bw);
 
       for (j = 0; j < ncoeffs; j++)
         {
-          double Bj = gsl_vector_get (B, j);
+          double Bj = gsl_vector_get(B, j);
           int s = (Bj < 0 || Bj > 1);
-          gsl_test (s,
-                    "basis-spline coefficient %u is in range [0,1] for x=%g",
-                    j, xi);
+          gsl_test(s,
+                   "basis-spline coefficient %u is in range [0,1] for x=%g",
+                   j, xi);
           sum += Bj;
         }
 
-      gsl_test_rel (sum, 1.0, order * GSL_DBL_EPSILON,
-                    "basis-spline order %u is normalized for x=%g", order,
-                    xi);
+      gsl_test_rel(sum, 1.0, order * GSL_DBL_EPSILON,
+                   "basis-spline order %u is normalized for x=%g", order,
+                   xi);
     }
 
-  gsl_vector_free (B);
-}
+  /* Ensure B-splines 0th derivatives agree with regular evaluation */
+  for (i = 0; i < n; i++)
+    {
+      double xi = a + (b - a) * (i / (n - 1.0));
+      gsl_bspline_eval(xi, B, bw);
+      gsl_bspline_eval_deriv(xi, 0, dB, bw);
 
+      for (j = 0; j < ncoeffs; j++)
+        {
+          gsl_test_abs(gsl_matrix_get(dB, j, 0), gsl_vector_get(B, j),
+                       GSL_DBL_EPSILON,
+                       "b-spline order %d basis #%d evaluation and 0th derivative consistent for x=%g",
+                       order, j, xi);
+        }
 
+    }
+
+  gsl_vector_free(B);
+  gsl_matrix_free(dB);
+}
 
 int
-main (int argc, char **argv)
+main(int argc, char **argv)
 {
   size_t order, breakpoints, i;
 
-  gsl_ieee_env_setup ();
+  gsl_ieee_env_setup();
 
   argc = 0;                     /* prevent warnings about unused parameters */
   argv = 0;
@@ -78,10 +99,10 @@ main (int argc, char **argv)
       for (breakpoints = 2; breakpoints < 100; breakpoints++)
         {
           double a = -1.23 * order, b = 45.6 * order;
-          gsl_bspline_workspace *bw = gsl_bspline_alloc (order, breakpoints);
-          gsl_bspline_knots_uniform (a, b, bw);
-          test_bspline (bw);
-          gsl_bspline_free (bw);
+          gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
+          gsl_bspline_knots_uniform(a, b, bw);
+          test_bspline(bw);
+          gsl_bspline_free(bw);
         }
     }
 
@@ -91,21 +112,152 @@ main (int argc, char **argv)
       for (breakpoints = 2; breakpoints < 100; breakpoints++)
         {
           double a = -1.23 * order, b = 45.6 * order;
-          gsl_bspline_workspace *bw = gsl_bspline_alloc (order, breakpoints);
-          gsl_vector *k = gsl_vector_alloc (breakpoints);
+          gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
+          gsl_vector *k = gsl_vector_alloc(breakpoints);
           for (i = 0; i < breakpoints; i++)
             {
               double f, x;
-              f = sqrt (i / (breakpoints - 1.0));
+              f = sqrt(i / (breakpoints - 1.0));
               x = (1 - f) * a + f * b;
-              gsl_vector_set (k, i, x);
+              gsl_vector_set(k, i, x);
             };
-          gsl_bspline_knots (k, bw);
-          test_bspline (bw);
-          gsl_vector_free (k);
-          gsl_bspline_free (bw);
+          gsl_bspline_knots(k, bw);
+          test_bspline(bw);
+          gsl_vector_free(k);
+          gsl_bspline_free(bw);
         }
     }
 
-  exit (gsl_test_summary ());
+  /* Spot check known 0th, 1st, 2nd derivative
+     evaluations for a particular k = 2 case.  */
+  {
+    size_t i, j; /* looping */
+
+    const double xloc[4]     =  { 0.0,  1.0,  6.0,  7.0};
+    const double deriv[4][3] =
+    {
+      { -1.0/2.0,  1.0/2.0, 0.0     },
+      { -1.0/2.0,  1.0/2.0, 0.0     },
+      {      0.0, -1.0/5.0, 1.0/5.0 },
+      {      0.0, -1.0/5.0, 1.0/5.0 }
+    };
+
+    gsl_bspline_workspace *bw = gsl_bspline_alloc(2, 3);
+    gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
+                                      gsl_bspline_order(bw) + 1);
+
+    gsl_vector *breakpts = gsl_vector_alloc(3);
+    gsl_vector_set(breakpts, 0, 0.0);
+    gsl_vector_set(breakpts, 1, 2.0);
+    gsl_vector_set(breakpts, 2, 7.0);
+    gsl_bspline_knots(breakpts, bw);
+
+
+    for (i = 0; i < 4; ++i)  /* at each location */
+      {
+        /* Initialize dB with poison to ensure we overwrite it */
+        gsl_matrix_set_all(dB, GSL_NAN);
+
+        gsl_bspline_eval_deriv(xloc[i], gsl_bspline_order(bw), dB, bw);
+        for (j = 0; j < gsl_bspline_ncoeffs(bw) ; ++j)
+          {
+            /* check basis function 1st deriv */
+            gsl_test_abs(gsl_matrix_get(dB, j, 1), deriv[i][j], GSL_DBL_EPSILON,
+                         "b-spline k=%d basis #%d derivative %d at x = %f",
+                         gsl_bspline_order(bw), j, 1, xloc[i]);
+          }
+        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
+          {
+            /* check k order basis function has k-th deriv equal to 0 */
+            gsl_test_abs(gsl_matrix_get(dB, j, gsl_bspline_order(bw)), 0.0,
+                         GSL_DBL_EPSILON,
+                         "b-spline k=%d basis #%d derivative %d at x = %f",
+                         gsl_bspline_order(bw), j, gsl_bspline_order(bw),
+                         xloc[i]);
+          }
+      }
+
+    gsl_matrix_free(dB);
+    gsl_bspline_free(bw);
+    gsl_vector_free(breakpts);
+  }
+
+  /* Spot check known 0th, 1st, 2nd derivative
+     evaluations for a particular k = 3 case.  */
+  {
+    size_t i, j; /* looping */
+    const double xloc[5]     =  { 0.0, 5.0, 9.0, 12.0, 15.0 };
+    const double eval[5][6]  =
+    {
+      { 4./25.,  69./100.,   3./ 20. ,  0.    , 0.   , 0.    },
+      { 0.     ,  4./21. , 143./210. ,  9./70., 0.   , 0.    },
+      { 0.     ,  0.     ,   3./ 10. ,  7./10., 0.   , 0.    },
+      { 0.     ,  0.     ,   0.      ,  3./4. , 1./4., 0.    },
+      { 0.     ,  0.     ,   0.      ,  1./3. , 5./9., 1./9. }
+    };
+    const double deriv[5][6] =
+    {
+      { -4./25.,  3./50.,   1./ 10.,  0.    , 0.    , 0.      },
+      {  0.    , -2./21.,   1./105.,  3./35., 0.    , 0.      },
+      {  0.    ,  0.    ,  -1./5.  ,  1./ 5., 0.    , 0.      },
+      {  0.    ,  0.    ,   0.     , -1./ 6., 1./6. , 0.      },
+      {  0.    ,  0.    ,   0.     , -1./ 9., 1./27., 2./27. }
+    };
+    const double deriv2[5][6] =
+    {
+      { 2./25., -17./150.,   1.0/30.0 ,  0.0     ,  0.     , 0.     },
+      { 0.    ,   1./ 42., -11.0/210.0,  1.0/35.0,  0.     , 0.     },
+      { 0.    ,   0.     ,   1.0/15.0 ,-11.0/90.0,  1./18. , 0.     },
+      { 0.    ,   0.     ,   0.0      ,  1.0/54.0, -7./162., 2./81. },
+      { 0.    ,   0.     ,   0.0      ,  1.0/54.0, -7./162., 2./81. }
+    };
+
+    gsl_bspline_workspace *bw = gsl_bspline_alloc(3, 5);
+
+    gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
+                                      gsl_bspline_order(bw) + 1);
+
+    gsl_vector *breakpts = gsl_vector_alloc(5);
+    gsl_vector_set(breakpts, 0, -3.0);
+    gsl_vector_set(breakpts, 1,  2.0);
+    gsl_vector_set(breakpts, 2,  9.0);
+    gsl_vector_set(breakpts, 3, 12.0);
+    gsl_vector_set(breakpts, 4, 21.0);
+    gsl_bspline_knots(breakpts, bw);
+
+    for (i = 0; i < 5; ++i)  /* at each location */
+      {
+        /* Initialize dB with poison to ensure we overwrite it */
+        gsl_matrix_set_all(dB, GSL_NAN);
+        gsl_bspline_eval_deriv(xloc[i], gsl_bspline_order(bw), dB, bw);
+
+        /* check basis function evaluation */
+        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
+          {
+            gsl_test_abs(gsl_matrix_get(dB, j, 0), eval[i][j], GSL_DBL_EPSILON,
+                         "b-spline k=%d basis #%d derivative %d at x = %f",
+                         gsl_bspline_order(bw), j, 0, xloc[i]);
+          }
+        /* check 1st derivative evaluation */
+        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
+          {
+            gsl_test_abs(gsl_matrix_get(dB, j, 1), deriv[i][j], GSL_DBL_EPSILON,
+                         "b-spline k=%d basis #%d derivative %d at x = %f",
+                         gsl_bspline_order(bw), j, 1, xloc[i]);
+          }
+        /* check 2nd derivative evaluation */
+        for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
+          {
+            gsl_test_abs(gsl_matrix_get(dB, j, 2), deriv2[i][j], GSL_DBL_EPSILON,
+                         "b-spline k=%d basis #%d derivative %d at x = %f",
+                         gsl_bspline_order(bw), j, 2, xloc[i]);
+          }
+      }
+
+    gsl_matrix_free(dB);
+    gsl_bspline_free(bw);
+    gsl_vector_free(breakpts);
+  }
+
+  exit(gsl_test_summary());
 }
diff --git a/doc/bspline.texi b/doc/bspline.texi
index d106fdb..8d70a05 100644
--- a/doc/bspline.texi
+++ b/doc/bspline.texi
@@ -10,6 +10,7 @@ contains prototypes for the bspline functions and related declarations.
 * Initializing the B-splines solver::
 * Constructing the knots vector::
 * Evaluation of B-spline basis functions::
+* Evaluation of B-spline basis function derivatives::
 * Example programs for B-splines::
 * References and Further Reading::
 @end menu
@@ -120,11 +121,51 @@ once than to compute them individually, due to the nature of the
 defining recurrence relation.
 @end deftypefun
 
+@deftypefun int gsl_bspline_eval_nonzero (const double @var{x}, gsl_vector * @var{Bk}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w})
+This function evaluates all potentially nonzero B-spline basis
+functions at the position @var{x} and stores them in @var{Bk}, so
+that the @math{i}th element of @var{Bk} is @math{B_(istart+i)(x)}.
+The last element of @var{Bk} is @math{B_(iend)(x)}.  @var{Bk} must be
+of length @math{k}.  Only returning nonzero basis functions allows users
+to more cheaply perform tasks requiring linear combinations of the basis
+functions, e.g.  when evaluating an interpolated function.
+@end deftypefun
+
 @deftypefun size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * @var{w})
 This function returns the number of B-spline coefficients given by
 @math{n = nbreak + k - 2}.
 @end deftypefun
 
+@node Evaluation of B-spline basis function derivatives
+@section Evaluation of B-spline derivatives
+@cindex basis splines, derivatives
+
+@deftypefun int gsl_bspline_eval_deriv (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, gsl_bspline_workspace * @var{w})
+This function evaluates all B-spline basis function derivatives of orders
+@math{0} through @math{nderiv} (inclusive) at the position @var{x}
+and stores them in @var{dB}.  The @math{(i,j)}th element of @var{dB}
+is @math{d^j/dx^j B_i(x)}.  @var{dB} must be of size @math{n = nbreak +
+k - 2} by @math{nderiv + 1}.  The value @math{n} may also be obtained
+by calling @code{gsl_bspline_ncoeffs}.  Note that function evaluations
+are included as the @math{0}th order derivatives in @var{dB}.
+It is far more efficient to compute all of the basis functions derivatives
+at once than to compute them individually, due to the nature of the
+defining recurrence relation.
+@end deftypefun
+
+@deftypefun int gsl_bspline_eval_deriv_nonzero (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w})
+This function evaluates all potentially nonzero B-spline basis function
+derivatives of orders @math{0} through @math{nderiv} (inclusive) at
+the position @var{x} and stores them in @var{dB}.  The @math{(i,j)}th
+element of @var{dB} is @math{d^j/dx^j B_(istart+i)(x)}.  The last element
+of @var{dB} is @math{B_(iend)(x)}.  @var{dB} must be of size @math{k} by
+at least @math{nderiv + 1}.  Note that function evaluations are included
+as the @math{0}th order derivatives in @var{dB}.  Only returning nonzero
+basis functions allows users to more cheaply perform tasks requiring
+linear combinations of the basis functions, e.g.  when evaluating an
+interpolated function.
+@end deftypefun
+
 @node Example programs for B-splines
 @section Example programs for B-splines
 @cindex basis splines, examples
