multi_bspline_create.cpp File Reference

#include "multi_bspline_create.h"
#include <stdlib.h>
#include <stdio.h>

Go to the source code of this file.

Macros
#define	__USE_XOPEN2K
#define	_XOPEN_SOURCE   600

Functions
multi_UBspline_1d_c *	create_multi_UBspline_1d_c (Ugrid x_grid, BCtype_c xBC, int num_splines)
multi_UBspline_1d_d *	create_multi_UBspline_1d_d (Ugrid x_grid, BCtype_d xBC, int num_splines)
multi_UBspline_1d_s *	create_multi_UBspline_1d_s (Ugrid x_grid, BCtype_s xBC, int num_splines)
multi_UBspline_1d_z *	create_multi_UBspline_1d_z (Ugrid x_grid, BCtype_z 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_2d_d *	create_multi_UBspline_2d_d (Ugrid x_grid, Ugrid y_grid, BCtype_d xBC, BCtype_d yBC, 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_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_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)
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)
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)
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	destroy_multi_UBspline (Bspline *spline)
void	find_coefs_1d_d (Ugrid grid, BCtype_d bc, double *data, intptr_t dstride, double *coefs, intptr_t cstride)
void	find_coefs_1d_s (Ugrid grid, BCtype_s bc, float *data, intptr_t dstride, float *coefs, intptr_t cstride)
void	init_sse_data ()
int	posix_memalign (void **memptr, size_t alignment, size_t size)
void	set_multi_UBspline_1d_c (multi_UBspline_1d_c *spline, int num, complex_float *data)
void	set_multi_UBspline_1d_d (multi_UBspline_1d_d *spline, int num, double *data)
void	set_multi_UBspline_1d_d_BC (multi_UBspline_1d_d *spline, int num, double *data, BCtype_d xBC)
void	set_multi_UBspline_1d_s (multi_UBspline_1d_s *spline, int num, float *data)
void	set_multi_UBspline_1d_z (multi_UBspline_1d_z *spline, int num, complex_double *data)
void	set_multi_UBspline_1d_z_BC (multi_UBspline_1d_z *spline, int num, complex_double *data, BCtype_z xBC)
void	set_multi_UBspline_2d_c (multi_UBspline_2d_c *spline, int num, complex_float *data)
void	set_multi_UBspline_2d_d (multi_UBspline_2d_d *spline, int num, double *data)
void	set_multi_UBspline_2d_s (multi_UBspline_2d_s *spline, int num, float *data)
void	set_multi_UBspline_2d_z (multi_UBspline_2d_z *spline, int num, complex_double *data)
void	set_multi_UBspline_3d_c (multi_UBspline_3d_c *spline, int num, complex_float *data)
void	set_multi_UBspline_3d_d (multi_UBspline_3d_d *spline, int num, double *data)
void	set_multi_UBspline_3d_s (multi_UBspline_3d_s *spline, int num, float *data)
void	set_multi_UBspline_3d_z (multi_UBspline_3d_z *spline, int num, complex_double *data)
void	solve_antiperiodic_interp_1d_s (float bands[], float coefs[], int M, int cstride)
void	solve_deriv_interp_1d_d (double bands[], double coefs[], int M, int cstride)
void	solve_deriv_interp_1d_s (float bands[], float coefs[], int M, int cstride)
void	solve_periodic_interp_1d_d (double bands[], double coefs[], int M, int cstride)
void	solve_periodic_interp_1d_s (float bands[], float coefs[], int M, int cstride)

Macro Definition Documentation

__USE_XOPEN2K

#define __USE_XOPEN2K

Definition at line 26 of file multi_bspline_create.cpp.

_XOPEN_SOURCE

#define _XOPEN_SOURCE   600

Definition at line 23 of file multi_bspline_create.cpp.

Function Documentation

create_multi_UBspline_1d_c()

multi_UBspline_1d_c * create_multi_UBspline_1d_c	(	Ugrid	x_grid,
		BCtype_c	xBC,
		int	num_splines )

Definition at line 388 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_1d_c* restrict spline = static_cast<multi_UBspline_1d_c*>(malloc(sizeof(multi_UBspline_1d_c)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_c.\n");
    abort();
  }
  spline->spcode = MULTI_U1D;
  spline->tcode  = SINGLE_COMPLEX;
  spline->xBC = xBC; 
  spline->num_splines = num_splines;
  // Setup internal variables
  int M = x_grid.num;
  int N;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num);
    N = M+3;
  }
  else {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
    N = M+2;
  }
 
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
  spline->x_stride = num_splines;
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (complex_float*)malloc (2*sizeof(float)*N*num_splines);
#else
  posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(float)*N*num_splines);
#endif
#ifdef HAVE_SSE
  init_sse_data();    
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_c.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_c::lCode, MULTI_U1D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_COMPLEX, and Ugrid::start.

create_multi_UBspline_1d_d()

multi_UBspline_1d_d * create_multi_UBspline_1d_d	(	Ugrid	x_grid,
		BCtype_d	xBC,
		int	num_splines )

Definition at line 766 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_1d_d* restrict spline = static_cast<multi_UBspline_1d_d*>(malloc(sizeof(multi_UBspline_1d_d)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_d.\n");
    abort();
  }
  spline->spcode = MULTI_U1D;
  spline->tcode  = DOUBLE_REAL;
  spline->xBC = xBC; 
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int Nx;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num);
    Nx = Mx+3;
  }
  else {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
    Nx = Mx+2;
  }
 
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  int N = num_splines;
#ifdef HAVE_SSE2
  // We must pad to keep data aligned for SSE operations
  if (N & 1)
    N++;
#endif
  spline->x_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (double*)malloc (sizeof(double)*Nx*N);
 
#else
  posix_memalign ((void**)&spline->coefs, 64, sizeof(double)*Nx*N);
#endif
#ifdef HAVE_SSE2
  init_sse_data();
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_d.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_REAL, Ugrid::end, init_sse_data(), BCtype_d::lCode, MULTI_U1D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

create_multi_UBspline_1d_s()

multi_UBspline_1d_s * create_multi_UBspline_1d_s	(	Ugrid	x_grid,
		BCtype_s	xBC,
		int	num_splines )

Definition at line 87 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_1d_s* restrict spline = static_cast<multi_UBspline_1d_s*>(malloc(sizeof(multi_UBspline_1d_s)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_s.\n");
    abort();
  }
  spline->spcode = MULTI_U1D;
  spline->tcode  = SINGLE_REAL;
  spline->xBC = xBC; spline->x_grid = x_grid;
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int Nx;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num);
    Nx = Mx+3;
  }
  else {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
    Nx = Mx+2;
  }
 
  int N = num_splines;
#ifdef HAVE_SSE
  if (N % 4) 
    N += 4 - (N % 4);
#endif 
 
  spline->x_stride = N;
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (float*)malloc (sizeof(float)*Nx*N);
#else
  posix_memalign ((void**)&spline->coefs, 64, (sizeof(float)*Nx*N));
#endif
#ifdef HAVE_SSE
  init_sse_data();    
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficient in create_multi_UBspline_1d_s.\n");
    abort();
  }
 
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_s::lCode, MULTI_U1D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_REAL, and Ugrid::start.

create_multi_UBspline_1d_z()

multi_UBspline_1d_z * create_multi_UBspline_1d_z	(	Ugrid	x_grid,
		BCtype_z	xBC,
		int	num_splines )

Definition at line 1090 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_1d_z* restrict spline = static_cast<multi_UBspline_1d_z*>(malloc(sizeof(multi_UBspline_1d_z)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_z.\n");
    abort();
  }
  spline->spcode = MULTI_U1D;
  spline->tcode  = DOUBLE_COMPLEX;
  spline->xBC = xBC; 
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int Nx;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num);
    Nx = Mx+3;
  }
  else {
    x_grid.delta     = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
    Nx = Mx+2;
  }
 
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
  spline->x_stride = num_splines;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (complex_double*)malloc (2*sizeof(double)*Nx*num_splines);
#else
  posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(double)*Nx*num_splines);
#endif
#ifdef HAVE_SSE2
  init_sse_data();   
#endif
   if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_z.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_COMPLEX, Ugrid::end, init_sse_data(), BCtype_z::lCode, MULTI_U1D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

create_multi_UBspline_2d_c()

multi_UBspline_2d_c * create_multi_UBspline_2d_c	(	Ugrid	x_grid,
		Ugrid	y_grid,
		BCtype_c	xBC,
		BCtype_c	yBC,
		int	num_splines )

Definition at line 455 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_2d_c* restrict spline = static_cast<multi_UBspline_2d_c*>(malloc(sizeof(multi_UBspline_2d_c)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_c.\n");
    abort();
  }
  spline->spcode = MULTI_U2D;
  spline->tcode  = SINGLE_COMPLEX;
  spline->xBC = xBC; 
  spline->yBC = yBC;
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int My = y_grid.num;
  int Nx, Ny;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  int N = num_splines;
#ifdef HAVE_SSE
  if (N % 2)
    N++;
#endif
 
  spline->x_stride = Ny*N;
  spline->y_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (complex_float*)malloc (2*sizeof(float)*Nx*Ny*N);
  spline->lapl2 = (complex_float*)malloc (4*sizeof(float)*N);
#else
  posix_memalign ((void**)&spline->coefs, 64, 
          2*sizeof(float)*Nx*Ny*N);
  posix_memalign ((void**)&spline->lapl2, 64,
          4*sizeof(float)*N);
#endif
#ifdef HAVE_SSE
  init_sse_data();
#endif
  if (!spline->coefs || !spline->lapl2) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_c.\n");
    abort();
  }
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_c::lCode, MULTI_U2D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_COMPLEX, and Ugrid::start.

create_multi_UBspline_2d_d()

multi_UBspline_2d_d * create_multi_UBspline_2d_d	(	Ugrid	x_grid,
		Ugrid	y_grid,
		BCtype_d	xBC,
		BCtype_d	yBC,
		int	num_splines )

Definition at line 839 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_2d_d* restrict spline = static_cast<multi_UBspline_2d_d*>(malloc(sizeof(multi_UBspline_2d_d)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_d.\n");
    abort();
  }
  spline->spcode = MULTI_U2D;
  spline->tcode  = DOUBLE_REAL;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int My = y_grid.num;
  int Nx, Ny;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  int N = num_splines;
 
#ifdef HAVE_SSE2
  // We must pad to keep data align for SSE operations
  if (num_splines & 1)
    N++;
#endif
  spline->x_stride = Ny*N;
  spline->y_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny*N);
#else
  posix_memalign ((void**)&spline->coefs, 64, (sizeof(double)*Nx*Ny*N));
#endif
#ifdef HAVE_SSE2
  init_sse_data();
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_d.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_REAL, Ugrid::end, init_sse_data(), BCtype_d::lCode, MULTI_U2D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

create_multi_UBspline_2d_s()

multi_UBspline_2d_s * create_multi_UBspline_2d_s	(	Ugrid	x_grid,
		Ugrid	y_grid,
		BCtype_s	xBC,
		BCtype_s	yBC,
		int	num_splines )

Definition at line 151 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_2d_s* restrict spline = static_cast<multi_UBspline_2d_s*>(malloc(sizeof(multi_UBspline_2d_s)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_s.\n");
    abort();
  }
  spline->spcode = MULTI_U2D;
  spline->tcode  = SINGLE_REAL;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->num_splines = num_splines;
  // Setup internal variables
  int Mx = x_grid.num;
  int My = y_grid.num;
  int Nx, Ny;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
 
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  int N = num_splines;
#ifdef HAVE_SSE
  if (N % 4) 
    N += 4 - (N % 4);
#endif
 
 
  spline->x_stride = Ny*N;
  spline->y_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny*N);
#else
  posix_memalign ((void**)&spline->coefs, 64, 
          sizeof(float)*Nx*Ny*N);
#endif
#ifdef HAVE_SSE
  init_sse_data();
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_s.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_s::lCode, MULTI_U2D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_REAL, and Ugrid::start.

create_multi_UBspline_2d_z()

multi_UBspline_2d_z * create_multi_UBspline_2d_z	(	Ugrid	x_grid,
		Ugrid	y_grid,
		BCtype_z	xBC,
		BCtype_z	yBC,
		int	num_splines )

Definition at line 1199 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_2d_z* restrict spline = static_cast<multi_UBspline_2d_z*>(malloc(sizeof(multi_UBspline_2d_z)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_z.\n");
    abort();
  }
  spline->spcode = MULTI_U2D;
  spline->tcode  = DOUBLE_COMPLEX;
  spline->xBC = xBC; 
  spline->yBC = yBC;
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;
  int My = y_grid.num;
  int Nx, Ny;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
  spline->x_stride = Ny*num_splines;
  spline->y_stride = num_splines;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (complex_double*)malloc (2*sizeof(double)*Nx*Ny*num_splines);
  spline->lapl2 = (complex_double*)malloc (4*sizeof(double)*num_splines);
#else
  posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(double)*Nx*Ny*num_splines);
  posix_memalign ((void**)&spline->lapl2, 64, 4*sizeof(double)*num_splines);
#endif
#ifdef HAVE_SSE2
  init_sse_data();
#endif
  if (!spline->coefs || !spline->lapl2) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_z.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_COMPLEX, Ugrid::end, init_sse_data(), BCtype_z::lCode, MULTI_U2D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

create_multi_UBspline_3d_c()

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 )

Definition at line 576 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_3d_c* restrict spline = static_cast<multi_UBspline_3d_c*>(malloc(sizeof(multi_UBspline_3d_c)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_c.\n");
    abort();
  }
  spline->spcode = MULTI_U3D;
  spline->tcode  = SINGLE_COMPLEX;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->zBC = zBC; 
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;  int My = y_grid.num; int Mz = z_grid.num;
  int Nx, Ny, Nz;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                           
    Nz = Mz+2;
  z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
  z_grid.delta_inv = 1.0/z_grid.delta;
  spline->z_grid   = z_grid;
 
  int N = spline->num_splines;
#ifdef HAVE_SSE
  if (N % 2)
    N++;
#endif
 
  spline->x_stride = Ny*Nz*N;
  spline->y_stride = Nz*N;
  spline->z_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs = (complex_float*)malloc (2*sizeof(float)*Nx*Ny*Nz*N);
  spline->lapl3 = (complex_float*)malloc (6*sizeof(float)*N);
#else
  posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(float)*Nx*Ny*Nz*N);
  posix_memalign ((void**)&spline->lapl3, 64, 6*sizeof(float)*N);  
#endif
#ifdef HAVE_SSE
  init_sse_data();
#endif
  if (!spline->coefs || !spline->lapl3) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_c.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_c::lCode, MULTI_U3D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_COMPLEX, and Ugrid::start.

create_multi_UBspline_3d_d()

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 )

Definition at line 941 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_3d_d* restrict spline;
#ifdef HAVE_POSIX_MEMALIGN
  posix_memalign ((void**)&spline, 64, (size_t)sizeof(multi_UBspline_3d_d));
#else
  spline = static_cast<multi_UBspline_3d_d*>(malloc(sizeof(multi_UBspline_3d_d)));
#endif
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_d.\n");
    abort();
  }
  spline->spcode = MULTI_U3D;
  spline->tcode  = DOUBLE_REAL;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->zBC = zBC; 
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;  int My = y_grid.num; int Mz = z_grid.num;
  int Nx, Ny, Nz;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                           
    Nz = Mz+2;
  z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
  z_grid.delta_inv = 1.0/z_grid.delta;
  spline->z_grid   = z_grid;
 
 
  int N = num_splines;
#if defined HAVE_SSE2 || defined HAVE_VSX
  // We must pad to keep data align for SSE operations
  if (N & 1)
    N++;
#endif
  
  spline->x_stride = Ny*Nz*N;
  spline->y_stride = Nz*N;
  spline->z_stride = N;
  
#ifdef HAVE_POSIX_MEMALIGN
  posix_memalign ((void**)&spline->coefs, 64, ((size_t)sizeof(double)*Nx*Ny*Nz*N));
#else
  spline->coefs      = new double[Nx*Ny*Nz*N];
#endif
#ifdef HAVE_SSE2
  init_sse_data();
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_d.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_REAL, Ugrid::end, init_sse_data(), BCtype_d::lCode, MULTI_U3D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

create_multi_UBspline_3d_s()

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 )

Definition at line 251 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_3d_s* restrict spline = static_cast<multi_UBspline_3d_s*>(malloc(sizeof(multi_UBspline_3d_s)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_s.\n");
    abort();
  }
  spline->spcode = MULTI_U3D;
  spline->tcode  = SINGLE_REAL;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->zBC = zBC; 
  spline->num_splines = num_splines;
  // Setup internal variables
  int Mx = x_grid.num;  int My = y_grid.num; int Mz = z_grid.num;
  int Nx, Ny, Nz;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                           
    Nz = Mz+2;
  z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
  z_grid.delta_inv = 1.0/z_grid.delta;
  spline->z_grid   = z_grid;
 
  int N = num_splines;
#ifdef HAVE_SSE
  if (N % 4) 
    N += 4 - (N % 4);
#endif
 
  spline->x_stride      = Ny*Nz*N;
  spline->y_stride      = Nz*N;
  spline->z_stride      = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs      = new float[Nx*Ny*Nz*N];
#else
  posix_memalign ((void**)&spline->coefs, 64, 
          (sizeof(float)*Nx*Ny*Nz*N));
#endif
#ifdef HAVE_SSE
  init_sse_data();
#endif
  if (!spline->coefs) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_s.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, Ugrid::end, init_sse_data(), BCtype_s::lCode, MULTI_U3D, Ugrid::num, PERIODIC, posix_memalign(), restrict, SINGLE_REAL, and Ugrid::start.

create_multi_UBspline_3d_z()

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 )

Definition at line 1318 of file multi_bspline_create.cpp.

{
  // Create new spline
  multi_UBspline_3d_z* restrict spline = static_cast<multi_UBspline_3d_z*>(malloc(sizeof(multi_UBspline_3d_z)));
  if (!spline) {
    fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_z.\n");
    abort();
  }
  spline->spcode = MULTI_U3D;
  spline->tcode  = DOUBLE_COMPLEX;
  spline->xBC = xBC; 
  spline->yBC = yBC; 
  spline->zBC = zBC;
  spline->num_splines = num_splines;
 
  // Setup internal variables
  int Mx = x_grid.num;  int My = y_grid.num; int Mz = z_grid.num;
  int Nx, Ny, Nz;
 
  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
  x_grid.delta_inv = 1.0/x_grid.delta;
  spline->x_grid   = x_grid;
 
  if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
  y_grid.delta_inv = 1.0/y_grid.delta;
  spline->y_grid   = y_grid;
 
  if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                           
    Nz = Mz+2;
  z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
  z_grid.delta_inv = 1.0/z_grid.delta;
  spline->z_grid   = z_grid;
 
  int N = num_splines;
#ifdef HAVE_SSE2
  if (N & 3)
    N += 4-(N & 3);
#endif
 
  spline->x_stride = (intptr_t)Ny*(intptr_t)Nz*N;
  spline->y_stride = Nz*N;
  spline->z_stride = N;
 
#ifndef HAVE_POSIX_MEMALIGN
  spline->coefs      = new complex_double[Nx*Ny*Nz*N];
  spline->lapl3 = new complex_double[3*N];
#else
  posix_memalign ((void**)&spline->coefs, 64, (size_t)2*sizeof(double)*Nx*Ny*Nz*N);
  posix_memalign ((void**)&spline->lapl3, 64, 6*sizeof(double)*N);
#endif
 
#ifdef HAVE_SSE2
  init_sse_data();
#endif
  if (!spline->coefs || !spline->lapl3) {
    fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_z.\n");
    abort();
  }
 
  return spline;
}

References ANTIPERIODIC, Ugrid::delta, Ugrid::delta_inv, DOUBLE_COMPLEX, Ugrid::end, init_sse_data(), BCtype_z::lCode, MULTI_U3D, Ugrid::num, PERIODIC, posix_memalign(), restrict, and Ugrid::start.

destroy_multi_UBspline()

void destroy_multi_UBspline

(

Bspline *

spline

)

Definition at line 1481 of file multi_bspline_create.cpp.

{
  free(spline->coefs);
  free(spline);
}

References Bspline::coefs.

find_coefs_1d_d()

void find_coefs_1d_d	(	Ugrid	grid,
		BCtype_d	bc,
		double *	data,
		intptr_t	dstride,
		double *	coefs,
		intptr_t	cstride )

Definition at line 1118 of file bspline_create.cpp.

{
  int M = grid.num;
  double basis[4] = {1.0/6.0, 2.0/3.0, 1.0/6.0, 0.0};
  if (bc.lCode == PERIODIC || bc.lCode == ANTIPERIODIC) {
#ifdef HAVE_C_VARARRAYS
    double bands[M*4];
#else
    double *bands = new double[4*M];
#endif
    for (int i=0; i<M; i++) {
      bands[4*i+0] = basis[0];
      bands[4*i+1] = basis[1];
      bands[4*i+2] = basis[2];
      bands[4*i+3] = data[i*dstride];
    }
    if (bc.lCode == ANTIPERIODIC) 
      solve_antiperiodic_interp_1d_d (bands, coefs, M, cstride);
    else
      solve_periodic_interp_1d_d (bands, coefs, M, cstride);
 
 
#ifndef HAVE_C_VARARRAYS
    delete[] bands;
#endif
  }
  else {
    // Setup boundary conditions
    double abcd_left[4];
    double abcd_right[4];
    for (int i = 0; i < 4; i++) {
        abcd_left[i] = 0.0;
        abcd_right[i] = 0.0;
    }
    // Left boundary
    if (bc.lCode == FLAT || bc.lCode == NATURAL)
      bc.lVal = 0.0;
    if (bc.lCode == FLAT || bc.lCode == DERIV1) {
      abcd_left[0] = -0.5 * grid.delta_inv;
      abcd_left[1] =  0.0 * grid.delta_inv; 
      abcd_left[2] =  0.5 * grid.delta_inv;
      abcd_left[3] =  bc.lVal;
    }
    if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
      abcd_left[0] = 1.0 * grid.delta_inv * grid.delta_inv;
      abcd_left[1] =-2.0 * grid.delta_inv * grid.delta_inv;
      abcd_left[2] = 1.0 * grid.delta_inv * grid.delta_inv;
      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) {
      abcd_right[0] = -0.5 * grid.delta_inv;
      abcd_right[1] =  0.0 * grid.delta_inv; 
      abcd_right[2] =  0.5 * grid.delta_inv;
      abcd_right[3] =  bc.rVal;
    }
    if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
      abcd_right[0] = 1.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[1] =-2.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[2] = 1.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[3] = bc.rVal;
    }
#ifdef HAVE_C_VARARRAYS
    double bands[(M+2)*4];
#else
    double *bands = new double[(M+2)*4];
#endif
    for (int i=0; i<4; i++) {
      bands[4*( 0 )+i] = abcd_left[i];
      bands[4*(M+1)+i] = abcd_right[i];
    }
    for (int i=0; i<M; i++) {
      for (int j=0; j<3; j++)
    bands[4*(i+1)+j] = basis[j];
      bands[4*(i+1)+3] = data[i*dstride];
    }   
    // Now, solve for coefficients
    solve_deriv_interp_1d_d (bands, coefs, M, cstride);
#ifndef HAVE_C_VARARRAYS
    delete[] bands;
#endif
  }
}

References ANTIPERIODIC, Ugrid::delta_inv, DERIV1, DERIV2, FLAT, BCtype_d::lCode, BCtype_d::lVal, NATURAL, Ugrid::num, PERIODIC, BCtype_d::rCode, BCtype_d::rVal, solve_antiperiodic_interp_1d_d(), solve_deriv_interp_1d_d(), and solve_periodic_interp_1d_d().

find_coefs_1d_s()

void find_coefs_1d_s	(	Ugrid	grid,
		BCtype_s	bc,
		float *	data,
		intptr_t	dstride,
		float *	coefs,
		intptr_t	cstride )

Definition at line 240 of file bspline_create.cpp.

{
  int M = grid.num;
  float basis[4] = {1.0/6.0, 2.0/3.0, 1.0/6.0, 0.0};
  if (bc.lCode == PERIODIC || bc.lCode == ANTIPERIODIC) {
 
    float *bands = new float[4*M];
 
    for (int i=0; i<M; i++) {
      bands[4*i+0] = basis[0];
      bands[4*i+1] = basis[1];
      bands[4*i+2] = basis[2];
      bands[4*i+3] = data[i*dstride];
    }
    if (bc.lCode == PERIODIC)
      solve_periodic_interp_1d_s (bands, coefs, M, cstride);
    else
      solve_antiperiodic_interp_1d_s (bands, coefs, M, cstride);
 
    delete[] bands;
  }
  else {
    // Setup boundary conditions
    float abcd_left[4];
    float abcd_right[4];
    for (int i = 0; i < 4; i++) {
        abcd_left[i] = 0.0;
        abcd_right[i] = 0.0;
    }
 
    // Left boundary
    if (bc.lCode == FLAT || bc.lCode == NATURAL)
      bc.lVal = 0.0;
    if (bc.lCode == FLAT || bc.lCode == DERIV1) {
      abcd_left[0] = -0.5 * grid.delta_inv;
      abcd_left[1] =  0.0 * grid.delta_inv; 
      abcd_left[2] =  0.5 * grid.delta_inv;
      abcd_left[3] =  bc.lVal;
    }
    if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
      abcd_left[0] = 1.0 * grid.delta_inv * grid.delta_inv;
      abcd_left[1] =-2.0 * grid.delta_inv * grid.delta_inv;
      abcd_left[2] = 1.0 * grid.delta_inv * grid.delta_inv;
      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) {
      abcd_right[0] = -0.5 * grid.delta_inv;
      abcd_right[1] =  0.0 * grid.delta_inv; 
      abcd_right[2] =  0.5 * grid.delta_inv;
      abcd_right[3] =  bc.rVal;
    }
    if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
      abcd_right[0] = 1.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[1] =-2.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[2] = 1.0 *grid.delta_inv * grid.delta_inv;
      abcd_right[3] = bc.rVal;
    }
 
    float *bands = new float[4*(M+2)];
    for (int i=0; i<4; i++) {
      bands[4*( 0 )+i]   = abcd_left[i];
      bands[4*(M+1)+i] = abcd_right[i];
    }
    for (int i=0; i<M; i++) {
      for (int j=0; j<3; j++)
    bands[4*(i+1)+j] = basis[j];
      bands[4*(i+1)+3] = data[i*dstride];
    }   
    // Now, solve for coefficients
    solve_deriv_interp_1d_s (bands, coefs, M, cstride);
    delete[] bands;
  }
}

References ANTIPERIODIC, Ugrid::delta_inv, DERIV1, DERIV2, FLAT, BCtype_s::lCode, BCtype_s::lVal, NATURAL, Ugrid::num, PERIODIC, BCtype_s::rCode, BCtype_s::rVal, solve_antiperiodic_interp_1d_s(), solve_deriv_interp_1d_s(), and solve_periodic_interp_1d_s().

init_sse_data()

void init_sse_data

(

)

Definition at line 62 of file bspline_data.cpp.

{
#ifdef HAVE_SSE
  if (A_s == 0) {
    posix_memalign ((void**)&A_s, 16, (sizeof(__m128)*12));
    A_s[0]  = _mm_setr_ps ( 1.0/6.0, -3.0/6.0,  3.0/6.0, -1.0/6.0 );
    A_s[0]  = _mm_setr_ps ( 1.0/6.0, -3.0/6.0,  3.0/6.0, -1.0/6.0 );      
    A_s[1]  = _mm_setr_ps ( 4.0/6.0,  0.0/6.0, -6.0/6.0,  3.0/6.0 );      
    A_s[2]  = _mm_setr_ps ( 1.0/6.0,  3.0/6.0,  3.0/6.0, -3.0/6.0 );      
    A_s[3]  = _mm_setr_ps ( 0.0/6.0,  0.0/6.0,  0.0/6.0,  1.0/6.0 );      
    A_s[4]  = _mm_setr_ps ( -0.5,  1.0, -0.5, 0.0  );         
    A_s[5]  = _mm_setr_ps (  0.0, -2.0,  1.5, 0.0  );         
    A_s[6]  = _mm_setr_ps (  0.5,  1.0, -1.5, 0.0  );         
    A_s[7]  = _mm_setr_ps (  0.0,  0.0,  0.5, 0.0  );         
    A_s[8]  = _mm_setr_ps (  1.0, -1.0,  0.0, 0.0  );         
    A_s[9]  = _mm_setr_ps ( -2.0,  3.0,  0.0, 0.0  );         
    A_s[10] = _mm_setr_ps (  1.0, -3.0,  0.0, 0.0  );         
    A_s[11] = _mm_setr_ps (  0.0,  1.0,  0.0, 0.0  );                  
  }
                 
#endif
#ifdef HAVE_SSE2
  if (A_d == 0) {
    posix_memalign ((void**)&A_d, 16, (sizeof(__m128d)*24));
    A_d[ 0] = _mm_setr_pd (  3.0/6.0, -1.0/6.0 );      
    A_d[ 1] = _mm_setr_pd (  1.0/6.0, -3.0/6.0 );      
    A_d[ 2] = _mm_setr_pd ( -6.0/6.0,  3.0/6.0 );      
    A_d[ 3] = _mm_setr_pd (  4.0/6.0,  0.0/6.0 );      
    A_d[ 4] = _mm_setr_pd (  3.0/6.0, -3.0/6.0 );      
    A_d[ 5] = _mm_setr_pd (  1.0/6.0,  3.0/6.0 );      
    A_d[ 6] = _mm_setr_pd (  0.0/6.0,  1.0/6.0 );      
    A_d[ 7] = _mm_setr_pd (  0.0/6.0,  0.0/6.0 );      
    A_d[ 8] = _mm_setr_pd ( -0.5,  0.0 );          
    A_d[ 9] = _mm_setr_pd ( -0.5,  1.0 );          
    A_d[10] = _mm_setr_pd (  1.5,  0.0 );          
    A_d[11] = _mm_setr_pd (  0.0, -2.0 );          
    A_d[12] = _mm_setr_pd ( -1.5,  0.0 );          
    A_d[13] = _mm_setr_pd (  0.5,  1.0 );          
    A_d[14] = _mm_setr_pd (  0.5,  0.0 );          
    A_d[15] = _mm_setr_pd (  0.0,  0.0 );          
    A_d[16] = _mm_setr_pd (  0.0,  0.0 );          
    A_d[17] = _mm_setr_pd (  1.0, -1.0 );          
    A_d[18] = _mm_setr_pd (  0.0,  0.0 );          
    A_d[19] = _mm_setr_pd ( -2.0,  3.0 );          
    A_d[20] = _mm_setr_pd (  0.0,  0.0 );          
    A_d[21] = _mm_setr_pd (  1.0, -3.0 );          
    A_d[22] = _mm_setr_pd (  0.0,  0.0 );          
    A_d[23] = _mm_setr_pd (  0.0,  1.0 );   
  }                
#endif
}

posix_memalign()

int posix_memalign	(	void **	memptr,
		size_t	alignment,
		size_t	size )

set_multi_UBspline_1d_c()

void set_multi_UBspline_1d_c	(	multi_UBspline_1d_c *	spline,
		int	num,
		complex_float *	data )

Definition at line 433 of file multi_bspline_create.cpp.

{
  complex_float *coefs = spline->coefs + num;
 
  BCtype_s xBC_r, xBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
 
  int xs = spline->x_stride;
  // Real part
  find_coefs_1d_s (spline->x_grid, xBC_r, 
           (float*)data, (intptr_t)2, (float*)coefs, (intptr_t)2*xs);
  // Imaginarty part
  find_coefs_1d_s (spline->x_grid, xBC_i, 
           ((float*)data)+1, (intptr_t)2, ((float*)coefs+1), (intptr_t)2*xs);
}

References multi_UBspline_1d_c::coefs, find_coefs_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, BCtype_c::rVal_r, multi_UBspline_1d_c::x_grid, multi_UBspline_1d_c::x_stride, and multi_UBspline_1d_c::xBC.

set_multi_UBspline_1d_d()

void set_multi_UBspline_1d_d	(	multi_UBspline_1d_d *	spline,
		int	num,
		double *	data )

Definition at line 821 of file multi_bspline_create.cpp.

{
  double *coefs = spline->coefs + num;
  int xs = spline->x_stride;
  find_coefs_1d_d (spline->x_grid, spline->xBC, data, 1, coefs, xs);
}

References multi_UBspline_1d_d::coefs, find_coefs_1d_d(), multi_UBspline_1d_d::x_grid, multi_UBspline_1d_d::x_stride, and multi_UBspline_1d_d::xBC.

set_multi_UBspline_1d_d_BC()

void set_multi_UBspline_1d_d_BC	(	multi_UBspline_1d_d *	spline,
		int	num,
		double *	data,
		BCtype_d	xBC )

Definition at line 829 of file multi_bspline_create.cpp.

{
  double *coefs = spline->coefs + num;
  int xs = spline->x_stride;
  find_coefs_1d_d (spline->x_grid, xBC, data, 1, coefs, xs);
}

References multi_UBspline_1d_d::coefs, find_coefs_1d_d(), multi_UBspline_1d_d::x_grid, and multi_UBspline_1d_d::x_stride.

set_multi_UBspline_1d_s()

void set_multi_UBspline_1d_s	(	multi_UBspline_1d_s *	spline,
		int	num,
		float *	data )

Definition at line 140 of file multi_bspline_create.cpp.

{
  float *coefs = spline->coefs + num;
  int xs = spline->x_stride;
  find_coefs_1d_s (spline->x_grid, spline->xBC, data, 1, 
           coefs, xs);
}

References multi_UBspline_1d_s::coefs, find_coefs_1d_s(), multi_UBspline_1d_s::x_grid, multi_UBspline_1d_s::x_stride, and multi_UBspline_1d_s::xBC.

set_multi_UBspline_1d_z()

void set_multi_UBspline_1d_z	(	multi_UBspline_1d_z *	spline,
		int	num,
		complex_double *	data )

Definition at line 1137 of file multi_bspline_create.cpp.

{
//  int Mx = spline->x_grid.num;
// Set but not used
//  int Nx;
 
  complex_double *coefs = spline->coefs + num;
 
//  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
//    Nx = Mx+3;
//  else
//    Nx = Mx+2;
 
  BCtype_d xBC_r, xBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
  int xs = spline->x_stride;
  // Real part
  find_coefs_1d_d (spline->x_grid, xBC_r, 
           (double*)data,      (intptr_t)2, 
           ((double*)coefs),   (intptr_t)2*xs);
  // Imaginary part
  find_coefs_1d_d (spline->x_grid, xBC_i, 
           ((double*)data)+1,  (intptr_t)2, 
           ((double*)coefs)+1, (intptr_t)2*xs);
 
}

References multi_UBspline_1d_z::coefs, find_coefs_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, BCtype_z::rVal_r, multi_UBspline_1d_z::x_grid, multi_UBspline_1d_z::x_stride, and multi_UBspline_1d_z::xBC.

set_multi_UBspline_1d_z_BC()

void set_multi_UBspline_1d_z_BC	(	multi_UBspline_1d_z *	spline,
		int	num,
		complex_double *	data,
		BCtype_z	xBC )

Definition at line 1168 of file multi_bspline_create.cpp.

{
//  int Mx = spline->x_grid.num;
//  int Nx;
 
  complex_double *coefs = spline->coefs + num;
 
//  if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
//    Nx = Mx+3;
//  else
//    Nx = Mx+2;
 
  BCtype_d xBC_r, xBC_i;
  xBC_r.lCode = xBC.lCode;  xBC_r.rCode = xBC.rCode;
  xBC_r.lVal  = xBC.lVal_r; xBC_r.rVal  = xBC.rVal_r;
  xBC_i.lCode = xBC.lCode;  xBC_i.rCode = xBC.rCode;
  xBC_i.lVal  = xBC.lVal_i; xBC_i.rVal  = xBC.rVal_i;
  int xs = spline->x_stride;
  // Real part
  find_coefs_1d_d (spline->x_grid, xBC_r, 
           (double*)data,    (intptr_t)2, 
           ((double*)coefs), (intptr_t)2*xs);
  // Imaginary part
  find_coefs_1d_d (spline->x_grid, xBC_i, 
           ((double*)data)+1,  (intptr_t)2, 
           ((double*)coefs)+1, (intptr_t)2*xs);
}

References multi_UBspline_1d_z::coefs, find_coefs_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, BCtype_z::rVal_r, multi_UBspline_1d_z::x_grid, and multi_UBspline_1d_z::x_stride.

set_multi_UBspline_2d_c()

void set_multi_UBspline_2d_c	(	multi_UBspline_2d_c *	spline,
		int	num,
		complex_float *	data )

Definition at line 521 of file multi_bspline_create.cpp.

{
  // Setup internal variables
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Nx, Ny;
 
  complex_float* coefs = spline->coefs + num;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
 
  BCtype_s xBC_r, xBC_i, yBC_r, yBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
  yBC_r.lCode = spline->yBC.lCode;  yBC_r.rCode = spline->yBC.rCode;
  yBC_r.lVal  = spline->yBC.lVal_r; yBC_r.rVal  = spline->yBC.rVal_r;
  yBC_i.lCode = spline->yBC.lCode;  yBC_i.rCode = spline->yBC.rCode;
  yBC_i.lVal  = spline->yBC.lVal_i; yBC_i.rVal  = spline->yBC.rVal_i;
 
  int ys = spline->y_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    intptr_t doffset = (2*iy);
    intptr_t coffset = (2*iy)*ys;
    // Real part
    find_coefs_1d_s (spline->x_grid, xBC_r, ((float*)data)+doffset, (intptr_t)2*My,
             (float*)coefs+coffset, (intptr_t)2*Ny*ys);
    // Imag part
    find_coefs_1d_s (spline->x_grid, xBC_i, ((float*)data)+doffset+1, (intptr_t)2*My,
             ((float*)coefs)+coffset+1, (intptr_t)2*Ny*ys);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    intptr_t doffset = (2*ix*Ny)*ys;
    intptr_t coffset = (2*ix*Ny)*ys;
    // Real part
    find_coefs_1d_s (spline->y_grid, yBC_r, ((float*)coefs)+doffset, 
             (intptr_t)2*ys, ((float*)coefs)+coffset, (intptr_t)2*ys);
    // Imag part
    find_coefs_1d_s (spline->y_grid, yBC_i, ((float*)coefs)+doffset+1, 
             (intptr_t)2*ys, ((float*)coefs)+coffset+1, (intptr_t)2*ys);
  }  
}

References ANTIPERIODIC, multi_UBspline_2d_c::coefs, find_coefs_1d_s(), BCtype_s::lCode, BCtype_c::lCode, BCtype_s::lVal, BCtype_c::lVal_i, BCtype_c::lVal_r, Ugrid::num, PERIODIC, BCtype_s::rCode, BCtype_c::rCode, BCtype_s::rVal, BCtype_c::rVal_i, BCtype_c::rVal_r, multi_UBspline_2d_c::x_grid, multi_UBspline_2d_c::xBC, multi_UBspline_2d_c::y_grid, multi_UBspline_2d_c::y_stride, and multi_UBspline_2d_c::yBC.

set_multi_UBspline_2d_d()

void set_multi_UBspline_2d_d	(	multi_UBspline_2d_d *	spline,
		int	num,
		double *	data )

Definition at line 902 of file multi_bspline_create.cpp.

{
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Nx, Ny;
  double *coefs = spline->coefs + num;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
 
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
 
  int ys = spline->y_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    intptr_t doffset = iy;
    intptr_t coffset = iy*ys;
    find_coefs_1d_d (spline->x_grid, spline->xBC, 
             data+doffset,  (intptr_t)My,
             coefs+coffset, (intptr_t)Ny*ys);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    intptr_t doffset = ix*Ny*ys;
    intptr_t coffset = ix*Ny*ys;
    find_coefs_1d_d (spline->y_grid, spline->yBC, 
             coefs+doffset, (intptr_t)ys, 
             coefs+coffset, (intptr_t)ys);
  }
}

References ANTIPERIODIC, multi_UBspline_2d_d::coefs, find_coefs_1d_d(), BCtype_d::lCode, Ugrid::num, PERIODIC, multi_UBspline_2d_d::x_grid, multi_UBspline_2d_d::xBC, multi_UBspline_2d_d::y_grid, multi_UBspline_2d_d::y_stride, and multi_UBspline_2d_d::yBC.

set_multi_UBspline_2d_s()

void set_multi_UBspline_2d_s	(	multi_UBspline_2d_s *	spline,
		int	num,
		float *	data )

Definition at line 215 of file multi_bspline_create.cpp.

{
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Nx, Ny;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
 
  float *coefs = spline->coefs + num;
  int ys = spline->y_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    intptr_t doffset = iy;
    intptr_t coffset = iy*ys;
    find_coefs_1d_s (spline->x_grid, spline->xBC, data+doffset, (intptr_t)My,
             coefs+coffset, (intptr_t)Ny*ys);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    intptr_t doffset = ix*Ny*ys;
    intptr_t coffset = ix*Ny*ys;
    find_coefs_1d_s (spline->y_grid, spline->yBC, coefs+doffset, (intptr_t)ys, 
             coefs+coffset, (intptr_t)ys);
  }
}

References ANTIPERIODIC, multi_UBspline_2d_s::coefs, find_coefs_1d_s(), BCtype_s::lCode, Ugrid::num, PERIODIC, multi_UBspline_2d_s::x_grid, multi_UBspline_2d_s::xBC, multi_UBspline_2d_s::y_grid, multi_UBspline_2d_s::y_stride, and multi_UBspline_2d_s::yBC.

set_multi_UBspline_2d_z()

void set_multi_UBspline_2d_z	(	multi_UBspline_2d_z *	spline,
		int	num,
		complex_double *	data )

Definition at line 1257 of file multi_bspline_create.cpp.

{
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Nx, Ny;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
 
  BCtype_d xBC_r, xBC_i, yBC_r, yBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
  yBC_r.lCode = spline->yBC.lCode;  yBC_r.rCode = spline->yBC.rCode;
  yBC_r.lVal  = spline->yBC.lVal_r; yBC_r.rVal  = spline->yBC.rVal_r;
  yBC_i.lCode = spline->yBC.lCode;  yBC_i.rCode = spline->yBC.rCode;
  yBC_i.lVal  = spline->yBC.lVal_i; yBC_i.rVal  = spline->yBC.rVal_i;
 
  complex_double *coefs = spline->coefs + num;
  int ys = spline->y_stride;
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) {
    intptr_t doffset = 2*iy;
    intptr_t coffset = 2*iy*ys;
    // Real part
    find_coefs_1d_d (spline->x_grid, xBC_r, 
             ((double*)data+doffset), (intptr_t)2*My, 
             (double*)coefs+coffset,  (intptr_t)2*Ny*ys);
    // Imag part
    find_coefs_1d_d (spline->x_grid, xBC_i, 
             ((double*)data)+doffset+1,  (intptr_t)2*My, 
             ((double*)coefs)+coffset+1, (intptr_t)2*Ny*ys);
  }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) {
    intptr_t doffset = 2*ix*Ny*ys;
    intptr_t coffset = 2*ix*Ny*ys;
    // Real part
    find_coefs_1d_d (spline->y_grid, yBC_r, 
             ((double*)coefs)+doffset, (intptr_t)2*ys, 
             (double*)coefs+coffset,   (intptr_t)2*ys);
    // Imag part
    find_coefs_1d_d (spline->y_grid, yBC_i, 
             (double*)coefs+doffset+1,   (intptr_t)2*ys, 
             ((double*)coefs)+coffset+1, (intptr_t)2*ys);
  }
}

References ANTIPERIODIC, multi_UBspline_2d_z::coefs, find_coefs_1d_d(), BCtype_d::lCode, BCtype_z::lCode, BCtype_d::lVal, BCtype_z::lVal_i, BCtype_z::lVal_r, Ugrid::num, PERIODIC, BCtype_d::rCode, BCtype_z::rCode, BCtype_d::rVal, BCtype_z::rVal_i, BCtype_z::rVal_r, multi_UBspline_2d_z::x_grid, multi_UBspline_2d_z::xBC, multi_UBspline_2d_z::y_grid, multi_UBspline_2d_z::y_stride, and multi_UBspline_2d_z::yBC.

set_multi_UBspline_3d_c()

void set_multi_UBspline_3d_c	(	multi_UBspline_3d_c *	spline,
		int	num,
		complex_float *	data )

Definition at line 650 of file multi_bspline_create.cpp.

{
  // Setup internal variables
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Mz = spline->z_grid.num;
  int Nx, Ny, Nz;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                           
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                           
    Ny = My+2;
  if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                          
    Nz = Mz+2;
 
  BCtype_s xBC_r, xBC_i, yBC_r, yBC_i, zBC_r, zBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
  yBC_r.lCode = spline->yBC.lCode;  yBC_r.rCode = spline->yBC.rCode;
  yBC_r.lVal  = spline->yBC.lVal_r; yBC_r.rVal  = spline->yBC.rVal_r;
  yBC_i.lCode = spline->yBC.lCode;  yBC_i.rCode = spline->yBC.rCode;
  yBC_i.lVal  = spline->yBC.lVal_i; yBC_i.rVal  = spline->yBC.rVal_i;
  zBC_r.lCode = spline->zBC.lCode;  zBC_r.rCode = spline->zBC.rCode;
  zBC_r.lVal  = spline->zBC.lVal_r; zBC_r.rVal  = spline->zBC.rVal_r;
  zBC_i.lCode = spline->zBC.lCode;  zBC_i.rCode = spline->zBC.rCode;
  zBC_i.lVal  = spline->zBC.lVal_i; zBC_i.rVal  = spline->zBC.rVal_i;
 
  complex_float *coefs = spline->coefs + num;
  int zs = spline->z_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      intptr_t doffset = 2*(iy*Mz+iz);
      intptr_t coffset = 2*(iy*Nz+iz)*zs;
      // Real part
      find_coefs_1d_s (spline->x_grid, xBC_r, 
               ((float*)data)+doffset,  (intptr_t)2*My*Mz,
               ((float*)coefs)+coffset, (intptr_t)2*Ny*Nz*zs);
      // Imag part
      find_coefs_1d_s (spline->x_grid, xBC_i, 
               ((float*)data)+doffset+1,  (intptr_t)2*My*Mz,
               ((float*)coefs)+coffset+1, (intptr_t)2*Ny*Nz*zs);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      intptr_t doffset = 2*(ix*Ny*Nz + iz)*zs;
      intptr_t coffset = 2*(ix*Ny*Nz + iz)*zs;
      // Real part
      find_coefs_1d_s (spline->y_grid, yBC_r, 
               ((float*)coefs)+doffset, (intptr_t)2*Nz*zs, 
               ((float*)coefs)+coffset, (intptr_t)2*Nz*zs);
      // Imag part
      find_coefs_1d_s (spline->y_grid, yBC_i, 
               ((float*)coefs)+doffset+1, (intptr_t)2*Nz*zs, 
               ((float*)coefs)+coffset+1, (intptr_t)2*Nz*zs);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      intptr_t doffset = 2*((ix*Ny+iy)*Nz)*zs;
      intptr_t coffset = 2*((ix*Ny+iy)*Nz)*zs;
      // Real part
      find_coefs_1d_s (spline->z_grid, zBC_r, 
               ((float*)coefs)+doffset, (intptr_t)2*zs, 
               ((float*)coefs)+coffset, (intptr_t)2*zs);
      // Imag part
      find_coefs_1d_s (spline->z_grid, zBC_i, 
               ((float*)coefs)+doffset+1, (intptr_t)2*zs, 
               ((float*)coefs)+coffset+1, (intptr_t)2*zs);
    }
}

References ANTIPERIODIC, multi_UBspline_3d_c::coefs, find_coefs_1d_s(), BCtype_s::lCode, BCtype_c::lCode, BCtype_s::lVal, BCtype_c::lVal_i, BCtype_c::lVal_r, Ugrid::num, PERIODIC, BCtype_s::rCode, BCtype_c::rCode, BCtype_s::rVal, BCtype_c::rVal_i, BCtype_c::rVal_r, multi_UBspline_3d_c::x_grid, multi_UBspline_3d_c::xBC, multi_UBspline_3d_c::y_grid, multi_UBspline_3d_c::yBC, multi_UBspline_3d_c::z_grid, multi_UBspline_3d_c::z_stride, and multi_UBspline_3d_c::zBC.

set_multi_UBspline_3d_d()

void set_multi_UBspline_3d_d	(	multi_UBspline_3d_d *	spline,
		int	num,
		double *	data )

Definition at line 1020 of file multi_bspline_create.cpp.

{
  int Mx = spline->x_grid.num;  
  int My = spline->y_grid.num; 
  int Mz = spline->z_grid.num;
  int Nx, Ny, Nz;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
  if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                                   
    Nz = Mz+2;
 
  double *coefs = spline->coefs + num;
  intptr_t zs = spline->z_stride;
 
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      intptr_t doffset = iy*Mz+iz;
      intptr_t coffset = (iy*Nz+iz)*zs;
      find_coefs_1d_d (spline->x_grid, spline->xBC, 
               data+doffset,  (intptr_t)My*Mz,
               coefs+coffset, (intptr_t)Ny*Nz*zs);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      intptr_t doffset = (ix*Ny*Nz + iz)*zs;
      intptr_t coffset = (ix*Ny*Nz + iz)*zs;
      find_coefs_1d_d (spline->y_grid, spline->yBC, 
               coefs+doffset, (intptr_t)Nz*zs, 
               coefs+coffset, (intptr_t)Nz*zs);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      intptr_t doffset = (ix*Ny+iy)*Nz*zs;
      intptr_t coffset = (ix*Ny+iy)*Nz*zs;
      find_coefs_1d_d (spline->z_grid, spline->zBC, 
               coefs+doffset, (intptr_t)zs, 
               coefs+coffset, (intptr_t)zs);
    }
}

References ANTIPERIODIC, multi_UBspline_3d_d::coefs, find_coefs_1d_d(), BCtype_d::lCode, Ugrid::num, PERIODIC, multi_UBspline_3d_d::x_grid, multi_UBspline_3d_d::xBC, multi_UBspline_3d_d::y_grid, multi_UBspline_3d_d::yBC, multi_UBspline_3d_d::z_grid, multi_UBspline_3d_d::z_stride, and multi_UBspline_3d_d::zBC.

set_multi_UBspline_3d_s()

void set_multi_UBspline_3d_s	(	multi_UBspline_3d_s *	spline,
		int	num,
		float *	data )

Definition at line 323 of file multi_bspline_create.cpp.

{
  int Mx = spline->x_grid.num;
  int My = spline->y_grid.num;
  int Mz = spline->z_grid.num;
  int Nx, Ny, Nz;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
  if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                                   
    Nz = Mz+2;
 
  float *coefs = spline->coefs + num;
 
  int zs = spline->z_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      intptr_t doffset = iy*Mz+iz;
      intptr_t coffset = (iy*Nz+iz)*zs;
      find_coefs_1d_s (spline->x_grid, spline->xBC, data+doffset, (intptr_t)My*Mz,
               coefs+coffset, (intptr_t)Ny*Nz*zs);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      intptr_t doffset = (ix*Ny*Nz + iz)*zs;
      intptr_t coffset = (ix*Ny*Nz + iz)*zs;
      find_coefs_1d_s (spline->y_grid, spline->yBC, coefs+doffset, (intptr_t)Nz*zs, 
               coefs+coffset, (intptr_t)Nz*zs);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      intptr_t doffset = ((ix*Ny+iy)*Nz)*zs;
      intptr_t coffset = ((ix*Ny+iy)*Nz)*zs;
      find_coefs_1d_s (spline->z_grid, spline->zBC, coefs+doffset, 
               (intptr_t)zs, coefs+coffset, (intptr_t)zs);
    }
}

References ANTIPERIODIC, multi_UBspline_3d_s::coefs, find_coefs_1d_s(), BCtype_s::lCode, Ugrid::num, PERIODIC, multi_UBspline_3d_s::x_grid, multi_UBspline_3d_s::xBC, multi_UBspline_3d_s::y_grid, multi_UBspline_3d_s::yBC, multi_UBspline_3d_s::z_grid, multi_UBspline_3d_s::z_stride, and multi_UBspline_3d_s::zBC.

set_multi_UBspline_3d_z()

void set_multi_UBspline_3d_z	(	multi_UBspline_3d_z *	spline,
		int	num,
		complex_double *	data )

Definition at line 1393 of file multi_bspline_create.cpp.

{
  // Setup internal variables
  int Mx = spline->x_grid.num;  
  int My = spline->y_grid.num; 
  int Mz = spline->z_grid.num;
  int Nx, Ny, Nz;
 
  if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)     
    Nx = Mx+3;
  else                                   
    Nx = Mx+2;
  if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)     
    Ny = My+3;
  else                                   
    Ny = My+2;
  if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)     
    Nz = Mz+3;
  else                                   
    Nz = Mz+2;
 
  BCtype_d xBC_r, xBC_i, yBC_r, yBC_i, zBC_r, zBC_i;
  xBC_r.lCode = spline->xBC.lCode;  xBC_r.rCode = spline->xBC.rCode;
  xBC_r.lVal  = spline->xBC.lVal_r; xBC_r.rVal  = spline->xBC.rVal_r;
  xBC_i.lCode = spline->xBC.lCode;  xBC_i.rCode = spline->xBC.rCode;
  xBC_i.lVal  = spline->xBC.lVal_i; xBC_i.rVal  = spline->xBC.rVal_i;
  yBC_r.lCode = spline->yBC.lCode;  yBC_r.rCode = spline->yBC.rCode;
  yBC_r.lVal  = spline->yBC.lVal_r; yBC_r.rVal  = spline->yBC.rVal_r;
  yBC_i.lCode = spline->yBC.lCode;  yBC_i.rCode = spline->yBC.rCode;
  yBC_i.lVal  = spline->yBC.lVal_i; yBC_i.rVal  = spline->yBC.rVal_i;
  zBC_r.lCode = spline->zBC.lCode;  zBC_r.rCode = spline->zBC.rCode;
  zBC_r.lVal  = spline->zBC.lVal_r; zBC_r.rVal  = spline->zBC.rVal_r;
  zBC_i.lCode = spline->zBC.lCode;  zBC_i.rCode = spline->zBC.rCode;
  zBC_i.lVal  = spline->zBC.lVal_i; zBC_i.rVal  = spline->zBC.rVal_i;
 
  complex_double *coefs = spline->coefs + num;
 
// unused variable
//  int N = spline->num_splines;
  int zs = spline->z_stride;
  // First, solve in the X-direction 
  for (int iy=0; iy<My; iy++) 
    for (int iz=0; iz<Mz; iz++) {
      intptr_t doffset = 2*(iy*Mz+iz);
      intptr_t coffset = 2*(iy*Nz+iz)*zs;
      // Real part
      find_coefs_1d_d (spline->x_grid, xBC_r, 
               ((double*)data)+doffset,  (intptr_t)2*My*Mz,
               ((double*)coefs)+coffset, (intptr_t)2*Ny*Nz*zs);
      // Imag part
      find_coefs_1d_d (spline->x_grid, xBC_i, 
               ((double*)data)+doffset+1,  (intptr_t)2*My*Mz,
               ((double*)coefs)+coffset+1, (intptr_t)2*Ny*Nz*zs);
    }
  
  // Now, solve in the Y-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iz=0; iz<Nz; iz++) {
      intptr_t doffset = 2*(ix*Ny*Nz + iz)*zs;
      intptr_t coffset = 2*(ix*Ny*Nz + iz)*zs;
      // Real part
      find_coefs_1d_d (spline->y_grid, yBC_r, 
               ((double*)coefs)+doffset, (intptr_t)2*Nz*zs, 
               ((double*)coefs)+coffset, (intptr_t)2*Nz*zs);
      // Imag part
      find_coefs_1d_d (spline->y_grid, yBC_i, 
               ((double*)coefs)+doffset+1, (intptr_t)2*Nz*zs, 
               ((double*)coefs)+coffset+1, (intptr_t)2*Nz*zs);
    }
 
  // Now, solve in the Z-direction
  for (int ix=0; ix<Nx; ix++) 
    for (int iy=0; iy<Ny; iy++) {
      intptr_t doffset = 2*((ix*Ny+iy)*Nz)*zs;
      intptr_t coffset = 2*((ix*Ny+iy)*Nz)*zs;
      // Real part
      find_coefs_1d_d (spline->z_grid, zBC_r, 
               ((double*)coefs)+doffset, (intptr_t)2*zs, 
               ((double*)coefs)+coffset, (intptr_t)2*zs);
      // Imag part
      find_coefs_1d_d (spline->z_grid, zBC_i, 
               ((double*)coefs)+doffset+1, (intptr_t)2*zs, 
               ((double*)coefs)+coffset+1, (intptr_t)2*zs);
    }
}

References ANTIPERIODIC, multi_UBspline_3d_z::coefs, find_coefs_1d_d(), BCtype_d::lCode, BCtype_z::lCode, BCtype_d::lVal, BCtype_z::lVal_i, BCtype_z::lVal_r, Ugrid::num, PERIODIC, BCtype_d::rCode, BCtype_z::rCode, BCtype_d::rVal, BCtype_z::rVal_i, BCtype_z::rVal_r, multi_UBspline_3d_z::x_grid, multi_UBspline_3d_z::xBC, multi_UBspline_3d_z::y_grid, multi_UBspline_3d_z::yBC, multi_UBspline_3d_z::z_grid, multi_UBspline_3d_z::z_stride, and multi_UBspline_3d_z::zBC.

solve_antiperiodic_interp_1d_s()

void solve_antiperiodic_interp_1d_s	(	float	bands[],
		float	coefs[],
		int	M,
		int	cstride )

Definition at line 156 of file bspline_create.cpp.

{
  bands[4*0+0]     *= -1.0;
  bands[4*(M-1)+2] *= -1.0;
 
  std::vector<float> lastCol(M);
 
  // 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];
  coefs[M*cstride] = bands[4*(M-1)+3];
  for (int row=M-2; row>=0; row--) 
    coefs[(row+1)*cstride] = 
      bands[4*(row)+3] - bands[4*(row)+2]*coefs[(row+2)*cstride] - lastCol[row]*coefs[M*cstride];
  
  coefs[0*cstride]     = -coefs[M*cstride];
  coefs[(M+1)*cstride] = -coefs[1*cstride];
  coefs[(M+2)*cstride] = -coefs[2*cstride];
 
 
}

solve_deriv_interp_1d_d()

void solve_deriv_interp_1d_d	(	double	bands[],
		double	coefs[],
		int	M,
		int	cstride )

Definition at line 959 of file bspline_create.cpp.

{
  // Solve interpolating equations
  // 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 < (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];
    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;
 
  coefs[(M+1)*cstride] = bands[4*(M+1)+3];
  // Now back substitute up
  for (int row=M; row>0; row--)
    coefs[row*cstride] = bands[4*(row)+3] - bands[4*(row)+2]*coefs[cstride*(row+1)];
  
  // Finish with first row
  coefs[0] = bands[4*(0)+3] - bands[4*(0)+1]*coefs[1*cstride] - bands[4*(0)+2]*coefs[2*cstride];
}

solve_deriv_interp_1d_s()

void solve_deriv_interp_1d_s	(	float	bands[],
		float	coefs[],
		int	M,
		int	cstride )

Definition at line 47 of file bspline_create.cpp.

{
  // Solve interpolating equations
  // 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 < (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];
    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;
 
  coefs[(M+1)*cstride] = bands[4*(M+1)+3];
  // Now back substitute up
  for (int row=M; row>0; row--)
    coefs[row*cstride] = bands[4*(row)+3] - bands[4*(row)+2]*coefs[cstride*(row+1)];
  
  // Finish with first row
  coefs[0] = bands[4*(0)+3] - bands[4*(0)+1]*coefs[1*cstride] - bands[4*(0)+2]*coefs[2*cstride];
}

solve_periodic_interp_1d_d()

void solve_periodic_interp_1d_d	(	double	bands[],
		double	coefs[],
		int	M,
		int	cstride )

Definition at line 1010 of file bspline_create.cpp.

{
  std::vector<double> lastCol(M);
 
  // 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];
  coefs[M*cstride] = bands[4*(M-1)+3];
  for (int row=M-2; row>=0; row--) 
    coefs[(row+1)*cstride] = 
      bands[4*(row)+3] - bands[4*(row)+2]*coefs[(row+2)*cstride] - lastCol[row]*coefs[M*cstride];
  
  coefs[0*cstride] = coefs[M*cstride];
  coefs[(M+1)*cstride] = coefs[1*cstride];
  coefs[(M+2)*cstride] = coefs[2*cstride];
}

solve_periodic_interp_1d_s()

void solve_periodic_interp_1d_s	(	float	bands[],
		float	coefs[],
		int	M,
		int	cstride )

Definition at line 98 of file bspline_create.cpp.

{
  std::vector<float> lastCol(M);
 
  // 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];
  coefs[M*cstride] = bands[4*(M-1)+3];
  for (int row=M-2; row>=0; row--) 
    coefs[(row+1)*cstride] = 
      bands[4*(row)+3] - bands[4*(row)+2]*coefs[(row+2)*cstride] - lastCol[row]*coefs[M*cstride];
  
  coefs[0*cstride] = coefs[M*cstride];
  coefs[(M+1)*cstride] = coefs[1*cstride];
  coefs[(M+2)*cstride] = coefs[2*cstride];
 
 
}