nubspline_create.cpp File Reference

#include "nubspline_create.h"
#include <math.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>

Go to the source code of this file.

Macros
#define	__USE_XOPEN2K
#define	_XOPEN_SOURCE   600

Functions
NUBspline_1d_c *	create_NUBspline_1d_c (NUgrid *x_grid, BCtype_c xBC, complex_float *data)
NUBspline_1d_d *	create_NUBspline_1d_d (NUgrid *x_grid, BCtype_d xBC, double *data)
NUBspline_1d_s *	create_NUBspline_1d_s (NUgrid *x_grid, BCtype_s xBC, float *data)
NUBspline_1d_z *	create_NUBspline_1d_z (NUgrid *x_grid, BCtype_z xBC, complex_double *data)
NUBspline_2d_c *	create_NUBspline_2d_c (NUgrid *x_grid, NUgrid *y_grid, BCtype_c xBC, BCtype_c yBC, complex_float *data)
NUBspline_2d_d *	create_NUBspline_2d_d (NUgrid *x_grid, NUgrid *y_grid, BCtype_d xBC, BCtype_d yBC, double *data)
NUBspline_2d_s *	create_NUBspline_2d_s (NUgrid *x_grid, NUgrid *y_grid, BCtype_s xBC, BCtype_s yBC, float *data)
NUBspline_2d_z *	create_NUBspline_2d_z (NUgrid *x_grid, NUgrid *y_grid, BCtype_z xBC, BCtype_z yBC, complex_double *data)
NUBspline_3d_c *	create_NUBspline_3d_c (NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_c xBC, BCtype_c yBC, BCtype_c zBC, complex_float *data)
NUBspline_3d_d *	create_NUBspline_3d_d (NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_d xBC, BCtype_d yBC, BCtype_d zBC, double *data)
NUBspline_3d_s *	create_NUBspline_3d_s (NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_s xBC, BCtype_s yBC, BCtype_s zBC, float *data)
NUBspline_3d_z *	create_NUBspline_3d_z (NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_z xBC, BCtype_z yBC, BCtype_z zBC, complex_double *data)
void	destroy_NUBspline (Bspline *spline)
void	find_NUBcoefs_1d_c (NUBasis *restrict basis, BCtype_c bc, complex_float *data, int dstride, complex_float *coefs, int cstride)
void	find_NUBcoefs_1d_d (NUBasis *restrict basis, BCtype_d bc, double *data, int dstride, double *coefs, int cstride)
void	find_NUBcoefs_1d_s (NUBasis *restrict basis, BCtype_s bc, float *data, int dstride, float *coefs, int cstride)
void	find_NUBcoefs_1d_z (NUBasis *restrict basis, BCtype_z bc, complex_double *data, int dstride, complex_double *coefs, int cstride)
void	solve_NUB_deriv_interp_1d_d (NUBasis *restrict basis, double *restrict data, int datastride, double *restrict p, int pstride, double abcdInitial[4], double abcdFinal[4])
void	solve_NUB_deriv_interp_1d_s (NUBasis *restrict basis, float *restrict data, int datastride, float *restrict p, int pstride, float abcdInitial[4], float abcdFinal[4])
void	solve_NUB_periodic_interp_1d_d (NUBasis *restrict basis, double *restrict data, int datastride, double *restrict p, int pstride)
void	solve_NUB_periodic_interp_1d_s (NUBasis *restrict basis, float *restrict data, int datastride, float *restrict p, int pstride)

Macro Definition Documentation

__USE_XOPEN2K

#define __USE_XOPEN2K

Definition at line 28 of file nubspline_create.cpp.

_XOPEN_SOURCE

#define _XOPEN_SOURCE   600

Definition at line 25 of file nubspline_create.cpp.

Function Documentation

create_NUBspline_1d_c()

NUBspline_1d_c * create_NUBspline_1d_c	(	NUgrid *	x_grid,
		BCtype_c	xBC,
		complex_float *	data )

Definition at line 763 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_1d_c* spline = static_cast<NUBspline_1d_c*>(malloc(sizeof(NUBspline_1d_c)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU1D;
  spline->t_code  = SINGLE_COMPLEX;
 
  // Next, create the basis
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  // M is the number of data points
//  int M;
//  if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
//  else                       M = x_grid->num_points;
  int N = x_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->coefs = (complex_float*)malloc (sizeof(complex_float)*N);
  find_NUBcoefs_1d_c (spline->x_basis, xBC, data, 1, spline->coefs, 1);
    
  return spline;
}

References NUBspline_1d_c::coefs, create_NUBasis(), find_NUBcoefs_1d_c(), BCtype_c::lCode, NU1D, NUgrid::num_points, PERIODIC, SINGLE_COMPLEX, NUBspline_1d_c::sp_code, NUBspline_1d_c::t_code, and NUBspline_1d_c::x_basis.

create_NUBspline_1d_d()

NUBspline_1d_d * create_NUBspline_1d_d	(	NUgrid *	x_grid,
		BCtype_d	xBC,
		double *	data )

Definition at line 591 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_1d_d* spline = static_cast<NUBspline_1d_d*>(malloc(sizeof(NUBspline_1d_d)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU1D;
  spline->t_code  = DOUBLE_REAL;
 
  // Next, create the basis
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  // M is the number of data points (set but unused)
//  int M;
//  if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
//  else                       M = x_grid->num_points;
  int N = x_grid->num_points + 2;
 
  // Allocate coefficients and solve
  spline->coefs = (double*)malloc (sizeof(double)*N);
  find_NUBcoefs_1d_d (spline->x_basis, xBC, data, 1, spline->coefs, 1);
    
  return spline;
}

References NUBspline_1d_d::coefs, create_NUBasis(), DOUBLE_REAL, find_NUBcoefs_1d_d(), BCtype_d::lCode, NU1D, NUgrid::num_points, PERIODIC, NUBspline_1d_d::sp_code, NUBspline_1d_d::t_code, and NUBspline_1d_d::x_basis.

create_NUBspline_1d_s()

NUBspline_1d_s * create_NUBspline_1d_s	(	NUgrid *	x_grid,
		BCtype_s	xBC,
		float *	data )

Definition at line 247 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_1d_s* spline = static_cast<NUBspline_1d_s*>(malloc(sizeof(NUBspline_1d_s)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU1D;
  spline->t_code  = SINGLE_REAL;
 
  // Next, create the basis
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  // M is the number of data points (but is unused)
//  int M;
//  if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
//  else                       M = x_grid->num_points;
  int N = x_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->coefs = (float*)malloc (sizeof(float)*N);
  find_NUBcoefs_1d_s (spline->x_basis, xBC, data, 1, spline->coefs, 1);
    
  return spline;
}

References NUBspline_1d_s::coefs, create_NUBasis(), find_NUBcoefs_1d_s(), BCtype_s::lCode, NU1D, NUgrid::num_points, PERIODIC, SINGLE_REAL, NUBspline_1d_s::sp_code, NUBspline_1d_s::t_code, and NUBspline_1d_s::x_basis.

create_NUBspline_1d_z()

NUBspline_1d_z * create_NUBspline_1d_z	(	NUgrid *	x_grid,
		BCtype_z	xBC,
		complex_double *	data )

Definition at line 932 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_1d_z* spline = static_cast<NUBspline_1d_z*>(malloc(sizeof(NUBspline_1d_z)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU1D;
  spline->t_code  = DOUBLE_COMPLEX;
 
  // Next, create the basis
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  // M is the number of data points
//  int M;
//  if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
//  else                       M = x_grid->num_points;
  int N = x_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->coefs = (complex_double*)malloc (sizeof(complex_double)*N);
  find_NUBcoefs_1d_z (spline->x_basis, xBC, data, 1, spline->coefs, 1);
    
  return spline;
}

References NUBspline_1d_z::coefs, create_NUBasis(), DOUBLE_COMPLEX, find_NUBcoefs_1d_z(), BCtype_z::lCode, NU1D, NUgrid::num_points, PERIODIC, NUBspline_1d_z::sp_code, NUBspline_1d_z::t_code, and NUBspline_1d_z::x_basis.

create_NUBspline_2d_c()

NUBspline_2d_c * create_NUBspline_2d_c	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		BCtype_c	xBC,
		BCtype_c	yBC,
		complex_float *	data )

Definition at line 788 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_2d_c* spline = static_cast<NUBspline_2d_c*>(malloc(sizeof(NUBspline_2d_c)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU2D;
  spline->t_code  = SINGLE_COMPLEX;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
//  int Mx,
  int My, Nx, Ny;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
    
  spline->x_stride = Ny;
#ifndef HAVE_SSE2
  spline->coefs = (complex_float*)malloc (sizeof(complex_float)*Nx*Ny);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_float)*Nx*Ny);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    int doffset = iy;
    int coffset = iy;
    find_NUBcoefs_1d_c (spline->x_basis, xBC, data+doffset, My,
            spline->coefs+coffset, Ny);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    int doffset = ix*Ny;
    int coffset = ix*Ny;
    find_NUBcoefs_1d_c (spline->y_basis, yBC, spline->coefs+doffset, 1, 
            spline->coefs+coffset, 1);
  }
    
  return spline;
}

References NUBspline_2d_c::coefs, create_NUBasis(), find_NUBcoefs_1d_c(), BCtype_c::lCode, NU2D, NUgrid::num_points, PERIODIC, posix_memalign(), SINGLE_COMPLEX, NUBspline_2d_c::sp_code, NUBspline_2d_c::t_code, NUBspline_2d_c::x_basis, NUBspline_2d_c::x_stride, and NUBspline_2d_c::y_basis.

create_NUBspline_2d_d()

NUBspline_2d_d * create_NUBspline_2d_d	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		BCtype_d	xBC,
		BCtype_d	yBC,
		double *	data )

Definition at line 616 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_2d_d* spline = static_cast<NUBspline_2d_d*>(malloc(sizeof(NUBspline_2d_d)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU2D;
  spline->t_code  = DOUBLE_REAL;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
 
  // int Mx, (set but unused)
  int My, Nx, Ny;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
  
  spline->x_stride = Ny;
#ifndef HAVE_SSE2
  spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(double)*Nx*Ny);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    int doffset = iy;
    int coffset = iy;
    find_NUBcoefs_1d_d (spline->x_basis, xBC, data+doffset, My,
            spline->coefs+coffset, Ny);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    int doffset = ix*Ny;
    int coffset = ix*Ny;
    find_NUBcoefs_1d_d (spline->y_basis, yBC, spline->coefs+doffset, 1, 
            spline->coefs+coffset, 1);
  }
    
  return spline;
}

References NUBspline_2d_d::coefs, create_NUBasis(), DOUBLE_REAL, find_NUBcoefs_1d_d(), BCtype_d::lCode, NU2D, NUgrid::num_points, PERIODIC, posix_memalign(), NUBspline_2d_d::sp_code, NUBspline_2d_d::t_code, NUBspline_2d_d::x_basis, NUBspline_2d_d::x_stride, and NUBspline_2d_d::y_basis.

create_NUBspline_2d_s()

NUBspline_2d_s * create_NUBspline_2d_s	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		BCtype_s	xBC,
		BCtype_s	yBC,
		float *	data )

Definition at line 272 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_2d_s* spline = static_cast<NUBspline_2d_s*>(malloc(sizeof(NUBspline_2d_s)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU2D;
  spline->t_code  = SINGLE_REAL;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
  // set but unused
  //int Mx,
  int My, Nx, Ny;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
    
  spline->x_stride = Ny;
#ifndef HAVE_SSE2
  spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(float)*Nx*Ny);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    int doffset = iy;
    int coffset = iy;
    find_NUBcoefs_1d_s (spline->x_basis, xBC, data+doffset, My,
            spline->coefs+coffset, Ny);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    int doffset = ix*Ny;
    int coffset = ix*Ny;
    find_NUBcoefs_1d_s (spline->y_basis, yBC, spline->coefs+doffset, 1, 
            spline->coefs+coffset, 1);
  }
    
  return spline;
}

References NUBspline_2d_s::coefs, create_NUBasis(), find_NUBcoefs_1d_s(), BCtype_s::lCode, NU2D, NUgrid::num_points, PERIODIC, posix_memalign(), SINGLE_REAL, NUBspline_2d_s::sp_code, NUBspline_2d_s::t_code, NUBspline_2d_s::x_basis, NUBspline_2d_s::x_stride, and NUBspline_2d_s::y_basis.

create_NUBspline_2d_z()

NUBspline_2d_z * create_NUBspline_2d_z	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		BCtype_z	xBC,
		BCtype_z	yBC,
		complex_double *	data )

Definition at line 957 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_2d_z* spline = static_cast<NUBspline_2d_z*>(malloc(sizeof(NUBspline_2d_z)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU2D;
  spline->t_code  = DOUBLE_COMPLEX;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
//  int Mx,
  int My, Nx, Ny;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
    
  spline->x_stride = Ny;
#ifndef HAVE_SSE2
  spline->coefs = (complex_double*)malloc (sizeof(complex_double)*Nx*Ny);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_double)*Nx*Ny);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    int doffset = iy;
    int coffset = iy;
    find_NUBcoefs_1d_z (spline->x_basis, xBC, data+doffset, My,
            spline->coefs+coffset, Ny);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    int doffset = ix*Ny;
    int coffset = ix*Ny;
    find_NUBcoefs_1d_z (spline->y_basis, yBC, spline->coefs+doffset, 1, 
            spline->coefs+coffset, 1);
  }
    
  return spline;
}

References NUBspline_2d_z::coefs, create_NUBasis(), DOUBLE_COMPLEX, find_NUBcoefs_1d_z(), BCtype_z::lCode, NU2D, NUgrid::num_points, PERIODIC, posix_memalign(), NUBspline_2d_z::sp_code, NUBspline_2d_z::t_code, NUBspline_2d_z::x_basis, NUBspline_2d_z::x_stride, and NUBspline_2d_z::y_basis.

create_NUBspline_3d_c()

NUBspline_3d_c * create_NUBspline_3d_c	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		NUgrid *	z_grid,
		BCtype_c	xBC,
		BCtype_c	yBC,
		BCtype_c	zBC,
		complex_float *	data )

Definition at line 839 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_3d_c* spline = static_cast<NUBspline_3d_c*>(malloc(sizeof(NUBspline_3d_c)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU3D;
  spline->t_code  = SINGLE_COMPLEX;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
  spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
//  int Mx,
  int My, Mz, Nx, Ny, Nz;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
  if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
  else                       Mz = z_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
  Nz = z_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->x_stride = Ny*Nz;
  spline->y_stride = Nz;
#ifndef HAVE_SSE2
  spline->coefs = (complex_float*)malloc (sizeof(complex_float)*Nx*Ny*Nz);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_float)*Nx*Ny*Nz);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      int doffset = iy*Mz+iz;
      int coffset = iy*Nz+iz;
      find_NUBcoefs_1d_c (spline->x_basis, xBC, data+doffset, My*Mz,
              spline->coefs+coffset, Ny*Nz);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      int doffset = ix*Ny*Nz + iz;
      int coffset = ix*Ny*Nz + iz;
      find_NUBcoefs_1d_c (spline->y_basis, yBC, spline->coefs+doffset, Nz, 
              spline->coefs+coffset, Nz);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      int doffset = (ix*Ny+iy)*Nz;
      int coffset = (ix*Ny+iy)*Nz;
      find_NUBcoefs_1d_c (spline->z_basis, zBC, spline->coefs+doffset, 1, 
              spline->coefs+coffset, 1);
    }
  return spline;
}

References NUBspline_3d_c::coefs, create_NUBasis(), find_NUBcoefs_1d_c(), BCtype_c::lCode, NU3D, NUgrid::num_points, PERIODIC, posix_memalign(), SINGLE_COMPLEX, NUBspline_3d_c::sp_code, NUBspline_3d_c::t_code, NUBspline_3d_c::x_basis, NUBspline_3d_c::x_stride, NUBspline_3d_c::y_basis, NUBspline_3d_c::y_stride, and NUBspline_3d_c::z_basis.

create_NUBspline_3d_d()

NUBspline_3d_d * create_NUBspline_3d_d	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		NUgrid *	z_grid,
		BCtype_d	xBC,
		BCtype_d	yBC,
		BCtype_d	zBC,
		double *	data )

Definition at line 668 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_3d_d* spline = static_cast<NUBspline_3d_d*>(malloc(sizeof(NUBspline_3d_d)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU3D;
  spline->t_code  = DOUBLE_REAL;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
  spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
 
  // set but unused
  //int Mx,
  int My, Mz, Nx, Ny, Nz;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
  if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
  else                       Mz = z_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
  Nz = z_grid->num_points + 2;
  
  spline->x_stride = Ny*Nz;
  spline->y_stride = Nz;
#ifndef HAVE_SSE2
  spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny*Nz);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(double)*Nx*Ny*Nz);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      int doffset = iy*Mz+iz;
      int coffset = iy*Nz+iz;
      find_NUBcoefs_1d_d (spline->x_basis, xBC, data+doffset, My*Mz,
              spline->coefs+coffset, Ny*Nz);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      int doffset = ix*Ny*Nz + iz;
      int coffset = ix*Ny*Nz + iz;
      find_NUBcoefs_1d_d (spline->y_basis, yBC, spline->coefs+doffset, Nz, 
              spline->coefs+coffset, Nz);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      int doffset = (ix*Ny+iy)*Nz;
      int coffset = (ix*Ny+iy)*Nz;
      find_NUBcoefs_1d_d (spline->z_basis, zBC, spline->coefs+doffset, 1, 
              spline->coefs+coffset, 1);
    }
  return spline;
}

References NUBspline_3d_d::coefs, create_NUBasis(), DOUBLE_REAL, find_NUBcoefs_1d_d(), BCtype_d::lCode, NU3D, NUgrid::num_points, PERIODIC, posix_memalign(), NUBspline_3d_d::sp_code, NUBspline_3d_d::t_code, NUBspline_3d_d::x_basis, NUBspline_3d_d::x_stride, NUBspline_3d_d::y_basis, NUBspline_3d_d::y_stride, and NUBspline_3d_d::z_basis.

create_NUBspline_3d_s()

NUBspline_3d_s * create_NUBspline_3d_s	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		NUgrid *	z_grid,
		BCtype_s	xBC,
		BCtype_s	yBC,
		BCtype_s	zBC,
		float *	data )

Definition at line 324 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_3d_s* spline = static_cast<NUBspline_3d_s*>(malloc(sizeof(NUBspline_3d_s)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU3D;
  spline->t_code  = SINGLE_REAL;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
  spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
  // set but unused
  // int Mx,
  int My, Mz, Nx, Ny, Nz;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
  if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
  else                       Mz = z_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
  Nz = z_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->x_stride = Ny*Nz;
  spline->y_stride = Nz;
#ifndef HAVE_SSE2
  spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny*Nz);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(float)*Nx*Ny*Nz);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      int doffset = iy*Mz+iz;
      int coffset = iy*Nz+iz;
      find_NUBcoefs_1d_s (spline->x_basis, xBC, data+doffset, My*Mz,
              spline->coefs+coffset, Ny*Nz);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      int doffset = ix*Ny*Nz + iz;
      int coffset = ix*Ny*Nz + iz;
      find_NUBcoefs_1d_s (spline->y_basis, yBC, spline->coefs+doffset, Nz, 
              spline->coefs+coffset, Nz);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      int doffset = (ix*Ny+iy)*Nz;
      int coffset = (ix*Ny+iy)*Nz;
      find_NUBcoefs_1d_s (spline->z_basis, zBC, spline->coefs+doffset, 1, 
              spline->coefs+coffset, 1);
    }
  return spline;
}

References NUBspline_3d_s::coefs, create_NUBasis(), find_NUBcoefs_1d_s(), BCtype_s::lCode, NU3D, NUgrid::num_points, PERIODIC, posix_memalign(), SINGLE_REAL, NUBspline_3d_s::sp_code, NUBspline_3d_s::t_code, NUBspline_3d_s::x_basis, NUBspline_3d_s::x_stride, NUBspline_3d_s::y_basis, NUBspline_3d_s::y_stride, and NUBspline_3d_s::z_basis.

create_NUBspline_3d_z()

NUBspline_3d_z * create_NUBspline_3d_z	(	NUgrid *	x_grid,
		NUgrid *	y_grid,
		NUgrid *	z_grid,
		BCtype_z	xBC,
		BCtype_z	yBC,
		BCtype_z	zBC,
		complex_double *	data )

Definition at line 1008 of file nubspline_create.cpp.

{
  // First, create the spline structure
  NUBspline_3d_z* spline = static_cast<NUBspline_3d_z*>(malloc(sizeof(NUBspline_3d_z)));
  if (spline == NULL)
    return spline;
  spline->sp_code = NU3D;
  spline->t_code  = DOUBLE_COMPLEX;
  spline->x_grid = x_grid;
  spline->y_grid = y_grid;
  spline->z_grid = z_grid;
 
  // Next, create the bases
  spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
  spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
  spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
//  int Mx,
  int My, Mz, Nx, Ny, Nz;
//  if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
//  else                       Mx = x_grid->num_points;
  if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
  else                       My = y_grid->num_points;
  if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
  else                       Mz = z_grid->num_points;
 
  Nx = x_grid->num_points + 2;
  Ny = y_grid->num_points + 2;
  Nz = z_grid->num_points + 2;
 
  // Allocate coefficients and solve  
  spline->x_stride = Ny*Nz;
  spline->y_stride = Nz;
#ifndef HAVE_SSE2
  spline->coefs = (complex_double*)malloc (sizeof(complex_double)*Nx*Ny*Nz);
#else
  posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_double)*Nx*Ny*Nz);
#endif
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      int doffset = iy*Mz+iz;
      int coffset = iy*Nz+iz;
      find_NUBcoefs_1d_z (spline->x_basis, xBC, data+doffset, My*Mz,
              spline->coefs+coffset, Ny*Nz);
      /* for (int ix=0; ix<Nx; ix++) {
    complex_double z = spline->coefs[coffset+ix*spline->x_stride];
    if (isnan(creal(z)))
      fprintf (stderr, "NAN encountered in create_NUBspline_3d_z at real part of (%d,%d,%d)\n",
           ix,iy,iz);
    if (isnan(cimag(z)))
      fprintf (stderr, "NAN encountered in create_NUBspline_3d_z at imag part of (%d,%d,%d)\n",
           ix,iy,iz);
       } */
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      int doffset = ix*Ny*Nz + iz;
      int coffset = ix*Ny*Nz + iz;
      find_NUBcoefs_1d_z (spline->y_basis, yBC, spline->coefs+doffset, Nz, 
              spline->coefs+coffset, Nz);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      int doffset = (ix*Ny+iy)*Nz;
      int coffset = (ix*Ny+iy)*Nz;
      find_NUBcoefs_1d_z (spline->z_basis, zBC, spline->coefs+doffset, 1, 
              spline->coefs+coffset, 1);
    }
  return spline;
}

References NUBspline_3d_z::coefs, create_NUBasis(), DOUBLE_COMPLEX, find_NUBcoefs_1d_z(), BCtype_z::lCode, NU3D, NUgrid::num_points, PERIODIC, posix_memalign(), NUBspline_3d_z::sp_code, NUBspline_3d_z::t_code, NUBspline_3d_z::x_basis, NUBspline_3d_z::x_grid, NUBspline_3d_z::x_stride, NUBspline_3d_z::y_basis, NUBspline_3d_z::y_grid, NUBspline_3d_z::y_stride, NUBspline_3d_z::z_basis, and NUBspline_3d_z::z_grid.

destroy_NUBspline()

void destroy_NUBspline

(

Bspline *

spline

)

Definition at line 1087 of file nubspline_create.cpp.

{
    free(spline->coefs);
  switch (spline->sp_code) {
  case NU1D:
    destroy_NUBasis (((NUBspline_1d*)spline)->x_basis);
    break;
  case NU2D:
    destroy_NUBasis (((NUBspline_2d*)spline)->x_basis);
    destroy_NUBasis (((NUBspline_2d*)spline)->y_basis);
    break;
    
  case NU3D:
    destroy_NUBasis (((NUBspline_3d*)spline)->x_basis);
    destroy_NUBasis (((NUBspline_3d*)spline)->y_basis);
    destroy_NUBasis (((NUBspline_3d*)spline)->z_basis);
    break;
  case U1D:
  case U2D:
  case MULTI_U1D:
  case MULTI_U2D:
  case MULTI_U3D:
  case MULTI_NU1D:
  case MULTI_NU2D:
  case MULTI_NU3D:
  default:
      ;
  }
  free(spline);
}

References Bspline::coefs, destroy_NUBasis(), MULTI_NU1D, MULTI_NU2D, MULTI_NU3D, MULTI_U1D, MULTI_U2D, MULTI_U3D, NU1D, NU2D, NU3D, Bspline::sp_code, U1D, and U2D.

find_NUBcoefs_1d_c()

void find_NUBcoefs_1d_c	(	NUBasis *restrict	basis,
		BCtype_c	bc,
		complex_float *	data,
		int	dstride,
		complex_float *	coefs,
		int	cstride )

Definition at line 742 of file nubspline_create.cpp.

{
  BCtype_s bc_r, bc_i;
  bc_r.lCode = bc.lCode;   bc_i.lCode = bc.lCode;
  bc_r.rCode = bc.rCode;   bc_i.rCode = bc.rCode;
  bc_r.lVal  = bc.lVal_r;  bc_r.rVal  = bc.rVal_r;
  bc_i.lVal  = bc.lVal_i;  bc_i.rVal  = bc.rVal_i;
 
  float *data_r  = ((float*)data );
  float *data_i  = ((float*)data )+1;
  float *coefs_r = ((float*)coefs);
  float *coefs_i = ((float*)coefs)+1;
 
  find_NUBcoefs_1d_s (basis, bc_r, data_r, 2*dstride, coefs_r, 2*cstride);
  find_NUBcoefs_1d_s (basis, bc_i, data_i, 2*dstride, coefs_i, 2*cstride);
}

References find_NUBcoefs_1d_s(), BCtype_s::lCode, BCtype_c::lCode, BCtype_s::lVal, BCtype_c::lVal_i, BCtype_c::lVal_r, BCtype_s::rCode, BCtype_c::rCode, BCtype_s::rVal, BCtype_c::rVal_i, and BCtype_c::rVal_r.

find_NUBcoefs_1d_d()

void find_NUBcoefs_1d_d	(	NUBasis *restrict	basis,
		BCtype_d	bc,
		double *	data,
		int	dstride,
		double *	coefs,
		int	cstride )

Definition at line 543 of file nubspline_create.cpp.

{
  if (bc.lCode == PERIODIC) 
    solve_NUB_periodic_interp_1d_d (basis, data, dstride, coefs, cstride);
  else {
    int M = basis->grid->num_points;
    // Setup boundary conditions
    double bfuncs[4] {};
    double dbfuncs[4] {};
    double abcd_left[4] {};
    double abcd_right[4] {};
 
    // Left boundary
    if (bc.lCode == FLAT || bc.lCode == NATURAL)
      bc.lVal = 0.0;
    if (bc.lCode == FLAT || bc.lCode == DERIV1) {
      get_NUBasis_dfuncs_di (basis, 0, bfuncs, abcd_left);
      abcd_left[3] = bc.lVal;
    }
    if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
      get_NUBasis_d2funcs_di (basis, 0, bfuncs, dbfuncs, abcd_left);
      abcd_left[3] = bc.lVal;
    }
    
    // Right boundary
    if (bc.rCode == FLAT || bc.rCode == NATURAL)
      bc.rVal = 0.0;
    if (bc.rCode == FLAT || bc.rCode == DERIV1) {
      get_NUBasis_dfuncs_di (basis, M-1, bfuncs, abcd_right);
      abcd_right[3] = bc.rVal;
    }
    if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
      get_NUBasis_d2funcs_di (basis, M-1, bfuncs, dbfuncs, abcd_right);
      abcd_right[3] = bc.rVal;
    }
 
    // Now, solve for coefficients
    solve_NUB_deriv_interp_1d_d (basis, data, dstride, coefs, cstride,
                 abcd_left, abcd_right);
  }
}

References DERIV1, DERIV2, FLAT, get_NUBasis_d2funcs_di(), get_NUBasis_dfuncs_di(), BCtype_d::lCode, BCtype_d::lVal, NATURAL, PERIODIC, BCtype_d::rCode, BCtype_d::rVal, solve_NUB_deriv_interp_1d_d(), and solve_NUB_periodic_interp_1d_d().

find_NUBcoefs_1d_s()

void find_NUBcoefs_1d_s	(	NUBasis *restrict	basis,
		BCtype_s	bc,
		float *	data,
		int	dstride,
		float *	coefs,
		int	cstride )

Definition at line 201 of file nubspline_create.cpp.

{
  if (bc.lCode == PERIODIC) 
    solve_NUB_periodic_interp_1d_s (basis, data, dstride, coefs, cstride);
  else {
    int M = basis->grid->num_points;
    // Setup boundary conditions
    float bfuncs[4] = {};
    float dbfuncs[4] = {};
    float abcd_left[4] = {};
    float abcd_right[4] = {};
    // Left boundary
    if (bc.lCode == FLAT || bc.lCode == NATURAL)
      bc.lVal = 0.0;
    if (bc.lCode == FLAT || bc.lCode == DERIV1) {
      get_NUBasis_dfuncs_si (basis, 0, bfuncs, abcd_left);
      abcd_left[3] = bc.lVal;
    }
    if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
      get_NUBasis_d2funcs_si (basis, 0, bfuncs, dbfuncs, abcd_left);
      abcd_left[3] = bc.lVal;
    }
    
    // Right boundary
    if (bc.rCode == FLAT || bc.rCode == NATURAL)
      bc.rVal = 0.0;
    if (bc.rCode == FLAT || bc.rCode == DERIV1) {
      get_NUBasis_dfuncs_si (basis, M-1, bfuncs, abcd_right);
      abcd_right[3] = bc.rVal;
    }
    if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
      get_NUBasis_d2funcs_si (basis, M-1, bfuncs, dbfuncs, abcd_right);
      abcd_right[3] = bc.rVal;
    }
    // Now, solve for coefficients
    solve_NUB_deriv_interp_1d_s (basis, data, dstride, coefs, cstride,
                 abcd_left, abcd_right);
  }
}

References DERIV1, DERIV2, FLAT, get_NUBasis_d2funcs_si(), get_NUBasis_dfuncs_si(), BCtype_s::lCode, BCtype_s::lVal, NATURAL, PERIODIC, BCtype_s::rCode, BCtype_s::rVal, solve_NUB_deriv_interp_1d_s(), and solve_NUB_periodic_interp_1d_s().

find_NUBcoefs_1d_z()

void find_NUBcoefs_1d_z	(	NUBasis *restrict	basis,
		BCtype_z	bc,
		complex_double *	data,
		int	dstride,
		complex_double *	coefs,
		int	cstride )

Definition at line 911 of file nubspline_create.cpp.

{
  BCtype_d bc_r, bc_i;
  bc_r.lCode = bc.lCode;   bc_i.lCode = bc.lCode;
  bc_r.rCode = bc.rCode;   bc_i.rCode = bc.rCode;
  bc_r.lVal  = bc.lVal_r;  bc_r.rVal  = bc.rVal_r;
  bc_i.lVal  = bc.lVal_i;  bc_i.rVal  = bc.rVal_i;
 
  double *data_r  = ((double*)data );
  double *data_i  = ((double*)data )+1;
  double *coefs_r = ((double*)coefs);
  double *coefs_i = ((double*)coefs)+1;
 
  find_NUBcoefs_1d_d (basis, bc_r, data_r, 2*dstride, coefs_r, 2*cstride);
  find_NUBcoefs_1d_d (basis, bc_i, data_i, 2*dstride, coefs_i, 2*cstride);
}

References find_NUBcoefs_1d_d(), BCtype_d::lCode, BCtype_z::lCode, BCtype_d::lVal, BCtype_z::lVal_i, BCtype_z::lVal_r, BCtype_d::rCode, BCtype_z::rCode, BCtype_d::rVal, BCtype_z::rVal_i, and BCtype_z::rVal_r.

solve_NUB_deriv_interp_1d_d()

void solve_NUB_deriv_interp_1d_d	(	NUBasis *restrict	basis,
		double *restrict	data,
		int	datastride,
		double *restrict	p,
		int	pstride,
		double	abcdInitial[4],
		double	abcdFinal[4] )

Definition at line 396 of file nubspline_create.cpp.

{
  int M = basis->grid->num_points;
  int N = M+2;
  // Banded matrix storage.  The first three elements in the
  // tinyvector store the tridiagonal coefficients.  The last element
  // stores the RHS data.
#ifdef HAVE_C_VARARRAYS
  double bands[4*N];
#else
  double *bands = new double[4*N];
#endif
 
  // Fill up bands
  for (int i=0; i<4; i++) {
    bands[i]         = abcdInitial[i];
    bands[4*(N-1)+i] = abcdFinal[i];
  }
  for (int i=0; i<M; i++) {
    get_NUBasis_funcs_di (basis, i, &(bands[4*(i+1)]));
    bands[4*(i+1)+3] = data[datastride*i];
  }
  /* for (int i=0; i<4*N; i++)
     if (isnan(bands[i]))
     fprintf(stderr, "NAN at bands[%d].\n", i); */
    
  // Now solve:
  // First and last rows are different
  bands[4*0+1] /= bands[4*0+0];
  bands[4*0+2] /= bands[4*0+0];
  bands[4*0+3] /= bands[4*0+0];
  bands[4*0+0] = 1.0;
  bands[4*1+1] -= bands[4*1+0]*bands[4*0+1];
  bands[4*1+2] -= bands[4*1+0]*bands[4*0+2];
  bands[4*1+3] -= bands[4*1+0]*bands[4*0+3];
  bands[4*0+0] = 0.0;
  bands[4*1+2] /= bands[4*1+1];
  bands[4*1+3] /= bands[4*1+1];
  bands[4*1+1] = 1.0;
  
  // Now do rows 2 through M+1
  for (int row=2; row < N-1; row++) {
    bands[4*(row)+1] -= bands[4*(row)+0]*bands[4*(row-1)+2];
    bands[4*(row)+3] -= bands[4*(row)+0]*bands[4*(row-1)+3];
    bands[4*(row)+2] /= bands[4*(row)+1];
    bands[4*(row)+3] /= bands[4*(row)+1];
    bands[4*(row)+0] = 0.0;
    bands[4*(row)+1] = 1.0;
  }
 
  // Do last row
  bands[4*(M+1)+1] -= bands[4*(M+1)+0]*bands[4*(M-1)+2];
  bands[4*(M+1)+3] -= bands[4*(M+1)+0]*bands[4*(M-1)+3];
  bands[4*(M+1)+2] -= bands[4*(M+1)+1]*bands[4*(M)+2];
  bands[4*(M+1)+3] -= bands[4*(M+1)+1]*bands[4*(M)+3];
  bands[4*(M+1)+3] /= bands[4*(M+1)+2];
  bands[4*(M+1)+2] = 1.0;
 
  p[pstride*(M+1)] = bands[4*(M+1)+3];
  // Now back substitute up
  for (int row=M; row>0; row--)
    p[pstride*(row)] = bands[4*(row)+3] - bands[4*(row)+2]*p[pstride*(row+1)];
  
  // Finish with first row
  p[0] = bands[4*(0)+3] - bands[4*(0)+1]*p[pstride*1] - bands[4*(0)+2]*p[pstride*2];
#ifndef HAVE_C_VARARRAYS
  delete[] bands;
#endif
}

References get_NUBasis_funcs_di(), and p.

solve_NUB_deriv_interp_1d_s()

void solve_NUB_deriv_interp_1d_s	(	NUBasis *restrict	basis,
		float *restrict	data,
		int	datastride,
		float *restrict	p,
		int	pstride,
		float	abcdInitial[4],
		float	abcdFinal[4] )

Definition at line 52 of file nubspline_create.cpp.

{
  int M = basis->grid->num_points;
  int N = M+2;
  // Banded matrix storage.  The first three elements in the
  // tinyvector store the tridiagonal coefficients.  The last element
  // stores the RHS data.
#ifdef HAVE_C_VARARRAYS
  float bands[4*N];
#else
  float *bands = new float[4*N];
#endif
  
 
  // Fill up bands
  for (int i=0; i<4; i++) {
    bands[i]         = abcdInitial[i];
    bands[4*(N-1)+i] = abcdFinal[i];
  }
  for (int i=0; i<M; i++) {
    get_NUBasis_funcs_si (basis, i, &(bands[4*(i+1)]));
    bands[4*(i+1)+3] = data[datastride*i];
  }
    
  // Now solve:
  // First and last rows are different
  bands[4*0+1] /= bands[4*0+0];
  bands[4*0+2] /= bands[4*0+0];
  bands[4*0+3] /= bands[4*0+0];
  bands[4*0+0] = 1.0;
  bands[4*1+1] -= bands[4*1+0]*bands[4*0+1];
  bands[4*1+2] -= bands[4*1+0]*bands[4*0+2];
  bands[4*1+3] -= bands[4*1+0]*bands[4*0+3];
  bands[4*0+0] = 0.0;
  bands[4*1+2] /= bands[4*1+1];
  bands[4*1+3] /= bands[4*1+1];
  bands[4*1+1] = 1.0;
  
  // Now do rows 2 through M+1
  for (int row=2; row < N-1; row++) {
    bands[4*(row)+1] -= bands[4*(row)+0]*bands[4*(row-1)+2];
    bands[4*(row)+3] -= bands[4*(row)+0]*bands[4*(row-1)+3];
    bands[4*(row)+2] /= bands[4*(row)+1];
    bands[4*(row)+3] /= bands[4*(row)+1];
    bands[4*(row)+0] = 0.0;
    bands[4*(row)+1] = 1.0;
  }
 
  // Do last row
  bands[4*(M+1)+1] -= bands[4*(M+1)+0]*bands[4*(M-1)+2];
  bands[4*(M+1)+3] -= bands[4*(M+1)+0]*bands[4*(M-1)+3];
  bands[4*(M+1)+2] -= bands[4*(M+1)+1]*bands[4*(M)+2];
  bands[4*(M+1)+3] -= bands[4*(M+1)+1]*bands[4*(M)+3];
  bands[4*(M+1)+3] /= bands[4*(M+1)+2];
  bands[4*(M+1)+2] = 1.0;
 
  p[pstride*(M+1)] = bands[4*(M+1)+3];
  // Now back substitute up
  for (int row=M; row>0; row--)
    p[pstride*(row)] = bands[4*(row)+3] - bands[4*(row)+2]*p[pstride*(row+1)];
  
  // Finish with first row
  p[0] = bands[4*(0)+3] - bands[4*(0)+1]*p[pstride*1] - bands[4*(0)+2]*p[pstride*2];
 
#ifndef HAVE_C_VARARRAYS
  delete[] bands;
#endif
}

References get_NUBasis_funcs_si(), and p.

solve_NUB_periodic_interp_1d_d()

void solve_NUB_periodic_interp_1d_d	(	NUBasis *restrict	basis,
		double *restrict	data,
		int	datastride,
		double *restrict	p,
		int	pstride )

Definition at line 471 of file nubspline_create.cpp.

{
  int M = basis->grid->num_points-1;
 
  // Banded matrix storage.  The first three elements in the
  // tinyvector store the tridiagonal coefficients.  The last element
  // stores the RHS data.
#ifdef HAVE_C_VARARRAYS
  double bands[4*M], lastCol[M];
#else
  double *bands   = new double[4*M];
  double *lastCol = new double[  M];
#endif
 
  // Fill up bands
  for (int i=0; i<M; i++) {
    get_NUBasis_funcs_di (basis, i, &(bands[4*i])); 
    bands[4*(i)+3] = data[i*datastride];
  }
    
  // Now solve:
  // First and last rows are different
  bands[4*(0)+2] /= bands[4*(0)+1];
  bands[4*(0)+0] /= bands[4*(0)+1];
  bands[4*(0)+3] /= bands[4*(0)+1];
  bands[4*(0)+1]  = 1.0;
  bands[4*(M-1)+1] -= bands[4*(M-1)+2]*bands[4*(0)+0];
  bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(0)+3];
  bands[4*(M-1)+2]  = -bands[4*(M-1)+2]*bands[4*(0)+2];
  lastCol[0] = bands[4*(0)+0];
  
  for (int row=1; row < (M-1); row++) {
    bands[4*(row)+1] -= bands[4*(row)+0] * bands[4*(row-1)+2];
    bands[4*(row)+3] -= bands[4*(row)+0] * bands[4*(row-1)+3];
    lastCol[row]   = -bands[4*(row)+0] * lastCol[row-1];
    bands[4*(row)+0] = 0.0;
    bands[4*(row)+2] /= bands[4*(row)+1];
    bands[4*(row)+3] /= bands[4*(row)+1];
    lastCol[row]  /= bands[4*(row)+1];
    bands[4*(row)+1]  = 1.0;
    if (row < (M-2)) {
      bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(row)+3];
      bands[4*(M-1)+1] -= bands[4*(M-1)+2]*lastCol[row];
      bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(row)+2];
    }
  }
  
  // Now do last row
  // The [2] element and [0] element are now on top of each other 
  bands[4*(M-1)+0] += bands[4*(M-1)+2];
  bands[4*(M-1)+1] -= bands[4*(M-1)+0] * (bands[4*(M-2)+2]+lastCol[M-2]);
  bands[4*(M-1)+3] -= bands[4*(M-1)+0] *  bands[4*(M-2)+3];
  bands[4*(M-1)+3] /= bands[4*(M-1)+1];
  p[pstride*M] = bands[4*(M-1)+3];
  for (int row=M-2; row>=0; row--) 
    p[pstride*(row+1)] = bands[4*(row)+3] - 
      bands[4*(row)+2]*p[pstride*(row+2)] - lastCol[row]*p[pstride*M];
  
  p[pstride*  0  ] = p[pstride*M];
  p[pstride*(M+1)] = p[pstride*1];
  p[pstride*(M+2)] = p[pstride*2];
#ifndef HAVE_C_VARARRAYS
  delete[] bands;
  delete[] lastCol;
#endif
}

References get_NUBasis_funcs_di(), and p.

solve_NUB_periodic_interp_1d_s()

void solve_NUB_periodic_interp_1d_s	(	NUBasis *restrict	basis,
		float *restrict	data,
		int	datastride,
		float *restrict	p,
		int	pstride )

Definition at line 129 of file nubspline_create.cpp.

{
  int M = basis->grid->num_points-1;
 
  // Banded matrix storage.  The first three elements in each row
  // store the tridiagonal coefficients.  The last element
  // stores the RHS data.
#ifdef HAVE_C_VARARRAYS
  float bands[4*M], lastCol[M];
#else
  float *bands   = new float[4*M];
  float *lastCol = new float[  M];
#endif
 
  // Fill up bands
  for (int i=0; i<M; i++) {
    get_NUBasis_funcs_si (basis, i, &(bands[4*i])); 
    bands[4*(i)+3] = data[i*datastride];
  }
    
  // Now solve:
  // First and last rows are different
  bands[4*(0)+2] /= bands[4*(0)+1];
  bands[4*(0)+0] /= bands[4*(0)+1];
  bands[4*(0)+3] /= bands[4*(0)+1];
  bands[4*(0)+1]  = 1.0;
  bands[4*(M-1)+1] -= bands[4*(M-1)+2]*bands[4*(0)+0];
  bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(0)+3];
  bands[4*(M-1)+2]  = -bands[4*(M-1)+2]*bands[4*(0)+2];
  lastCol[0] = bands[4*(0)+0];
  
  for (int row=1; row < (M-1); row++) {
    bands[4*(row)+1] -= bands[4*(row)+0] * bands[4*(row-1)+2];
    bands[4*(row)+3] -= bands[4*(row)+0] * bands[4*(row-1)+3];
    lastCol[row]   = -bands[4*(row)+0] * lastCol[row-1];
    bands[4*(row)+0] = 0.0;
    bands[4*(row)+2] /= bands[4*(row)+1];
    bands[4*(row)+3] /= bands[4*(row)+1];
    lastCol[row]  /= bands[4*(row)+1];
    bands[4*(row)+1]  = 1.0;
    if (row < (M-2)) {
      bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(row)+3];
      bands[4*(M-1)+1] -= bands[4*(M-1)+2]*lastCol[row];
      bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(row)+2];
    }
  }
  
  // Now do last row
  // The [2] element and [0] element are now on top of each other 
  bands[4*(M-1)+0] += bands[4*(M-1)+2];
  bands[4*(M-1)+1] -= bands[4*(M-1)+0] * (bands[4*(M-2)+2]+lastCol[M-2]);
  bands[4*(M-1)+3] -= bands[4*(M-1)+0] *  bands[4*(M-2)+3];
  bands[4*(M-1)+3] /= bands[4*(M-1)+1];
  p[pstride*M] = bands[4*(M-1)+3];
  for (int row=M-2; row>=0; row--) 
    p[pstride*(row+1)] = bands[4*(row)+3] - 
      bands[4*(row)+2]*p[pstride*(row+2)] - lastCol[row]*p[pstride*M];
  
  p[pstride*  0  ] = p[pstride*M];
  p[pstride*(M+1)] = p[pstride*1];
  p[pstride*(M+2)] = p[pstride*2];
#ifndef HAVE_C_VARARRAYS
  delete[] bands;
  delete[] lastCol;
#endif
}

References get_NUBasis_funcs_si(), and p.