* Re: [Help-gsl] Submitted patch to add bspline derivative capabilities
[not found] ` <4a00655d0812051440v46865adbsfcc10bb7dce28d98@mail.gmail.com>
@ 2008-12-11 15:18 ` Rhys Ulerich
2008-12-12 11:14 ` Brian Gough
0 siblings, 1 reply; 13+ messages in thread
From: Rhys Ulerich @ 2008-12-11 15:18 UTC (permalink / raw)
To: Brian Gough, Patrick Alken; +Cc: gsl-discuss
[-- 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
^ permalink raw reply [flat|nested] 13+ messages in thread