Multiple uniform splines

Spline creation

When evaluating many splines at the same point (such as in filling in a row of a determinant matrix in quantum Monte Carlo calculations), a signficant speedup can be achieved by using the multiple spline version of the uniform B-spline routines. The creation routines are similar to those for the single splines, which the total number of splines is passed and the raw data is NOT passed in at time of creation. Rather, we can the create_multi_UBspline_* routine to allocate the memory for the splines. We then call the set_multi_UBspline_* routine for each spline, passing the the spline index and the raw data to be interpolated.

Base data types

Each spline creation and evaluation routines have four versions, corresponding to the four support base data types. Following the LAPACK naming convention, these types are specified by single-letter codes:

s: single-precision real (float)
d: double-precision real (double)
c: single-precision complex (complex_float)
z: double-precision complex (complex_double)

Grid dimensions

The uniform grid structure, has the following elements, which must be specified:

double start : the first grid point
double end : the last grid point
int num : the number of grid points

Boundary conditions

The boundary conditions at the first and last grid point must be specified. They are specifed with BCtype_x structures, where x is one of {s, d, c, z}, as described above. For the real types (s and d) the structure contains the following elements

BCtype_s and BCtype_d data structure elements
Type	Name	Description
bc_code	lCode	"left" boundary condition code
bc_code	lCode	"right" boundary condition code
d_type	lVal	"left" boundary condition value
d_type	rVal	"right" boundary condition value

For the complex types, we must specify both the real and imaginary parts:

BCtype_s and BCtype_d data structure elements
Type	Name	Description
bc_code	lCode	"left" boundary condition code
bc_code	lCode	"right" boundary condition code
d_type	lVal_r	"left" boundary condition real part
d_type	lVal_i	"left" boundary condition imag part
d_type	rVal_r	"right" boundary condition real part
d_type	rVal_i	"right" boundary condition imag part

lCode and lVal specify the boundary conditions at the first grid point (the "left" boundary), while rCode and rVal specify the boundary conditions on the last grid point (the "right" boundary).
bc_code is one of the enumerated value, {PERIODIC, DERIV1, DERIV2, FLAT, NATURAL, ANTIPERIODIC }.
d_type is the the C type corresponding to {s, d, c, z}, i.e. {float, double, complex_float, complex_double}

The codes have the following meaning

Code	Meaning
PERIODIC	Use periodic boundary conditions. The value, first derivative and second derivative at the left boundary match those at the right boundary.
DERIV1	The value of the first derivative is specified in lVal or rVal.
DERIV2	The value of the second derivative is specified in lVal or rVal.
FLAT	The value of the first derivative is set to zero at the boundary.
NATURAL	The value of the second derivative is set to zero at the boundary.
ANTIPERIODIC	Use anti-periodic boundary conditions. The value, first derivative and second derivative at the left boundary are the negative of those at the right boundary.

Data to be interpolated

The data to be interpolated should have N_{x N_{y N_{z contiguous
elements, arranged is row-order (C-style) format. That is, the offset
of the (ix,iy,iz) element is (ix*(Ny+iy)*Nz+iz). Complex numbers are
stored in the standard format of (real,imaginary) pairs, which the
real element first.

Function prototypes:
Single-precision real:

multi_UBspline_1d_s * create_multi_UBspline_1d_s (Ugrid x_grid, BCtype_s xBC, int num_splines);
multi_UBspline_2d_s * create_multi_UBspline_2d_s (Ugrid x_grid, Ugrid y_grid,
BCtype_s xBC, BCtype_s yBC, int num_splines);
multi_UBspline_3d_s * create_multi_UBspline_3d_s (Ugrid x_grid, Ugrid y_grid, Ugrid z_grid,
BCtype_s xBC, BCtype_s yBC, BCtype_s zBC, int num_splines);

void set_multi_UBspline_1d_s (multi_UBspline_1d_s* spline, int spline_num, float *data);
void set_multi_UBspline_2d_s (multi_UBspline_2d_s* spline, int spline_num, float *data);
void set_multi_UBspline_3d_s (multi_UBspline_3d_s* spline, int spline_num, float *data);

Single-precision complex:

multi_UBspline_1d_c * create_multi_UBspline_1d_c (Ugrid x_grid, BCtype_c xBC, int num_splines);
multi_UBspline_2d_c * create_multi_UBspline_2d_c (Ugrid x_grid, Ugrid y_grid,
BCtype_c xBC, BCtype_c yBC, int num_splines);
multi_UBspline_3d_c * create_multi_UBspline_3d_c (Ugrid x_grid, Ugrid y_grid, Ugrid z_grid,
BCtype_c xBC, BCtype_c yBC, BCtype_c zBC,
int num_splines);

void set_multi_UBspline_1d_c (multi_UBspline_1d_c* spline, int spline_num, complex_float *data);
void set_multi_UBspline_2d_c (multi_UBspline_2d_c* spline, int spline_num, complex_float *data);
void set_multi_UBspline_3d_c (multi_UBspline_3d_c* spline, int spline_num, complex_float *data);

Double-precision real:

multi_UBspline_1d_d * create_multi_UBspline_1d_d (Ugrid x_grid, BCtype_d xBC, int num_splines);
multi_UBspline_2d_d * create_multi_UBspline_2d_d (Ugrid x_grid, Ugrid y_grid,
BCtype_d xBC, BCtype_d yBC, int num_splines);
multi_UBspline_3d_d * create_multi_UBspline_3d_d (Ugrid x_grid, Ugrid y_grid, Ugrid z_grid,
BCtype_d xBC, BCtype_d yBC, BCtype_d zBC, int num_splines);

void set_multi_UBspline_1d_d (multi_UBspline_1d_d* spline, int spline_num, double *data);
void set_multi_UBspline_2d_d (multi_UBspline_2d_d* spline, int spline_num, double *data);
void set_multi_UBspline_3d_d (multi_UBspline_3d_d* spline, int spline_num, double *data);

Double-precision complex:

multi_UBspline_1d_z * create_multi_UBspline_1d_z (Ugrid x_grid, BCtype_z xBC, int num_splines);
multi_UBspline_2d_z * create_multi_UBspline_2d_z (Ugrid x_grid, Ugrid y_grid,
BCtype_z xBC, BCtype_z yBC, int num_splines);
multi_UBspline_3d_z * create_multi_UBspline_3d_z (Ugrid x_grid, Ugrid y_grid, Ugrid z_grid,
BCtype_z xBC, BCtype_z yBC, BCtype_z zBC,
int num_splines);

void set_multi_UBspline_1d_z (multi_UBspline_1d_z* spline, int spline_num, complex_double *data);
void set_multi_UBspline_2d_z (multi_UBspline_2d_z* spline, int spline_num, complex_double *data);
void set_multi_UBspline_3d_z (multi_UBspline_3d_z* spline, int spline_num, complex_double *data);

Spline destruction
The memory used for spline storage can be freed simply by a call to
void
destroy_Bspline (void *spline);

The spline parameter can be a spline of any type and dimension,
uniform or nonuniform.

Spline evaluation
In contrast to the single-spline routines, eval_multi_UBspline_* routines
evaluate all num_splines splines at a given point with a single call.
The prototypes for the functions are very similar to the single-spline prototypes, but the
value, gradient, hessian, and laplacian arguments are now arrays which much be
of sufficicient size to hold the values for all the splines. If the multi-spline object
was created with N splines, the value parameter must be an array of length N.
Similarly, the gradient array must have dimension N*D, where D is the
spline dimensionality (i.e. 1, 2, or 3). Finally the hessian argument must have
dimension N*D*D.

For each of the four datatypes, there are four evaluation routines, depending on
which quantities need to be computed:

Value only

Value and gradient

Value, gradient, and Laplacian

Value, gradient, and Hessian (matrix of 2nd derivatives)

For consistency, all results are returned through pointers passed to the evaluation
routines. Currently, no bounds checking is done for the sake of speed. The user is
responsible for ensuring that the points passed to the evaluation functions fall within
the grids specified at the time of spline creation.

Function prototypes:

Quick Jump Table
1D
2D
3D

Single-precision real
1ds
2ds
3ds

Single-precision complex
1dc
2dc
3dc

Double-precision real
1dd
2dd
3dd

Double-precision complex
1dz
2dz
3dz

Single-precision real:
1D
inline void
eval_multi_UBspline_1d_s (multi_UBspline_1d_s * restrict spline,
double x, float* restrict val);

inline void
eval_multi_UBspline_1d_s_vg (multi_UBspline_1d_s * restrict spline, double x,
float* restrict val, float* restrict grad);

inline void
eval_multi_UBspline_1d_s_vgl (multi_UBspline_1d_s * restrict spline, double x,
float* restrict val, float* restrict grad, float* restrict lapl);

inline void /* identical to above routine in 1D */
eval_multi_UBspline_1d_s_vgh (multi_UBspline_1d_s * restrict spline, double x,
float* restrict val, float* restrict grad, float* restrict hess);

2D
inline void
eval_multi_UBspline_2d_s (multi_UBspline_2d_s * restrict spline, double x, double y,
float* restrict val);

inline void
eval_multi_UBspline_2d_s_vg (multi_UBspline_2d_s * restrict spline, double x, double y,
float* restrict val, float* restrict grad);

inline void
eval_multi_UBspline_2d_s_vgl (multi_UBspline_2d_s * restrict spline, double x, double y,
float* restrict val, float* restrict grad, float* restrict lapl);

inline void
eval_multi_UBspline_2d_s_vgh (multi_UBspline_2d_s * restrict spline, double x, double y,
float* restrict val, float* restrict grad, float* restrict hess);

3D
inline void
eval_multi_UBspline_3d_s (multi_UBspline_3d_s * restrict spline, double x, double y, double z,
float* restrict val);

inline void
eval_multi_UBspline_3d_s_vg (multi_UBspline_3d_s * restrict spline, double x, double y, double z,
float* restrict val, float* restrict grad);

inline void
eval_multi_UBspline_3d_s_vgl (multi_UBspline_3d_s * restrict spline, double x, double y, double z,
float* restrict val, float* restrict grad, float* restrict lapl);

inline void
eval_multi_UBspline_3d_s_vgh (multi_UBspline_3d_s * restrict spline, double x, double y,
float* restrict val, float* restrict grad, float* restrict hess);

Single-precision complex:
1D
inline void
eval_multi_UBspline_1d_c (multi_UBspline_1d_c * restrict spline,
double x, complex_float* restrict val);

inline void
eval_multi_UBspline_1d_c_vg (multi_UBspline_1d_c * restrict spline, double x,
complex_float* restrict val, complex_float* restrict grad);

inline void
eval_multi_UBspline_1d_c_vgl (multi_UBspline_1d_c * restrict spline, double x,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict lapl);

inline void /* identical to above routine in 1D */
eval_multi_UBspline_1d_c_vgh (multi_UBspline_1d_c * restrict spline, double x,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict hess);

2D
inline void
eval_multi_UBspline_2d_c (multi_UBspline_2d_c * restrict spline, double x, double y,
complex_float* restrict val);

inline void
eval_multi_UBspline_2d_c_vg (multi_UBspline_2d_c * restrict spline, double x, double y,
complex_float* restrict val, complex_float* restrict grad);

inline void
eval_multi_UBspline_2d_c_vgl (multi_UBspline_2d_c * restrict spline, double x, double y,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict lapl);

inline void
eval_multi_UBspline_2d_c_vgh (multi_UBspline_2d_c * restrict spline, double x, double y,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict hess);

3D
inline void
eval_multi_UBspline_3d_c (multi_UBspline_3d_c * restrict spline, double x, double y, double z,
complex_float* restrict val);

inline void
eval_multi_UBspline_3d_c_vg (multi_UBspline_3d_c * restrict spline, double x, double y, double z,
complex_float* restrict val, complex_float* restrict grad);

inline void
eval_multi_UBspline_3d_c_vgl (multi_UBspline_3d_c * restrict spline, double x, double y, double z,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict lapl);

inline void
eval_multi_UBspline_3d_c_vgh (multi_UBspline_3d_c * restrict spline, double x, double y,
complex_float* restrict val, complex_float* restrict grad, complex_float* restrict hess);

Double-precision real:
1D
inline void
eval_multi_UBspline_1d_d (multi_UBspline_1d_d * restrict spline,
double x, double* restrict val);

inline void
eval_multi_UBspline_1d_d_vg (multi_UBspline_1d_d * restrict spline, double x,
double* restrict val, double* restrict grad);

inline void
eval_multi_UBspline_1d_d_vgl (multi_UBspline_1d_d * restrict spline, double x,
double* restrict val, double* restrict grad, double* restrict lapl);

inline void /* identical to above routine in 1D */
eval_multi_UBspline_1d_d_vgh (multi_UBspline_1d_d * restrict spline, double x,
double* restrict val, double* restrict grad, double* restrict hess);

2D
inline void
eval_multi_UBspline_2d_d (multi_UBspline_2d_d * restrict spline, double x, double y,
double* restrict val);

inline void
eval_multi_UBspline_2d_d_vg (multi_UBspline_2d_d * restrict spline, double x, double y,
double* restrict val, double* restrict grad);

inline void
eval_multi_UBspline_2d_d_vgl (multi_UBspline_2d_d * restrict spline, double x, double y,
double* restrict val, double* restrict grad, double* restrict lapl);

inline void
eval_multi_UBspline_2d_d_vgh (multi_UBspline_2d_d * restrict spline, double x, double y,
double* restrict val, double* restrict grad, double* restrict hess);

3D
inline void
eval_multi_UBspline_3d_d (multi_UBspline_3d_d * restrict spline, double x, double y, double z,
double* restrict val);

inline void
eval_multi_UBspline_3d_d_vg (multi_UBspline_3d_d * restrict spline, double x, double y, double z,
double* restrict val, double* restrict grad);

inline void
eval_multi_UBspline_3d_d_vgl (multi_UBspline_3d_d * restrict spline, double x, double y, double z,
double* restrict val, double* restrict grad, double* restrict lapl);

inline void
eval_multi_UBspline_3d_d_vgh (multi_UBspline_3d_d * restrict spline, double x, double y,
double* restrict val, double* restrict grad, double* restrict hess);

Double-precision complex:
1D
inline void
eval_multi_UBspline_1d_z (multi_UBspline_1d_z * restrict spline,
double x, complex_double* restrict val);

inline void
eval_multi_UBspline_1d_z_vg (multi_UBspline_1d_z * restrict spline, double x,
complex_double* restrict val, complex_double* restrict grad);

inline void
eval_multi_UBspline_1d_z_vgl (multi_UBspline_1d_z * restrict spline, double x,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict lapl);

inline void /* identical to above routine in 1D */
eval_multi_UBspline_1d_z_vgh (multi_UBspline_1d_z * restrict spline, double x,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict hess);

2D
inline void
eval_multi_UBspline_2d_z (multi_UBspline_2d_z * restrict spline, double x, double y,
complex_double* restrict val);

inline void
eval_multi_UBspline_2d_z_vg (multi_UBspline_2d_z * restrict spline, double x, double y,
complex_double* restrict val, complex_double* restrict grad);

inline void
eval_multi_UBspline_2d_z_vgl (multi_UBspline_2d_z * restrict spline, double x, double y,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict lapl);

inline void
eval_multi_UBspline_2d_z_vgh (multi_UBspline_2d_z * restrict spline, double x, double y,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict hess);

3D
inline void
eval_multi_UBspline_3d_z (multi_UBspline_3d_z * restrict spline, double x, double y, double z,
complex_double* restrict val);

inline void
eval_multi_UBspline_3d_z_vg (multi_UBspline_3d_z * restrict spline, double x, double y, double z,
complex_double* restrict val, complex_double* restrict grad);

inline void
eval_multi_UBspline_3d_z_vgl (multi_UBspline_3d_z * restrict spline, double x, double y, double z,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict lapl);

inline void
eval_multi_UBspline_3d_z_vgh (multi_UBspline_3d_z * restrict spline, double x, double y,
complex_double* restrict val, complex_double* restrict grad, complex_double* restrict hess);}}}

Quick Jump Table	1D	2D	3D
Single-precision real	1ds	2ds	3ds
Single-precision complex	1dc	2dc	3dc
Double-precision real	1dd	2dd	3dd
Double-precision complex	1dz	2dz	3dz