Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Patrick Alken]

I'm moving this discussion to the gsl-discuss list. Good work!
This is something several people have been asking for and I haven't
had any time to work on.

I will take a look at the patch in the next few days.

Patrick Alken

On Wed, Nov 19, 2008 at 04:12:51PM -0600, Rhys Ulerich wrote:
> Hi all,
> 
> I've just submitted a patch to Savannah to add bspline derivative
> capabilities to GSL.
> Please let me know how I can help get this into the 1.12 release.
> Full details at https://savannah.gnu.org/bugs/index.php?24881
> 
> Thanks,
> Rhys
> _______________________________________________
> Help-gsl mailing list
> Help-gsl@gnu.org
> http://lists.gnu.org/mailman/listinfo/help-gsl

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

At Wed, 19 Nov 2008 16:50:47 -0700,
Patrick Alken wrote:
> 
> I'm moving this discussion to the gsl-discuss list. Good work!
> This is something several people have been asking for and I haven't
> had any time to work on.
> 
> I will take a look at the patch in the next few days.

Patrick, I'm happy for you to take the decision on this wrt the 1.12
release.  The main thing we want to avoid at this stage is breaking
the build in ways that we might not detect in our own testing.

-- 
Brian Gough

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Patrick Alken]

The new code looks clean and very well written. I made a test
program comparing the gsl_bspline_eval_deriv output directly
with PPPACK's bvalue() routine, which is an independent implementation
of bspline derivatives (gsl_bspline_eval_deriv is based on bsplvd()).

The results look fine.  I also ran it through valgrind with good
results. I think it can be put into the next release.

Patrick

On Mon, Nov 24, 2008 at 03:37:57PM +0000, Brian Gough wrote:
> At Wed, 19 Nov 2008 16:50:47 -0700,
> Patrick Alken wrote:
> > 
> > I'm moving this discussion to the gsl-discuss list. Good work!
> > This is something several people have been asking for and I haven't
> > had any time to work on.
> > 
> > I will take a look at the patch in the next few days.
> 
> Patrick, I'm happy for you to take the decision on this wrt the 1.12
> release.  The main thing we want to avoid at this stage is breaking
> the build in ways that we might not detect in our own testing.
> 
> -- 
> Brian Gough

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

At Mon, 24 Nov 2008 10:58:54 -0700,
Patrick Alken wrote:
> 
> The new code looks clean and very well written. I made a test
> program comparing the gsl_bspline_eval_deriv output directly
> with PPPACK's bvalue() routine, which is an independent implementation
> of bspline derivatives (gsl_bspline_eval_deriv is based on bsplvd()).
> 
> The results look fine.  I also ran it through valgrind with good
> results. I think it can be put into the next release.

Sounds good, please go ahead and commit it to savannah.

-- 
Brian Gough

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Patrick Alken]

Ok the new code is in the repository!

On Tue, Nov 25, 2008 at 02:51:07PM +0000, Brian Gough wrote:
> At Mon, 24 Nov 2008 10:58:54 -0700,
> Patrick Alken wrote:
> > 
> > The new code looks clean and very well written. I made a test
> > program comparing the gsl_bspline_eval_deriv output directly
> > with PPPACK's bvalue() routine, which is an independent implementation
> > of bspline derivatives (gsl_bspline_eval_deriv is based on bsplvd()).
> > 
> > The results look fine.  I also ran it through valgrind with good
> > results. I think it can be put into the next release.
> 
> Sounds good, please go ahead and commit it to savannah.
> 
> -- 
> Brian Gough

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Patrick Alken]

Oops, yes it is now in there.

Patrick

On Sat, Nov 29, 2008 at 06:50:20PM -0800, Rhys Ulerich wrote:
>    Hi Patrick,
> 
>    The B-spline derivative commit in Savannah doesn't look like it included
>    the updates to doc/bspline.texi I'd included in the patch.
> 
>    Thank you for accepting this one,
>    Rhys
>    On Sat, Nov 29, 2008 at 1:42 PM, Patrick Alken
>    <[1]patrick.alken@colorado.edu> wrote:
> 
>      Ok the new code is in the repository!
>      On Tue, Nov 25, 2008 at 02:51:07PM +0000, Brian Gough wrote:
>      > At Mon, 24 Nov 2008 10:58:54 -0700,
>      > Patrick Alken wrote:
>      > >
>      > > The new code looks clean and very well written. I made a test
>      > > program comparing the gsl_bspline_eval_deriv output directly
>      > > with PPPACK's bvalue() routine, which is an independent
>      implementation
>      > > of bspline derivatives (gsl_bspline_eval_deriv is based on
>      bsplvd()).
>      > >
>      > > The results look fine.  I also ran it through valgrind with good
>      > > results. I think it can be put into the next release.
>      >
>      > Sounds good, please go ahead and commit it to savannah.
>      >
>      > --
>      > Brian Gough
> 
> References
> 
>    Visible links
>    1. mailto:patrick.alken@colorado.edu

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

I noticed a problem when preparing the new release.  The extra entries
in the bspline derivatives workspace break binary compatibility.

This isn't necessarily a huge problem, but it seems undesirable to
force an upgrade of all applications just for those two struct
entries.  If we're going to break binary compatibility we may as well
do it for all other functions that need it at the same time.

So I'd propose using a separate workspace for the derivative functions
in this release.

See e.g.  branch bspline-compat at http://src.network-theory.com/gsl.git/
which is basically finished, I just need to add some error checking.

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Patrick Alken]

Ah I hadn't considered this issue. If there are other areas of
GSL which have the same issues maybe we should keep a list so
when there is a major new release like 2.0 we could integrate
these "extra" workspaces? Of course one problem with that is
if we have a separate workspace for bspline derivatives it may
never be convenient to integrate it with the normal bspline
workspace.

On Fri, Dec 05, 2008 at 06:56:15PM +0000, Brian Gough wrote:
> I noticed a problem when preparing the new release.  The extra entries
> in the bspline derivatives workspace break binary compatibility.
> 
> This isn't necessarily a huge problem, but it seems undesirable to
> force an upgrade of all applications just for those two struct
> entries.  If we're going to break binary compatibility we may as well
> do it for all other functions that need it at the same time.
> 
> So I'd propose using a separate workspace for the derivative functions
> in this release.
> 
> See e.g.  branch bspline-compat at http://src.network-theory.com/gsl.git/
> which is basically finished, I just need to add some error checking.

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Rhys Ulerich]

[-- Attachment #1: Type: text/plain, Size: 2343 bytes --]

Hi Brian and Patrick,

The attached patch moves the derivative functionality into a separate
workspace.  Users allocate a gsl_bspline_workspace then a
gsl_bspline_deriv_workspace if they need the new capabilities.  All
the derivative-related functions are named like gsl_bspline_deriv_* to
help make the separation clear.  The patch includes updates to
doc/bspline.texi.

Applies cleanly to savannah and Brian's bspline-compat branch.  Brian,
I'm not sure if I missed running a specific compatibility test in your
branch or not...  If you had one in there, I missed it and haven't run
it.

- Rhys

On Fri, Dec 5, 2008 at 4:40 PM, Rhys Ulerich <rhys.ulerich@gmail.com> wrote:
> I'll break out the bspline derivative workspace and get back to you guys
> with a patch.
>
> FYI, the derivative workspace storage is O(order^2) in contrast to the
> evaluation workspace being only O(order).  Leading constant on the
> derivative storage term is small, but I could see where later joining the
> workspaces is a disadvantage for someone needing only evaluation
> capabilities.
>
> - Rhys
>
> On Fri, Dec 5, 2008 at 1:10 PM, Patrick Alken <patrick.alken@colorado.edu>
> wrote:
>>
>> Ah I hadn't considered this issue. If there are other areas of
>> GSL which have the same issues maybe we should keep a list so
>> when there is a major new release like 2.0 we could integrate
>> these "extra" workspaces? Of course one problem with that is
>> if we have a separate workspace for bspline derivatives it may
>> never be convenient to integrate it with the normal bspline
>> workspace.
>>
>> On Fri, Dec 05, 2008 at 06:56:15PM +0000, Brian Gough wrote:
>> > I noticed a problem when preparing the new release.  The extra entries
>> > in the bspline derivatives workspace break binary compatibility.
>> >
>> > This isn't necessarily a huge problem, but it seems undesirable to
>> > force an upgrade of all applications just for those two struct
>> > entries.  If we're going to break binary compatibility we may as well
>> > do it for all other functions that need it at the same time.
>> >
>> > So I'd propose using a separate workspace for the derivative functions
>> > in this release.
>> >
>> > See e.g.  branch bspline-compat at
>> > http://src.network-theory.com/gsl.git/
>> > which is basically finished, I just need to add some error checking.
>
>

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: move-bspline-derivative-data-i.patch --]
[-- Type: text/x-diff; name=move-bspline-derivative-data-i.patch, Size: 48412 bytes --]

Move bspline derivative data into separate workspace.

From: Rhys Ulerich <rhys.ulerich@gmail.com>

Required for binary compatibility with existing code.
---

 bspline/TODO          |    3 
 bspline/bspline.c     |  558 +++++++++++++++++++++++++++----------------------
 bspline/gsl_bspline.h |   66 ++++--
 bspline/test.c        |   20 +-
 doc/bspline.texi      |   21 ++
 5 files changed, 389 insertions(+), 279 deletions(-)

diff --git a/bspline/TODO b/bspline/TODO
index 41e9bf0..37a52c3 100644
--- a/bspline/TODO
+++ b/bspline/TODO
@@ -1,5 +1,8 @@
 Add functions:
 
+add bspline derivatve code to examples/bspline.c or produce a 
+derivative-related example.
+
 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 45bd628..f286707 100644
--- a/bspline/bspline.c
+++ b/bspline/bspline.c
@@ -36,12 +36,12 @@
  *
  */
 
+static inline size_t
+bspline_find_interval(const double x, int *flag, gsl_bspline_workspace *w);
 
-static inline size_t bspline_find_interval(const double x, int *flag,
-    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 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,
@@ -68,7 +68,7 @@ bspline_pppack_bsplvd(const gsl_vector *t,
 /*
 gsl_bspline_alloc()
   Allocate space for a bspline workspace. The size of the
-workspace is O(5k + nbreak + 2k^2)
+workspace is O(5k + nbreak)
 
 Inputs: k      - spline order (cubic = 4)
         nbreak - number of breakpoints
@@ -91,8 +91,7 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak)
     {
       gsl_bspline_workspace *w;
 
-      w = (gsl_bspline_workspace *)
-          calloc(1, sizeof(gsl_bspline_workspace));
+      w = (gsl_bspline_workspace *) calloc(1, sizeof(gsl_bspline_workspace));
       if (w == 0)
         {
           GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);
@@ -108,7 +107,8 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak)
       if (w->knots == 0)
         {
           gsl_bspline_free(w);
-          GSL_ERROR_NULL("failed to allocate space for knots vector", GSL_ENOMEM);
+          GSL_ERROR_NULL("failed to allocate space for knots vector",
+                         GSL_ENOMEM);
         }
 
       w->deltal = gsl_vector_alloc(k);
@@ -116,58 +116,120 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak)
       if (!w->deltal || !w->deltar)
         {
           gsl_bspline_free(w);
-          GSL_ERROR_NULL("failed to allocate space for delta vectors", GSL_ENOMEM);
+          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);
+          GSL_ERROR_NULL("failed to allocate space for derivative work matrix",
+                         GSL_ENOMEM);
         }
 
       w->B = gsl_vector_alloc(k);
       if (w->B == 0)
         {
           gsl_bspline_free(w);
-          GSL_ERROR_NULL("failed to allocate space for temporary spline vector", GSL_ENOMEM);
+          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);
+          GSL_ERROR_NULL(
+            "failed to allocate space for temporary derivative matrix",
+            GSL_ENOMEM);
         }
 
-
       return (w);
     }
 } /* gsl_bspline_alloc() */
 
+/*
+gsl_bspline_deriv_alloc()
+  Allocate space for a bspline derivative workspace. The size of the
+workspace is O(2k^2)
+
+Inputs: bspline - gsl_bspline_workspace of the associated bspline
+
+Return: pointer to workspace
+*/
+
+gsl_bspline_deriv_workspace *
+gsl_bspline_deriv_alloc(gsl_bspline_workspace * bspline)
+{
+  if (bspline == 0)
+    {
+      GSL_ERROR_NULL("bspline must not be null", GSL_EINVAL);
+    }
+  else if (bspline->k == 0)
+    {
+      /* Not strictly necessary, but a sanity check on the bspline */
+      GSL_ERROR_NULL("k must be at least 1", GSL_EINVAL);
+    }
+  else
+    {
+      gsl_bspline_deriv_workspace *w;
+
+      w = (gsl_bspline_deriv_workspace *) calloc(1,
+          sizeof(gsl_bspline_deriv_workspace));
+      if (w == 0)
+        {
+          GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);
+        }
+
+      w->A = gsl_matrix_alloc(bspline->k, bspline->k);
+      if (w->A == 0)
+        {
+          gsl_bspline_deriv_free(w);
+          GSL_ERROR_NULL("failed to allocate space for derivative work matrix",
+                         GSL_ENOMEM);
+        }
+
+      w->dB = gsl_matrix_alloc(bspline->k, bspline->k + 1);
+      if (w->dB == 0)
+        {
+          gsl_bspline_deriv_free(w);
+          GSL_ERROR_NULL(
+            "failed to allocate space for temporary derivative matrix",
+            GSL_ENOMEM);
+        }
+
+      w->bspline = bspline;
+
+      return (w);
+    }
+} /* gsl_bspline_deriv_alloc() */
+
 /* Return number of coefficients */
 size_t
-gsl_bspline_ncoeffs (gsl_bspline_workspace * w)
+gsl_bspline_ncoeffs(gsl_bspline_workspace * w)
 {
   return w->n;
 }
 
 /* Return order */
 size_t
-gsl_bspline_order (gsl_bspline_workspace * w)
+gsl_bspline_order(gsl_bspline_workspace * w)
 {
   return w->k;
 }
 
 /* Return number of breakpoints */
 size_t
-gsl_bspline_nbreak (gsl_bspline_workspace * w)
+gsl_bspline_nbreak(gsl_bspline_workspace * w)
 {
   return w->nbreak;
 }
 
+/* Return the location of the i-th breakpoint*/
 double
-gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
+gsl_bspline_breakpoint(size_t i, gsl_bspline_workspace * w)
 {
   size_t j = i + w->k - 1;
   return gsl_vector_get(w->knots, j);
@@ -175,7 +237,8 @@ gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
 
 /*
 gsl_bspline_free()
-  Free a bspline workspace
+  Free a gsl_bspline_workspace.
+Any associated gsl_bspline_deriv_workspace should be freed beforehand.
 
 Inputs: w - workspace to free
 
@@ -194,28 +257,46 @@ 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);
 
+  free(w);
+} /* gsl_bspline_free() */
+
+/*
+gsl_bspline_deriv_free()
+  Free a gsl_bspline_deriv_workspace.
+The associated gsl_bspline_workspace should be freed afterwards.
+
+Inputs: w - workspace to free
+
+Return: none
+*/
+
+void
+gsl_bspline_deriv_free(gsl_bspline_deriv_workspace *w)
+{
+  w->bspline = 0;
+
+  if (w->A)
+    gsl_matrix_free(w->A);
+
   if (w->dB)
     gsl_matrix_free(w->dB);
 
   free(w);
-} /* gsl_bspline_free() */
+} /* gsl_bspline_deriv_free() */
 
 /*
 gsl_bspline_knots()
   Compute the knots from the given breakpoints:
 
-knots(1:k) = breakpts(1)
-knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
-knots(n+1:n+k) = breakpts(l + 1)
+   knots(1:k) = breakpts(1)
+   knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
+   knots(n+1:n+k) = breakpts(l + 1)
 
 where l is the number of polynomial pieces (l = nbreak - 1) and
-n = k + l - 1
+   n = k + l - 1
 (using matlab syntax for the arrays)
 
 The repeated knots at the beginning and end of the interval
@@ -244,8 +325,7 @@ gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w)
 
       for (i = 1; i < w->l; ++i)
         {
-          gsl_vector_set(w->knots, w->k - 1 + i,
-                         gsl_vector_get(breakpts, i));
+          gsl_vector_set(w->knots, w->k - 1 + i, gsl_vector_get(breakpts, i));
         }
 
       for (i = w->n; i < w->n + w->k; ++i)
@@ -323,8 +403,7 @@ Notes: The w->knots vector must be initialized prior to calling
 */
 
 int
-gsl_bspline_eval(const double x, gsl_vector *B,
-                 gsl_bspline_workspace *w)
+gsl_bspline_eval(const double x, gsl_vector *B, gsl_bspline_workspace *w)
 {
   if (B->size != w->n)
     {
@@ -358,7 +437,6 @@ gsl_bspline_eval(const double x, gsl_vector *B,
     }
 } /* gsl_bspline_eval() */
 
-
 /*
 gsl_bspline_eval_nonzero()
   Evaluate all non-zero B-spline functions at point x.
@@ -368,7 +446,7 @@ 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
+                 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.
@@ -381,18 +459,15 @@ Notes: 1) the w->knots vector must be initialized before calling
 
        2) On output, B contains
 
-              [B_{istart,k}, B_{istart+1,k},
-               ..., B_{iend-1,k}, B_{iend,k}]
+             [B_{istart,k}, B_{istart+1,k},
+             ..., B_{iend-1,k}, B_{iend,k}]
 
           evaluated at the given x.
 */
 
 int
-gsl_bspline_eval_nonzero(const double x,
-                         gsl_vector *Bk,
-                         size_t *istart,
-                         size_t *iend,
-                         gsl_bspline_workspace *w)
+gsl_bspline_eval_nonzero(const double x, gsl_vector *Bk, size_t *istart,
+                         size_t *iend, gsl_bspline_workspace *w)
 {
   if (Bk->size != w->k)
     {
@@ -400,9 +475,9 @@ gsl_bspline_eval_nonzero(const double x,
     }
   else
     {
-      size_t i;  /* spline index */
-      size_t j;  /* looping */
-      int flag  = 0; /* interval search flag */
+      size_t i;      /* spline index */
+      size_t j;      /* looping */
+      int flag = 0;  /* interval search flag */
       int error = 0; /* error flag */
 
       i = bspline_find_interval(x, &flag, w);
@@ -415,27 +490,26 @@ gsl_bspline_eval_nonzero(const double x,
       *istart = i - w->k + 1;
       *iend = i;
 
-      bspline_pppack_bsplvb(w->knots, w->k, 1, x, *iend,
-                            &j, w->deltal, w->deltar, Bk);
+      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()
+gsl_bspline_deriv_eval()
   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()
+This is a wrapper function for gsl_bspline_deriv_eval_nonzero()
 which formats the output in a nice way.
 
 Inputs: x      - point for evaluation
-        nderiv - number of derivatives to compute,
-                 inclusive.
+        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
+        w      - bspline derivative workspace
 
 Return: success or error
 
@@ -446,30 +520,29 @@ Notes: 1) The w->knots vector must be initialized prior to calling
 */
 
 int
-gsl_bspline_eval_deriv(const double x,
-                       const size_t nderiv,
-                       gsl_matrix *dB,
-                       gsl_bspline_workspace *w)
+gsl_bspline_deriv_eval(const double x, const size_t nderiv, gsl_matrix *dB,
+                       gsl_bspline_deriv_workspace *w)
 {
-  if (dB->size1 != w->n)
+  if (dB->size1 != w->bspline->n)
     {
       GSL_ERROR("dB matrix first dimension not of length n", GSL_EBADLEN);
     }
-  else if (dB->size2 < nderiv+1)
+  else if (dB->size2 < nderiv + 1)
     {
-      GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN);
+      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 */
+      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);
+      error = gsl_bspline_deriv_eval_nonzero(x, nderiv, w->dB, &istart, &iend,
+                                             w);
       if (error)
         {
           return error;
@@ -484,16 +557,16 @@ gsl_bspline_eval_deriv(const double x,
           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)
+          for (i = iend + 1; i < w->bspline->n; ++i)
             gsl_matrix_set(dB, i, j, 0.0);
         }
 
       return GSL_SUCCESS;
     }
-} /* gsl_bspline_eval_deriv() */
+} /* gsl_bspline_deriv_eval() */
 
 /*
-gsl_bspline_eval_deriv_nonzero()
+gsl_bspline_deriv_eval_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].
@@ -508,7 +581,7 @@ Inputs: x      - point at which to evaluate splines
         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
+        w      - bspline derivative workspace
 
 Return: success or error
 
@@ -517,14 +590,14 @@ Notes: 1) the w->knots vector must be initialized before calling
 
        2) On output, dB contains
 
-       [[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}]]
+            [[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}]]
 
-        evaluated at x.  B_{istart, k} is stored in dB(0,0).
-        Each additional column contains an additional derivative.
+          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.
@@ -533,48 +606,46 @@ Notes: 1) the w->knots vector must be initialized before calling
 */
 
 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)
+gsl_bspline_deriv_eval_nonzero(const double x, const size_t nderiv,
+                               gsl_matrix *dB, size_t *istart, size_t *iend,
+                               gsl_bspline_deriv_workspace *w)
 {
-  if (dB->size1 != w->k)
+  if (dB->size1 != w->bspline->k)
     {
       GSL_ERROR("dB matrix first dimension not of length k", GSL_EBADLEN);
     }
-  else if (dB->size2 < nderiv+1)
+  else if (dB->size2 < nderiv + 1)
     {
-      GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN);
+      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 flag = 0;  /* interval search flag */
       int error = 0; /* error flag */
       size_t min_nderivk;
 
-      i = bspline_find_interval(x, &flag, w);
-      error = bspline_process_interval_for_eval(x, &i, flag, w);
+      i = bspline_find_interval(x, &flag, w->bspline);
+      error = bspline_process_interval_for_eval(x, &i, flag, w->bspline);
       if (error)
         {
           return error;
         }
 
-      *istart = i - w->k + 1;
+      *istart = i - w->bspline->k + 1;
       *iend = i;
 
-      bspline_pppack_bsplvd(w->knots, w->k, x, *iend,
-                            w->deltal, w->deltar, w->A, dB, nderiv);
+      bspline_pppack_bsplvd(w->bspline->knots, w->bspline->k, x, *iend,
+                            w->bspline->deltal, w->bspline->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);
+       so we need to zero all requested higher order derivatives */
+      min_nderivk = GSL_MIN_INT(nderiv, w->bspline->k - 1);
       for (j = min_nderivk + 1; j <= nderiv; ++j)
         {
-          for (i = 0; i < w->k; ++i)
+          for (i = 0; i < w->bspline->k; ++i)
             {
               gsl_matrix_set(dB, i, j, 0.0);
             }
@@ -582,8 +653,7 @@ gsl_bspline_eval_deriv_nonzero(const double x,
 
       return GSL_SUCCESS;
     }
-} /* gsl_bspline_eval_deriv_nonzero() */
-
+} /* gsl_bspline_deriv_eval_nonzero() */
 
 /****************************************
  *          INTERNAL ROUTINES           *
@@ -591,10 +661,7 @@ gsl_bspline_eval_deriv_nonzero(const double x,
 
 /*
 bspline_find_interval()
-  Find knot interval such that
-
-  t_i <= x < t_{i + 1}
-
+  Find knot interval such that t_i <= x < t_{i + 1}
 where the t_i are knot values.
 
 Inputs: x    - x value
@@ -605,17 +672,16 @@ Return: i (index in w->knots corresponding to left limit of interval)
 
 Notes: The error conditions are reported as follows:
 
-Condition                        Return value        Flag
----------                        ------------        ----
-x < t_0                               0               -1
-t_i <= x < t_{i+1}                    i                0
-t_i < x = t_{i+1} = t_{n+k-1}         i                0
-t_{n+k-1} < x                       l+k-1             +1
+       Condition                        Return value        Flag
+       ---------                        ------------        ----
+       x < t_0                               0               -1
+       t_i <= x < t_{i+1}                    i                0
+       t_i < x = t_{i+1} = t_{n+k-1}         i                0
+       t_{n+k-1} < x                       l+k-1             +1
 */
 
 static inline size_t
-bspline_find_interval(const double x, int *flag,
-                      gsl_bspline_workspace *w)
+bspline_find_interval(const double x, int *flag, gsl_bspline_workspace *w)
 {
   size_t i;
 
@@ -639,8 +705,8 @@ bspline_find_interval(const double x, int *flag,
       if (ti <= x && x < tip1)
         break;
 
-      if (ti < x && x == tip1 &&
-          tip1 == gsl_vector_get(w->knots, w->k + w->l - 1))
+      if (ti < x && x == tip1 && tip1 == gsl_vector_get(w->knots, w->k + w->l
+          - 1))
         break;
     }
 
@@ -652,17 +718,16 @@ 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
+  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)
+*/
+static inline int
+bspline_process_interval_for_eval(const double x, size_t *i, const int flag,
+                                  gsl_bspline_workspace *w)
 {
   if (flag == -1)
     {
@@ -689,7 +754,6 @@ static inline int bspline_process_interval_for_eval(
   return GSL_SUCCESS;
 } /* bspline_process_interval_for_eval */
 
-
 /********************************************************************
  * PPPACK ROUTINES
  *
@@ -705,66 +769,71 @@ bspline_pppack_bsplvb()
 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)) .
+   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].
+   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.
+   (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().
+   (2) This routine contains no error checking.  Please use routines
+       like gsl_bspline_eval().
 */
 
 static void
@@ -778,7 +847,7 @@ bspline_pppack_bsplvb(const gsl_vector *t,
                       gsl_vector *deltar,
                       gsl_vector *biatx)
 {
-  size_t i;   /* looping */
+  size_t i; /* looping */
   double saved;
   double term;
 
@@ -797,12 +866,10 @@ bspline_pppack_bsplvb(const gsl_vector *t,
 
       for (i = 0; i <= *j; ++i)
         {
-          term = gsl_vector_get(biatx, i) /
-                 (gsl_vector_get(deltar, i) +
-                  gsl_vector_get(deltal, *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);
+          gsl_vector_set(biatx, i, saved + gsl_vector_get(deltar, i) * term);
 
           saved = gsl_vector_get(deltal, *j - i) * term;
         }
@@ -815,48 +882,47 @@ bspline_pppack_bsplvb(const gsl_vector *t,
 
 /*
 bspline_pppack_bsplvd()
-calculates value and derivs of all b-splines which do not vanish at x
+  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.)
+   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 .
+   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.
+   (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.
+   (2) This routine contains no error checking.
 */
 
 static void
@@ -877,36 +943,36 @@ bspline_pppack_bsplvd(const gsl_vector *t,
   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);
+  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.  */
+       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 < (int)k; ++j, ++jp1mid)
+          for (j = ideriv, jp1mid = 0; j < (int) k; ++j, ++jp1mid)
             {
-              gsl_matrix_set(dbiatx, j, ideriv,
-                             gsl_matrix_get(dbiatx, jp1mid, 0));
+              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);
+          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. */
+       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 < (int)k; ++i)
+      for (i = 0; i < (int) k; ++i)
         {
-          for (j = jlow; j < (int)k; ++j)
+          for (j = jlow; j < (int) k; ++j)
             {
               gsl_matrix_set(a, j, i, 0.0);
             }
@@ -915,7 +981,7 @@ bspline_pppack_bsplvd(const gsl_vector *t,
         }
 
       /* at this point, a(.,j) contains the b-coeffs for the j-th of the
-         k b-splines of interest here. */
+       k b-splines of interest here. */
       for (m = 1; m <= mhigh; ++m)
         {
           kmm = k - m;
@@ -924,40 +990,38 @@ bspline_pppack_bsplvd(const gsl_vector *t,
           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. */
+           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)
-                       );
+              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)
+               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)));
+                  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 < (int)k; ++i)
+           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 < (int) k; ++i)
             {
               sum = 0;
               jlow = GSL_MAX_INT(i, m);
-              for (j = jlow; j < (int)k; ++j)
+              for (j = jlow; j < (int) k; ++j)
                 {
-                  sum += gsl_matrix_get(a, j, i)
-                         * gsl_matrix_get(dbiatx, j, m);
+                  sum += gsl_matrix_get(a, j, i) * gsl_matrix_get(dbiatx, j, m);
                 }
               gsl_matrix_set(dbiatx, i, m, sum);
             }
diff --git a/bspline/gsl_bspline.h b/bspline/gsl_bspline.h
index 8a2ffc5..f7f8175 100644
--- a/bspline/gsl_bspline.h
+++ b/bspline/gsl_bspline.h
@@ -39,41 +39,53 @@
 __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 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) */
+    size_t n; /* number of bspline basis functions (l + k - 1) */
 
-    gsl_vector *knots;  /* knots vector */
+    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_matrix *A; /* work matrix */
+    gsl_vector *B; /* temporary spline results */
+    gsl_matrix *dB; /* temporary derivative results */
+} gsl_bspline_workspace;
+
+typedef struct
+{
+    gsl_bspline_workspace *bspline; /* associated bspline */
+    gsl_matrix *A; /* work matrix */
+    gsl_matrix *dB; /* temporary derivative results */
+} gsl_bspline_deriv_workspace;
 
 gsl_bspline_workspace *
 gsl_bspline_alloc(const size_t k, const size_t nbreak);
 
-void gsl_bspline_free(gsl_bspline_workspace *w);
+void
+gsl_bspline_free(gsl_bspline_workspace *w);
 
-size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * w);
-size_t gsl_bspline_order (gsl_bspline_workspace * w);
-size_t gsl_bspline_nbreak (gsl_bspline_workspace * w);
-double gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w);
+size_t
+gsl_bspline_ncoeffs(gsl_bspline_workspace * w);
+size_t
+gsl_bspline_order(gsl_bspline_workspace * w);
+size_t
+gsl_bspline_nbreak(gsl_bspline_workspace * w);
+double
+gsl_bspline_breakpoint(size_t i, gsl_bspline_workspace * w);
 
 int
 gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w);
 
-int gsl_bspline_knots_uniform(const double a, const double b,
-                              gsl_bspline_workspace *w);
+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_workspace *w);
+gsl_bspline_eval(const double x, gsl_vector *B, gsl_bspline_workspace *w);
 
 int
 gsl_bspline_eval_nonzero(const double x,
@@ -82,19 +94,25 @@ gsl_bspline_eval_nonzero(const double x,
                          size_t *iend,
                          gsl_bspline_workspace *w);
 
+gsl_bspline_deriv_workspace *
+gsl_bspline_deriv_alloc(gsl_bspline_workspace * bspline);
+
+void
+gsl_bspline_deriv_free(gsl_bspline_deriv_workspace *w);
+
 int
-gsl_bspline_eval_deriv(const double x,
+gsl_bspline_deriv_eval(const double x,
                        const size_t nderiv,
                        gsl_matrix *dB,
-                       gsl_bspline_workspace *w);
+                       gsl_bspline_deriv_workspace *w);
 
 int
-gsl_bspline_eval_deriv_nonzero(const double x,
+gsl_bspline_deriv_eval_nonzero(const double x,
                                const size_t nderiv,
                                gsl_matrix *dB,
                                size_t *istart,
                                size_t *iend,
-                               gsl_bspline_workspace *w);
+                               gsl_bspline_deriv_workspace *w);
 
 __END_DECLS
 
diff --git a/bspline/test.c b/bspline/test.c
index 393fbf3..5b511d3 100644
--- a/bspline/test.c
+++ b/bspline/test.c
@@ -26,7 +26,7 @@
 #include <gsl/gsl_nan.h>
 
 void
-test_bspline(gsl_bspline_workspace * bw)
+test_bspline(gsl_bspline_workspace * bw, gsl_bspline_deriv_workspace * dbw)
 {
   gsl_vector *B;
   gsl_matrix *dB;
@@ -68,7 +68,7 @@ test_bspline(gsl_bspline_workspace * bw)
     {
       double xi = a + (b - a) * (i / (n - 1.0));
       gsl_bspline_eval(xi, B, bw);
-      gsl_bspline_eval_deriv(xi, 0, dB, bw);
+      gsl_bspline_deriv_eval(xi, 0, dB, dbw);
 
       for (j = 0; j < ncoeffs; j++)
         {
@@ -100,8 +100,10 @@ main(int argc, char **argv)
         {
           double a = -1.23 * order, b = 45.6 * order;
           gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
+          gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw);
           gsl_bspline_knots_uniform(a, b, bw);
-          test_bspline(bw);
+          test_bspline(bw, dbw);
+          gsl_bspline_deriv_free(dbw);
           gsl_bspline_free(bw);
         }
     }
@@ -113,6 +115,7 @@ main(int argc, char **argv)
         {
           double a = -1.23 * order, b = 45.6 * order;
           gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints);
+          gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw);
           gsl_vector *k = gsl_vector_alloc(breakpoints);
           for (i = 0; i < breakpoints; i++)
             {
@@ -122,8 +125,9 @@ main(int argc, char **argv)
               gsl_vector_set(k, i, x);
             };
           gsl_bspline_knots(k, bw);
-          test_bspline(bw);
+          test_bspline(bw, dbw);
           gsl_vector_free(k);
+          gsl_bspline_deriv_free(dbw);
           gsl_bspline_free(bw);
         }
     }
@@ -143,6 +147,7 @@ main(int argc, char **argv)
     };
 
     gsl_bspline_workspace *bw = gsl_bspline_alloc(2, 3);
+    gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw);
     gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
                                       gsl_bspline_order(bw) + 1);
 
@@ -158,7 +163,7 @@ main(int argc, char **argv)
         /* 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);
+        gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, dbw);
         for (j = 0; j < gsl_bspline_ncoeffs(bw) ; ++j)
           {
             /* check basis function 1st deriv */
@@ -178,6 +183,7 @@ main(int argc, char **argv)
       }
 
     gsl_matrix_free(dB);
+    gsl_bspline_deriv_free(dbw);
     gsl_bspline_free(bw);
     gsl_vector_free(breakpts);
   }
@@ -213,6 +219,7 @@ main(int argc, char **argv)
     };
 
     gsl_bspline_workspace *bw = gsl_bspline_alloc(3, 5);
+    gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw);
 
     gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw),
                                       gsl_bspline_order(bw) + 1);
@@ -229,7 +236,7 @@ main(int argc, char **argv)
       {
         /* 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);
+        gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, dbw);
 
         /* check basis function evaluation */
         for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j)
@@ -255,6 +262,7 @@ main(int argc, char **argv)
       }
 
     gsl_matrix_free(dB);
+    gsl_bspline_deriv_free(dbw);
     gsl_bspline_free(bw);
     gsl_vector_free(breakpts);
   }
diff --git a/doc/bspline.texi b/doc/bspline.texi
index 36b3dce..b2cdb1b 100644
--- a/doc/bspline.texi
+++ b/doc/bspline.texi
@@ -78,6 +78,8 @@ be readily obtained from a least-squares fit.
 @section Initializing the B-splines solver
 @cindex basis splines, initializing
 
+Using B-splines requires a gsl_bspline_workspace:
+
 @deftypefun {gsl_bspline_workspace *} gsl_bspline_alloc (const size_t @var{k}, const size_t @var{nbreak})
 This function allocates a workspace for computing B-splines of order
 @var{k}. The number of breakpoints is given by @var{nbreak}. This
@@ -88,6 +90,21 @@ are specified by @math{k = 4}. The size of the workspace is
 
 @deftypefun void gsl_bspline_free (gsl_bspline_workspace * @var{w})
 This function frees the memory associated with the workspace @var{w}.
+Any associated gsl_bspline_deriv_workspace should be freed beforehand.
+@end deftypefun
+
+Evaluation of B-spline basis function derivatives additionally requires
+a gsl_bspline_deriv_workspace:
+
+@deftypefun {gsl_bspline_deriv_workspace *} gsl_bspline_deriv_alloc (gsl_bspline_workspace * bspline)
+This function allocates a workspace for computing the B-spline basis
+function derivatives for @var{bspline}.  The size of the workspace is
+@math{O(2k^2)}.
+@end deftypefun
+
+@deftypefun void gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * @var{w})
+This function frees the memory associated with the workspace @var{w}.
+The associated gsl_bspline_workspace should be freed afterwards.
 @end deftypefun
 
 @node Constructing the knots vector
@@ -140,7 +157,7 @@ This function returns the number of B-spline coefficients given by
 @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})
+@deftypefun int gsl_bspline_deriv_eval (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}
@@ -153,7 +170,7 @@ 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})
+@deftypefun int gsl_bspline_deriv_eval_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

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

At Fri, 5 Dec 2008 12:10:58 -0700,
Patrick Alken wrote:
> Ah I hadn't considered this issue. If there are other areas of GSL
> which have the same issues maybe we should keep a list so when there
> is a major new release like 2.0 we could integrate these "extra"
> workspaces?

I've put the following in the toplevel TODO file. We should
probably make 2.0 the next release after this.

Changes for Release 2.0

Break binary compatibility, but keep source compatibility.

* Add a 'void *' to all workspaces, to allow for future changes.

* Disable deprecated functions

* Fix up the workspace_alloc functions so they have consistent names
  (add functions where needed, don't remove)


> Of course one problem with that is if we have a separate workspace
> for bspline derivatives it may never be convenient to integrate it
> with the normal bspline workspace.

True, unfortunately I don't see any easy way round that.

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

I've added this with some minor changes (I've moved the bspline
workspace out of the derivative workspace, to avoid any confusion
about the ownership of the memory -- so the derivative functions take
the normal bspline workspace and a derivative workspace as arguments).

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Rhys Ulerich]

[-- Attachment #1: Type: text/plain, Size: 598 bytes --]

Makes sense, and the edits/documentation cleanups look consistent throughout.

There a couple of lingering comments that'll be misleading to anyone
reading the bspline.c file.  One final patch, attached, will kill
those.

Thanks,
Rhys


On Mon, Dec 8, 2008 at 3:11 PM, Brian Gough <bjg@network-theory.co.uk> wrote:
> I've added this with some minor changes (I've moved the bspline
> workspace out of the derivative workspace, to avoid any confusion
> about the ownership of the memory -- so the derivative functions take
> the normal bspline workspace and a derivative workspace as arguments).
>
>

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: remove-associated-workspace-co.patch --]
[-- Type: text/x-diff; name=remove-associated-workspace-co.patch, Size: 1118 bytes --]

remove associated workspace comments

From: Rhys Ulerich <rhys.ulerich@gmail.com>


---

 bspline/bspline.c |    6 ++----
 1 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/bspline/bspline.c b/bspline/bspline.c
index 465943b..ab7c742 100644
--- a/bspline/bspline.c
+++ b/bspline/bspline.c
@@ -127,7 +127,7 @@ gsl_bspline_deriv_alloc()
   Allocate space for a bspline derivative workspace. The size of the
 workspace is O(2k^2)
 
-Inputs: bspline - gsl_bspline_workspace of the associated bspline
+Inputs: k      - spline order (cubic = 4)
 
 Return: pointer to workspace
 */
@@ -208,7 +208,6 @@ gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
 /*
 gsl_bspline_free()
   Free a gsl_bspline_workspace.
-Any associated gsl_bspline_deriv_workspace should be freed beforehand.
 
 Inputs: w - workspace to free
 
@@ -228,9 +227,8 @@ gsl_bspline_free (gsl_bspline_workspace * w)
 /*
 gsl_bspline_deriv_free()
   Free a gsl_bspline_deriv_workspace.
-The associated gsl_bspline_workspace should be freed afterwards.
 
-Inputs: w - workspace to free
+Inputs: dw - workspace to free
 
 Return: none
 */

Re: [Help-gsl] Submitted patch to add bspline derivative capabilities

[Brian Gough]

At Mon, 8 Dec 2008 15:43:05 -0600,
Rhys Ulerich wrote:
> 
> Makes sense, and the edits/documentation cleanups look consistent throughout.
> 
> There a couple of lingering comments that'll be misleading to anyone
> reading the bspline.c file.  One final patch, attached, will kill
> those.
> 

Thanks, I've added that.