/* 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;

}				/* bspline_process_interval_for_eval */

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

 * PPPACK ROUTINES

 *

 * The routines herein deliberately avoid using the bspline workspace,

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

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

 * parallel scenarios.

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

/*

bspline_pppack_bsplvb()

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

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

Parameters:

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

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

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

   jhigh  -

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

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

            are to be returned.  index  is used to avoid

            recalculations when several columns of the triangular

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

            or in  bsplvd ).  precisely,

            if  index = 1 ,

               the calculation starts from scratch and the entire

               triangular array of b-spline values of orders

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

               column by column .

            if  index = 2 ,

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

               are generated, the assumption being that biatx, j,

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

               at the previous call.

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

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

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

   left   - an integer chosen (usually) so that

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

   j      - (output) a working scalar for indexing

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

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

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

            containing the value at  x  of the polynomial of order

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

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

Method:

   the recurrence relation

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

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

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

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

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

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

   facts that

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

   and that

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

   are used. the particular organization of the calculations follows

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

Notes:

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

       j, deltal, deltar rewritten as input parameters and

       utilizing zero-based indexing.

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

       like gsl_bspline_eval().

*/

static void

bspline_pppack_bsplvb (const gsl_vector * t,

		       const size_t jhigh,

		       const size_t index,

		       const double x,

		       const size_t left,

		       size_t * j,

		       gsl_vector * deltal,

		       gsl_vector * deltar, gsl_vector * biatx)

{

  size_t i;			/* looping */

  double saved;

  double term;

  if (index == 1)

    {

      *j = 0;

      gsl_vector_set (biatx, 0, 1.0);

    }

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

    {

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

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

      saved = 0.0;

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

	{

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

					      + gsl_vector_get (deltal,

								*j - i));

	  gsl_vector_set (biatx, i,

			  saved + gsl_vector_get (deltar, i) * term);

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

	}

      gsl_vector_set (biatx, *j + 1, saved);

    }

  return;

}				/* gsl_bspline_pppack_bsplvb */

/*

bspline_pppack_bsplvd()

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

Parameters:

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

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

   x      - the point at which these values are sought

   left   - an integer indicating the left endpoint of the interval

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

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

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

            division by zero will result otherwise (in  bsplvb).

            also, the output is as advertised only if

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

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

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

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

            derivatives of a certain order of the k b-splines

            of interest.

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

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

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

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

            their derivatives up to AND INCLUDING the nderiv-th

            are asked for. (nderiv is replaced internally by the

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

Method:

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

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

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

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

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

   dbiatx  to produce the desired values .

Notes:

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

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

       to bspline_pppack_bsplvb) and utilizing zero-based indexing.

   (2) This routine contains no error checking.

*/

static void

bspline_pppack_bsplvd (const gsl_vector * t,

		       const size_t k,

		       const double x,

		       const size_t left,

		       gsl_vector * deltal,

		       gsl_vector * deltar,

		       gsl_matrix * a,

		       gsl_matrix * dbiatx, const size_t nderiv)

{

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

  double factor, fkmm, sum;

  size_t bsplvb_j;

  gsl_vector_view dbcol = gsl_matrix_column (dbiatx, 0);

  mhigh = GSL_MIN_INT (nderiv, k - 1);

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

			 &dbcol.vector);

  if (mhigh > 0)

    {

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

         values for the current order. these are stored in column

         k-current order before bsplvb is called to put values

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

      ideriv = mhigh;

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

	{

	  for (j = ideriv, jp1mid = 0; j < (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++)

	    {