/* bspline/bspline.c

*

* Copyright (C) 2006, 2007, 2008, 2009 Patrick Alken

* Copyright (C) 2008 Rhys Ulerich

*

* This program is free software; you can redistribute it and/or modify

* it under the terms of the GNU General Public License as published by

* the Free Software Foundation; either version 3 of the License, or (at

* your option) any later version.

*

* This program is distributed in the hope that it will be useful, but

* WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

* General Public License for more details.

*

* You should have received a copy of the GNU General Public License

* along with this program; if not, write to the Free Software

* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

*/

#include <config.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_bspline.h>

/*

* This module contains routines related to calculating B-splines.

* The algorithms used are described in

*

* [1] Carl de Boor, "A Practical Guide to Splines", Springer

* Verlag, 1978.

*

* The bspline_pppack_* internal routines contain code adapted from

*

* [2] "PPPACK - Piecewise Polynomial Package",

* http://www.netlib.org/pppack/

*

*/

#include "bspline.h"

/*

gsl_bspline_alloc()

Allocate space for a bspline workspace. The size of the

workspace is O(5k nbreak)

Inputs: k - spline order (cubic = 4)

nbreak - number of breakpoints

Return: pointer to workspace

*/

gsl_bspline_workspace *

gsl_bspline_alloc (const size_t k, const size_t nbreak)

{

if (k == 0)

{

GSL_ERROR_NULL ("k must be at least 1", GSL_EINVAL);

}

else if (nbreak < 2)

{

GSL_ERROR_NULL ("nbreak must be at least 2", GSL_EINVAL);

}

else

{

gsl_bspline_workspace *w;

w = calloc (1, sizeof (gsl_bspline_workspace));

if (w == 0)

{

GSL_ERROR_NULL ("failed to allocate space for workspace",

GSL_ENOMEM);

}

w->k = k;

w->km1 = k - 1;

w->nbreak = nbreak;

w->l = nbreak - 1;

w->n = w->l k - 1;

w->knots = gsl_vector_alloc (w->n k);

if (w->knots == 0)

{

gsl_bspline_free (w);

GSL_ERROR_NULL ("failed to allocate space for knots vector",

GSL_ENOMEM);

}

w->deltal = gsl_vector_alloc (k);

if (w->deltal == 0)

{

gsl_bspline_free (w);

GSL_ERROR_NULL ("failed to allocate space for deltal vector",

GSL_ENOMEM);

}

w->deltar = gsl_vector_alloc (k);

if (w->deltar == 0)

{

gsl_bspline_free (w);

GSL_ERROR_NULL ("failed to allocate space for deltar vector",

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);

}

w->A = gsl_matrix_alloc (k, k);

if (w->A == 0)

{

gsl_bspline_free (w);

GSL_ERROR_NULL

("failed to allocate space for derivative work matrix",

GSL_ENOMEM);

}

w->dB = gsl_matrix_alloc (k, k 1);

if (w->dB == 0)

{

gsl_bspline_free (w);

GSL_ERROR_NULL

("failed to allocate space for temporary derivative matrix",

GSL_ENOMEM);

}

return w;

}

} /* gsl_bspline_alloc() */

/*

gsl_bspline_free()

Free a gsl_bspline_workspace.

Inputs: w - workspace to free

Return: none

*/

void

gsl_bspline_free (gsl_bspline_workspace * w)

{

RETURN_IF_NULL (w);

if (w->knots)

gsl_vector_free (w->knots);

if (w->deltal)

gsl_vector_free (w->deltal);

if (w->deltar)

gsl_vector_free (w->deltar);

if (w->B)

gsl_vector_free (w->B);

if (w->A)

gsl_matrix_free(w->A);

if (w->dB)

gsl_matrix_free(w->dB);

free (w);

} /* gsl_bspline_free() */

/* Return number of coefficients */

size_t

gsl_bspline_ncoeffs (gsl_bspline_workspace * w)

{

return w->n;

}

/* Return order */

size_t

gsl_bspline_order (gsl_bspline_workspace * w)

{

return w->k;

}

/* Return number of breakpoints */

size_t

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)

{

size_t j = i w->k - 1;

return gsl_vector_get (w->knots, j);

}

/*

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)

where l is the number of polynomial pieces (l = nbreak - 1) and

n = k l - 1

(using matlab syntax for the arrays)

The repeated knots at the beginning and end of the interval

correspond to the continuity condition there. See pg. 119

of [1].

Inputs: breakpts - breakpoints

w - bspline workspace

Return: success or error

*/

int

gsl_bspline_knots (const gsl_vector * breakpts, gsl_bspline_workspace * w)

{

if (breakpts->size != w->nbreak)

{

GSL_ERROR ("breakpts vector has wrong size", GSL_EBADLEN);

}

else

{

size_t i; /* looping */

for (i = 0; i < w->k; i )

gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, 0));

for (i = 1; i < w->l; 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 )

gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, w->l));

return GSL_SUCCESS;

}

} /* gsl_bspline_knots() */

/*

gsl_bspline_knots_uniform()

Construct uniformly spaced knots on the interval [a,b] using

the previously specified number of breakpoints. 'a' is the position

of the first breakpoint and 'b' is the position of the last

breakpoint.

Inputs: a - left side of interval

b - right side of interval

w - bspline workspace

Return: success or error

Notes: 1) w->knots is modified to contain the uniformly spaced

knots

2) The knots vector is set up as follows (using octave syntax):

knots(1:k) = a

knots(k 1:k l-1) = a i*delta, i = 1 .. l - 1

knots(n 1:n k) = b

*/

int

gsl_bspline_knots_uniform (const double a, const double b,

gsl_bspline_workspace * w)

{

size_t i; /* looping */

double delta; /* interval spacing */

double x;

delta = (b - a) / (double) w->l;

for (i = 0; i < w->k; i )

gsl_vector_set (w->knots, i, a);

x = a delta;

for (i = 0; i < w->l - 1; i )

{

gsl_vector_set (w->knots, w->k i, x);

x = delta;

}

for (i = w->n; i < w->n w->k; i )

gsl_vector_set (w->knots, i, b);

return GSL_SUCCESS;

} /* gsl_bspline_knots_uniform() */

/*

gsl_bspline_eval()

Evaluate the basis functions B_i(x) for all i. This is

a wrapper function for gsl_bspline_eval_nonzero() which

formats the output in a nice way.

Inputs: x - point for evaluation

B - (output) where to store B_i(x) values

the length of this vector is

n = nbreak k - 2 = l k - 1 = w->n

w - bspline workspace

Return: success or error

Notes: The w->knots vector must be initialized prior to calling

this function (see gsl_bspline_knots())

*/

int

gsl_bspline_eval (const double x, gsl_vector * B, gsl_bspline_workspace * w)

{

if (B->size != w->n)

{

GSL_ERROR ("vector B not of length n", GSL_EBADLEN);

}

else

{

size_t i; /* looping */

size_t istart; /* first non-zero spline for x */

size_t iend; /* last non-zero spline for x, knot for x */

int error; /* error handling */

/* find all non-zero B_i(x) values */

error = gsl_bspline_eval_nonzero (x, w->B, &istart, &iend, w);

if (error)

return error;

/* store values in appropriate part of given vector */

for (i = 0; i < istart; i )

gsl_vector_set (B, i, 0.0);

for (i = istart; i <= iend; i )

gsl_vector_set (B, i, gsl_vector_get (w->B, i - istart));

for (i = iend 1; i < w->n; i )

gsl_vector_set (B, i, 0.0);

return GSL_SUCCESS;

}

} /* gsl_bspline_eval() */

/*

gsl_bspline_eval_nonzero()

Evaluate all non-zero B-spline functions at point x.

These are the B_i(x) for i in [istart, iend].

Always B_i(x) = 0 for i < istart and for i > iend.

Inputs: x - point at which to evaluate splines

Bk - (output) where to store B-spline values (length k)

istart - (output) B-spline function index of

first non-zero basis for given x

iend - (output) B-spline function index of

last non-zero basis for given x.

This is also the knot index corresponding to x.

w - bspline workspace

Return: success or error

Notes: 1) the w->knots vector must be initialized before calling

this function

2) On output, B contains

[B_{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)

{

if (Bk->size != w->k)

{

GSL_ERROR ("Bk vector length does not match order k", GSL_EBADLEN);

}

else

{

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);

error = bspline_process_interval_for_eval (x, &i, flag, w);

if (error)

return error;

*istart = i - w->k 1;

*iend = i;

bspline_pppack_bsplvb (w->knots, w->k, 1, x, *iend, &j, w->deltal,

w->deltar, Bk);

return GSL_SUCCESS;

}

} /* gsl_bspline_eval_nonzero() */

/*

gsl_bspline_deriv_eval()

Evaluate d^j/dx^j B_i(x) for all i, 0 <= j <= nderiv.

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.

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 derivative workspace

Return: success or error

Notes: 1) The w->knots vector must be initialized prior to calling

this function (see gsl_bspline_knots())

2) based on PPPACK's bsplvd

*/

int

gsl_bspline_deriv_eval (const double x, const size_t nderiv,

gsl_matrix * dB, gsl_bspline_workspace * w)

{

if (dB->size1 != w->n)

{

GSL_ERROR ("dB matrix first dimension not of length n", GSL_EBADLEN);

}

else if (dB->size2 < nderiv 1)

{

GSL_ERROR

("dB matrix second dimension must be at least length nderiv 1",

GSL_EBADLEN);

}

else

{

size_t i; /* looping */

size_t j; /* looping */

size_t istart; /* first non-zero spline for x */

size_t iend; /* last non-zero spline for x, knot for x */

int error; /* error handling */

/* find all non-zero d^j/dx^j B_i(x) values */

error =

gsl_bspline_deriv_eval_nonzero (x, nderiv, w->dB, &istart, &iend, w);

if (error)

return error;

/* store values in appropriate part of given matrix */

for (j = 0; j <= nderiv; j )

{

for (i = 0; i < istart; i )

gsl_matrix_set (dB, i, j, 0.0);

for (i = istart; i <= iend; i )

gsl_matrix_set (dB, i, j, gsl_matrix_get (w->dB, i - istart, j));

for (i = iend 1; i < w->n; i )

gsl_matrix_set (dB, i, j, 0.0);

}

return GSL_SUCCESS;

}

} /* gsl_bspline_deriv_eval() */

/*

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].

Always d^j/dx^j B_i(x) = 0 for i < istart and for i > iend.

Inputs: x - point at which to evaluate splines

nderiv - number of derivatives to request, inclusive

dB - (output) where to store dB-spline derivatives

(size k by nderiv 1)

istart - (output) B-spline function index of

first non-zero basis for given x

iend - (output) B-spline function index of

last non-zero basis for given x.

This is also the knot index corresponding to x.

w - bspline derivative workspace

Return: success or error

Notes: 1) the w->knots vector must be initialized before calling

this function

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}]]

evaluated at x. B_{istart, k} is stored in dB(0,0).

Each additional column contains an additional derivative.

3) Note that the zero-th column of the result contains the

0th derivative, which is simply a function evaluation.

4) based on PPPACK's bsplvd

*/

int

gsl_bspline_deriv_eval_nonzero (const double x, const size_t nderiv,

gsl_matrix * dB, size_t * istart,

size_t * iend, gsl_bspline_workspace * w)

{

if (dB->size1 != w->k)

{

GSL_ERROR ("dB matrix first dimension not of length k", GSL_EBADLEN);

}

else if (dB->size2 < nderiv 1)

{

GSL_ERROR

("dB matrix second dimension must be at least length nderiv 1",

GSL_EBADLEN);

}

else

{

size_t i; /* spline index */

size_t j; /* looping */

int flag = 0; /* interval search flag */

int error = 0; /* error flag */

size_t min_nderivk;

i = bspline_find_interval (x, &flag, w);

error = bspline_process_interval_for_eval (x, &i, flag, w);

if (error)

return error;

*istart = i - w->k 1;

*iend = i;

bspline_pppack_bsplvd (w->knots, w->k, x, *iend,

w->deltal, w->deltar, w->A, dB, nderiv);

/* An order k b-spline has at most k-1 nonzero derivatives

so we need to zero all requested higher order derivatives */

min_nderivk = GSL_MIN_INT (nderiv, w->k - 1);

for (j = min_nderivk 1; j <= nderiv; j )

{

for (i = 0; i < w->k; i )

gsl_matrix_set (dB, i, j, 0.0);

}

return GSL_SUCCESS;

}

} /* gsl_bspline_deriv_eval_nonzero() */

/****************************************

* INTERNAL ROUTINES *

****************************************/

/*

bspline_find_interval()

Find knot interval such that t_i <= x < t_{i 1}

where the t_i are knot values.

Inputs: x - x value

flag - (output) error flag

w - bspline workspace

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

*/

static inline size_t

bspline_find_interval (const double x, int *flag, gsl_bspline_workspace * w)

{

size_t i;

if (x < gsl_vector_get (w->knots, 0))

{

*flag = -1;

return 0;

}

/* find i such that t_i <= x < t_{i 1} */

for (i = w->k - 1; i < w->k w->l - 1; i )

{

const double ti = gsl_vector_get (w->knots, i);

const double tip1 = gsl_vector_get (w->knots, i 1);

if (tip1 < ti)

{

GSL_ERROR ("knots vector is not increasing", GSL_EINVAL);

}

if (ti <= x && x < tip1)

break;

if (ti < x && x == tip1 && tip1 == gsl_vector_get (w->knots, w->k w->l

- 1))

break;

}

if (i == w->k w->l - 1)

*flag = 1;

else

*flag = 0;

return i;

} /* bspline_find_interval() */

/*

bspline_process_interval_for_eval()

Consumes an x location, left knot from bspline_find_interval, flag

from bspline_find_interval, and a workspace. Checks that x lies within

the splines' knots, enforces some endpoint continuity requirements, and

avoids divide by zero errors in the underlying bspline_pppack_* functions.

*/

static inline int

bspline_process_interval_for_eval (const double x, size_t * i, const int flag,

gsl_bspline_workspace * w)

{

if (flag == -1)

{

GSL_ERROR ("x outside of knot interval", GSL_EINVAL);

}

else if (flag == 1)

{

if (x <= gsl_vector_get (w->knots, *i) GSL_DBL_EPSILON)

{

*i -= 1;

}

else

{

GSL_ERROR ("x outside of knot interval", GSL_EINVAL);

}

}

if (gsl_vector_get (w->knots, *i) == gsl_vector_get (w->knots, *i 1))

{

GSL_ERROR ("knot(i) = knot(i 1) will result in division by zero", GSL_EINVAL);

}

return GSL_SUCCESS;

}

/********************************************************************

* PPPACK ROUTINES

*

* The routines herein deliberately avoid using the bspline workspace,

* choosing instead to pass all work areas explicitly. This allows

* others to more easily adapt these routines to low memory or

* parallel scenarios.

********************************************************************/

/*

bspline_pppack_bsplvb()

calculates the value of all possibly nonzero b-splines at x of order

jout = max( jhigh , (j 1)*(index-1) ) with knot sequence t.

Parameters:

t - knot sequence, of length left jout , assumed to be

nondecreasing. assumption t(left).lt.t(left 1).

division by zero will result if t(left) = t(left 1)

jhigh -

index - integers which determine the order jout = max(jhigh,

(j 1)*(index-1)) of the b-splines whose values at x

are to be returned. index is used to avoid

recalculations when several columns of the triangular

array of b-spline values are needed (e.g., in bsplpp

or in bsplvd ). precisely,

if index = 1 ,

the calculation starts from scratch and the entire

triangular array of b-spline values of orders

1,2,...,jhigh is generated order by order , i.e.,

column by column .

if index = 2 ,

only the b-spline values of order j 1, j 2, ..., jout

are generated, the assumption being that biatx, j,

deltal, deltar are, on entry, as they were on exit

at the previous call.

in particular, if jhigh = 0, then jout = j 1, i.e.,

just the next column of b-spline values is generated.

x - the point at which the b-splines are to be evaluated.

left - an integer chosen (usually) so that

t(left) .le. x .le. t(left 1).

j - (output) a working scalar for indexing

deltal - (output) a working area which must be of length at least jout

deltar - (output) a working area which must be of length at least jout

biatx - (output) array of length jout, with biatx(i)

containing the value at x of the polynomial of order

jout which agrees with the b-spline b(left-jout i,jout,t)

on the interval (t(left), t(left 1)) .

Method:

the recurrence relation

x - t(i) t(i j 1) - x

b(i,j 1)(x) = -----------b(i,j)(x) ---------------b(i 1,j)(x)

t(i j)-t(i) t(i j 1)-t(i 1)

is used (repeatedly) to generate the (j 1)-vector b(left-j,j 1)(x),

...,b(left,j 1)(x) from the j-vector b(left-j 1,j)(x),...,

b(left,j)(x), storing the new values in biatx over the old. the

facts that

b(i,1) = 1 if t(i) .le. x .lt. t(i 1)

and that

b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i j)

are used. the particular organization of the calculations follows

algorithm (8) in chapter x of [1].

Notes:

(1) This is a direct translation of PPPACK's bsplvb routine with

j, deltal, deltar rewritten as input parameters and

utilizing zero-based indexing.

(2) This routine contains no error checking. Please use routines

like gsl_bspline_eval().

*/

static void

bspline_pppack_bsplvb (const gsl_vector * t,

gsl_vector * deltar, gsl_vector * biatx)

{

size_t i; /* looping */

double saved;

double term;

if (index == 1)

{

*j = 0;

gsl_vector_set (biatx, 0, 1.0);

}

for ( /* NOP */ ; *j < jhigh - 1; *j = 1)

{

gsl_vector_set (deltar, *j, gsl_vector_get (t, left *j 1) - x);

gsl_vector_set (deltal, *j, x - gsl_vector_get (t, left - *j));

saved = 0.0;

for (i = 0; i <= *j; i )

{

term = gsl_vector_get (biatx, i) / (gsl_vector_get (deltar, i)

gsl_vector_get (deltal, *j - i));

gsl_vector_set (biatx, i, saved gsl_vector_get (deltar, i) * term);

saved = gsl_vector_get (deltal, *j - i) * term;

}

gsl_vector_set (biatx, *j 1, saved);

}

return;

}

/*

bspline_pppack_bsplvd()

calculates value and derivs of all b-splines which do not vanish at x

Parameters:

t - the knot array, of length left k (at least)

k - the order of the b-splines to be evaluated

x - the point at which these values are sought

left - an integer indicating the left endpoint of the interval

of interest. the k b-splines whose support contains the

interval (t(left), t(left 1)) are to be considered.

it is assumed that t(left) .lt. t(left 1)

division by zero will result otherwise (in bsplvb).

also, the output is as advertised only if

t(left) .le. x .le. t(left 1) .

deltal - a working area which must be of length at least k

deltar - a working area which must be of length at least k

a - an array of order (k,k), to contain b-coeffs of the

derivatives of a certain order of the k b-splines

of interest.

dbiatx - an array of order (k,nderiv). its entry (i,m) contains

value of (m)th derivative of (left-k i)-th b-spline

of order k for knot sequence t, i=1,...,k, m=0,...,nderiv.

nderiv - an integer indicating that values of b-splines and

their derivatives up to AND INCLUDING the nderiv-th

are asked for. (nderiv is replaced internally by the

integer mhigh in (1,k) closest to it.)

Method:

values at x of all the relevant b-splines of order k,k-1,..., k 1-nderiv

are generated via bsplvb and stored temporarily in dbiatx. then, the

b-coeffs of the required derivatives of the b-splines of interest are

generated by differencing, each from the preceeding one of lower order,

and combined with the values of b-splines of corresponding order in

dbiatx to produce the desired values .

Notes:

(1) This is a direct translation of PPPACK's bsplvd routine with

deltal, deltar rewritten as input parameters (to later feed them

to bspline_pppack_bsplvb) and utilizing zero-based indexing.

(2) This routine contains no error checking.

*/

static void

bspline_pppack_bsplvd (const gsl_vector * t,

gsl_matrix * dbiatx, const size_t nderiv)

{

int i, ideriv, il, j, jlow, jp1mid, kmm, ldummy, m, mhigh;

double factor, fkmm, sum;

size_t bsplvb_j;

gsl_vector_view dbcol = gsl_matrix_column (dbiatx, 0);

mhigh = GSL_MIN_INT (nderiv, k - 1);

bspline_pppack_bsplvb (t, k - mhigh, 1, x, left, &bsplvb_j, deltal, deltar,

&dbcol.vector);

if (mhigh > 0)

{

/* the first column of dbiatx always contains the b-spline

values for the current order. these are stored in column

k-current order before bsplvb is called to put values

for the next higher order on top of it. */

ideriv = mhigh;

for (m = 1; m <= mhigh; m )

{

for (j = ideriv, jp1mid = 0; j < (int) k; j , jp1mid )

{

gsl_matrix_set (dbiatx, j, ideriv, gsl_matrix_get (dbiatx, jp1mid, 0));

}

ideriv--;

bspline_pppack_bsplvb (t, k - ideriv, 2, x, left, &bsplvb_j, deltal,

deltar, &dbcol.vector);

}

/* at this point, b(left-k i, k 1-j)(x) is in dbiatx(i,j)

for i=j,...,k-1 and j=0,...,mhigh. in particular, the

first column of dbiatx is already in final form. to obtain

corresponding derivatives of b-splines in subsequent columns,

generate their b-repr. by differencing, then evaluate at x. */

jlow = 0;

for (i = 0; i < (int) k; i )

{

for (j = jlow; j < (int) k; j )

{

gsl_matrix_set (a, j, i, 0.0);

}

jlow = i;

gsl_matrix_set (a, i, i, 1.0);

}

/* at this point, a(.,j) contains the b-coeffs for the j-th of the

k b-splines of interest here. */

for (m = 1; m <= mhigh; m )

{

kmm = k - m;

fkmm = (float) kmm;

il = left;

i = k - 1;

/* for j=1,...,k, construct b-coeffs of (m)th derivative

of b-splines from those for preceding derivative by

differencing and store again in a(.,j) . the fact that

a(i,j) = 0 for i .lt. j is used. */

for (ldummy = 0; ldummy < kmm; ldummy )

{

factor = fkmm / (gsl_vector_get (t, il kmm) -

gsl_vector_get (t, il));

/* the assumption that t(left).lt.t(left 1) makes

denominator in factor nonzero. */

for (j = 0; j <= i; j )

{

gsl_matrix_set (a, i, j,

factor * (gsl_matrix_get (a, i, j)

- gsl_matrix_get (a, i - 1, j)));

}

il--;

i--;

}

/* for i=1,...,k, combine b-coeffs a(.,i) with b-spline values

stored in dbiatx(.,m) to get value of (m)th derivative

of i-th b-spline (of interest here) at x, and store in

dbiatx(i,m). storage of this value over the value of a

b-spline of order m there is safe since the remaining

b-spline derivatives of the same order do not use this

value due to the fact that a(j,i) = 0 for j .lt. i . */

for (i = 0; i < (int) k; i )

{

sum = 0;

jlow = GSL_MAX_INT (i, m);

for (j = jlow; j < (int) k; j )

{

sum = gsl_matrix_get (a, j, i) * gsl_matrix_get (dbiatx, j, m);

}

gsl_matrix_set (dbiatx, i, m, sum);

}

}

}

return;

}