initial checkin of new LocalInterp thorn

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@2 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 1150 insertions, 0 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
new file mode 100644
index 0000000..bcf34e5
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -0,0 +1,1150 @@
+/*@@
+  @file      template.c
+  @date      23 Oct 2001
+  @author    Jonathan Thornburg <jthorn@aei.mpg.de>
+  @desc
+	This file is a template for generalized interpolation for
+	1-, 2-, or 3-d molecules.  It's used by defining various
+	C preprocessor macros, then #including this file.  Each
+	such inclusion defines a single interpolation function,
+	which does generalized interpolation for a single
+	combination of (N_dims, molecule_family, order, smoothing).
+
+	The following header files must be #included before
+	#including this file:
+		<math.h>
+		<limits.h>
+		<stdlib.h>
+		<string.h>
+		<stdio.h>
+		"util_ErrorCodes.h"
+		"cctk.h"
+		"InterpLocalUniform.h"
+
+	The following preprocessor macros must be defined before
+	#including this file (note the macros which are file names
+	must all include the proper quotes for a #include in their
+	expansion, e.g.
+		#define DATA_VAR_DCL_FILE_NAME	\
+			"coeffs/2d.cube.size3/data-var.dcl.c"
+  @enddesc
+
+  @var	    FUNCTION_NAME
+  @vdesc    The name of the interpolation function, e.g.
+		#define FUNCTION_NAME	LocalInterp_ILA_2d_cube_ord2_sm0
+  @endvar
+
+  @var	    N_DIMS
+  @vdesc    The number of dimensions in which to interpolate, e.g.
+		#define N_DIMS	2
+	    The present implementation restricts this to 1, 2,
+	    or 3, but this could easily be changed if needed.
+	    Note that MAX_N_DIMS (defined in "InterpLocalUniform.c")
+	    is a compile-time upper bound for N_DIMS, useful for
+	    sizing arrays etc.
+  @endvar
+
+  @var	    MOLECULE_SIZE
+  @vdesc    The diameter of (number of points in) the molecules to be used,
+	    e.g.
+		#define MOLECULE_SIZE	3
+  	    The present implementation takes this to be the same in each
+	    dimension, but this could easily be changed if needed.  
+  @endvar
+
+  @var	    HAVE_OP_{I,DX,DY,DXX,DXY,DYY,...}
+	    Each of these symbols should be defined or not defined
+	    according as if the corresponding derivative operator is
+	    to be supported or not supported by this function.
+  @endvar
+
+  @var	    DATA_VAR_DCL_FILE_NAME
+  @vdesc    The name of a file (presumably machine-generated) containing
+	    a sequence of one or more C declarations for a molecule-sized
+	    set of "data variables", eg (for 2D size-3 molecules)
+		fp data_m1_m1;
+		fp data_0_m1;
+		fp data_p1_m1;
+		fp data_m1_0;
+		fp data_0_0;
+		fp data_p1_0;
+		fp data_m1_p1;
+		fp data_0_p1;
+		fp data_p1_p1;
+  @endvar
+
+  @var	    COEFF_{I,DX,DY,DXX,DXY,DYY,...}_DCL_FILE_NAME
+  @vdesc    Each of these macros should be the name of a file
+	    (presumably machine-generated) containing a sequence
+	    of one or more C declarations for a molecule-sized set
+	    of coefficients for the corresponding derivative operator,
+	    eg (for dx with 2D size-3 molecules)
+		fp coeff_dx_m1_m1;
+		fp coeff_dx_0_m1;
+		fp coeff_dx_p1_m1;
+		fp coeff_dx_m1_0;
+		fp coeff_dx_0_0;
+		fp coeff_dx_p1_0;
+		fp coeff_dx_m1_p1;
+		fp coeff_dx_0_p1;
+		fp coeff_dx_p1_p1;
+  @endvar
+
+  @var	    DATA_VAR_ASSIGN_FILE_NAME
+  @vdesc    The name of a file (presumably machine-generated) containing
+	    a sequence of C assignment statements to assign DATA(...) to
+	    the corresponding data variables for each molecule component,
+	    eg (for 2D size-3 molecules)
+		data_m1_m1 = DATA(-1,-1);
+		data_0_m1 = DATA(0,-1);
+		data_p1_m1 = DATA(1,-1);
+		data_m1_0 = DATA(-1,0);
+		data_0_0 = DATA(0,0);
+		data_p1_0 = DATA(1,0);
+		data_m1_p1 = DATA(-1,+1);
+		data_0_p1 = DATA(0,+1);
+		data_p1_p1 = DATA(1,+1);
+  @endvar
+
+  @var	    INTERP_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME
+  @vdesc    Each of these macros should be the name of a file
+	    (presumably machine-generated) containing a C assignment
+	    statement (or a sequence of such statements) to compute
+	    the variable  result as the molecule-sized linear combination
+	    of the data variables (cf DATA_VAR_{DCL,ASSIGN}_FILE_NAME)
+	    and the interpolation coefficients
+	    (cf COEFF_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME),
+	    eg (for 2D size-3 molecules, dx operator)
+		result =
+			coeff_dx_m1_m1*data_m1_m1
+			+ coeff_dx_0_m1*data_0_m1
+			+ coeff_dx_p1_m1*data_p1_m1
+			+ coeff_dx_p2_m1*data_p2_m1
+			+ coeff_dx_m1_0*data_m1_0
+			+ coeff_dx_0_0*data_0_0
+			+ coeff_dx_p1_0*data_p1_0
+			+ coeff_dx_p2_0*data_p2_0
+			+ coeff_dx_m1_p1*data_m1_p1
+			+ coeff_dx_0_p1*data_0_p1
+			+ coeff_dx_p1_p1*data_p1_p1
+			+ coeff_dx_p2_p1*data_p2_p1
+			+ coeff_dx_m1_p2*data_m1_p2
+			+ coeff_dx_0_p2*data_0_p2
+			+ coeff_dx_p1_p2*data_p1_p2
+			+ coeff_dx_p2_p2*data_p2_p2;
+	    Note that this may *NOT* include any variable declarations;
+	    it should be valid to appear in the middle of a sequence of
+	    C statements.
+  @endvar
+
+  @var	    COEFF_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME
+  @vdesc    Each of these macros should be the name of a file
+	    (presumably machine-generated) containing a sequence of C
+	    assignment statements to compute all the (molecule-sized
+	    set of) interpolation coefficients for the specified
+	    derivative operator, as polynomials in the variables
+	    (x,y,z), eg (for 2D size-3 molecules, dx operator)
+		fp t35,
+		   t42,
+		   t41,
+		   t40,
+		   t39,
+		   t36;
+		      t35 = RATIONAL(1.0,3.0)*x;
+		      t42 = t35+RATIONAL(-1.0,4.0)*y;
+		      t41 = t35+RATIONAL(1.0,4.0)*y;
+		      t40 = RATIONAL(-1.0,6.0);
+		      t39 = RATIONAL(1.0,6.0);
+		      t36 = RATIONAL(-2.0,3.0)*x;
+		      coeff_dx_m1_m1 = t40+t41;
+		      coeff_dx_0_m1 = t36;
+		      coeff_dx_p1_m1 = t39+t42;
+		      coeff_dx_m1_0 = t40+t35;
+		      coeff_dx_0_0 = t36;
+		      coeff_dx_p1_0 = t35+t39;
+		      coeff_dx_m1_p1 = t40+t42;
+		      coeff_dx_0_p1 = t36;
+		      coeff_dx_p1_p1 = t39+t41;
+
+	    As illustrated, the code may use the macro RATIONAL
+	    (defined later in this file) to represent rational-number
+	    coefficients, and it may also declare temporary variables
+	    as needed.
+
+	    Note that these are the only coefficients which depend
+	    on the actual interpolation scheme used; all the others
+	    just depend on the interpolation dimension and molecule
+	    family and size.
+  @endvar
+
+  @version   $Id$
+@@*/
+
+/******************************************************************************/
+
+/*@@
+   @routine    FUNCTION_NAME
+   @date       23 Oct 2001
+   @author     Jonathan Thornburg <jthorn@aei.mpg.de>
+   @desc
+	This function does generalized interpolation of one or more
+	2d arrays to arbitrary points.  For details, see the header
+	comments for InterpLocalUniform() (in "InterpLocalUniform.c"
+	in this same directory).
+
+	This function's arguments are all a subset of those of
+	InterpLocalUniform() ; the only difference is that this function
+	takes all its arguments explicitly, whereas  InputLocalArrays()
+	takes some of them indirectly via a key/value parameter table.
+
+  @returntype   int
+  @returndesc	This function's return result is the same as that of
+		InterpLocalUniform():
+		0			success
+		UTIL_ERROR_BAD_INPUT	one of the input arguments is invalid
+		UTIL_ERROR_NO_MEMORY	unable to malloc() temporary memory
+		CCTK_ERROR_INTERP_POINT_X_RANGE
+					interpolation point is out of range
+  @endreturndesc
+
+  @@*/
+int FUNCTION_NAME(/***** coordinate system *****/
+		  const CCTK_REAL coord_system_origin[],
+		  const CCTK_REAL grid_spacing[],
+		  /***** interpolation points *****/
+		  int N_interp_points,
+		  int interp_coords_type_code,
+		  const void *const interp_coords[],
+		  const CCTK_REAL out_of_range_tolerance[],
+		  /***** input arrays *****/
+		  int N_input_arrays,
+		  const CCTK_INT input_array_offsets[],
+		  const CCTK_INT input_array_strides[],
+		  const CCTK_INT input_array_min_subscripts[],
+		  const CCTK_INT input_array_max_subscripts[],
+		  const CCTK_INT input_array_type_codes[],
+		  const void *const input_arrays[],
+		  /***** output arrays *****/
+		  int N_output_arrays,
+		  const CCTK_INT output_array_type_codes[],
+		  void *const output_arrays[],
+		  /***** operation info */
+		  const CCTK_INT operand_indices[],
+		  const CCTK_INT operation_codes[],
+		  /***** error-reporting *****/
+		  int* error_pt, int* error_axis, int* error_end)
+{
+/*
+ * Implementation notes:
+ * 
+ * The basic outline of this function is as follows:
+ *
+ *	compute "which derivatives are wanted" flags
+ *	precompute 1/dx factors
+ *		for each interpolation point
+ *		{
+ *		declare all the coefficients
+ *		declare all the data-values variables
+ *		***fetch*** interpolation point coordinates
+ *		compute coefficients for all derivatives which are wanted
+ *			for each output array
+ *			{
+ *			***decode*** the input/output array datatypes
+ *			   to determine whether they're real or complex
+ *			   (they must both be the same in this regard),
+ *			   and define
+ *			      int N_parts = data is complex ? 2 : 1;
+ *			int part;
+ *				for (part = 0 ; part < N_parts ; ++part)
+ *				{
+ *				if (this output array is computed
+ *				    using a different input array
+ *				    than the previous one || part != 0)
+ *				   then ***fetch*** the input array values
+ *					            into local "data" variables
+ *				  {
+ *				fp result;
+ *				switch	(operation_code)
+ *					{
+ *				case 0:
+ *					result = compute the interpolant
+ *						 itself as a linear combination
+ *						 of data variables & op=0 coeffs
+ *					break;
+ *				case 1:
+ *					result = compute the interpolant
+ *						 itself as a linear combination
+ *						 of data variables & op=1 coeffs
+ *					result *= inverse_dx;
+ *					break;
+ *				case ...
+ *					}
+ *				***store*** result in output array
+ *				  }
+ *				}
+ *			}
+ *		}
+ *
+ * Here "***decode***", "***fetch***", and "***store***" are all
+ * actually switches on the various array datatypes.  For complex
+ * datatypes "***fetch***" and "***store***" pointer-alias offset the
+ * N-dimensional input array to a 2-dimensional array where the 1st axis
+ * is the 1-D subscript corresponding to the N input axes, and the 2nd
+ * axes has 2 elements subscripted by [part] for the (real,imaginary)
+ * components of the complex values.
+ *
+ * At present we do all floating-point computations in type "fp"
+ * (typically a typedef for CCTK_REAL), so arrays of higher precision
+ * than this will incur extra rounding errors.  In practice these should
+ * be negligible compared to the "truncation" interpolation errors.
+ */
+
+/*
+ * Naming conventions:
+ * in, out = 0-origin indices each selecting an input/output array
+ * pt = 0-origin index selecting an interpolation point
+ */
+
+/* number of real "parts" in a complex number */
+#define COMPLEX_N_PARTS	2
+
+/* layout of axes in N_dims-element arrays */
+#define X_AXIS	0
+#define Y_AXIS	1
+#define Z_AXIS	2
+#if (N_DIMS > 3)
+  #error "N_DIMS may not be > 3!"
+#endif
+
+/* input array size, strides, and subscripting computation */
+#if (N_DIMS >= 1)
+  const int stride_i = input_array_strides[X_AXIS];
+#endif
+#if (N_DIMS >= 2)
+  const int stride_j = input_array_strides[Y_AXIS];
+#endif
+#if (N_DIMS >= 3)
+  const int stride_k = input_array_strides[Z_AXIS];
+#endif
+
+#if   (N_DIMS == 1)
+  #define SUB1(i)	(i*stride_i)
+#elif (N_DIMS == 2)
+  #define SUB2(i,j)	(i*stride_i + j*stride_j)
+#elif (N_DIMS == 3)
+  #define SUB3(i,j,k)	(i*stride_i + j*stride_j + k*stride_k)
+#else
+  #error "N_DIMS must be 1, 2, or 3!"
+#endif
+
+/* macros used by machine-generated interpolation coefficient expressions */
+/*
+ * FIXME: right now this is used as (eg) RATIONAL(1.0,2.0);
+ *	  it might be cleaner if it were RATIONAL(1,2) with the
+ *	  preprocessor ## operator used to glue on the .0
+ *	  (I _think_ this is portable, but is it really??)
+ */
+#define RATIONAL(num,den)	(num/den)
+
+
+/*
+ * compute flags specifying which derivatives are wanted
+ */
+#ifdef HAVE_OP_I
+  bool want_I = false;
+#endif
+#ifdef HAVE_OP_DX
+  bool want_dx = false;
+#endif
+#ifdef HAVE_OP_DY
+  bool want_dy = false;
+#endif
+#ifdef HAVE_OP_DZ
+  bool want_dz = false;
+#endif
+#ifdef HAVE_OP_DXX
+  bool want_dxx = false;
+#endif
+#ifdef HAVE_OP_DXY
+  bool want_dxy = false;
+#endif
+#ifdef HAVE_OP_DXZ
+  bool want_dxz = false;
+#endif
+#ifdef HAVE_OP_DYY
+  bool want_dyy = false;
+#endif
+#ifdef HAVE_OP_DYZ
+  bool want_dyz = false;
+#endif
+#ifdef HAVE_OP_DZZ
+  bool want_dzz = false;
+#endif
+  {
+int out;
+	for (out = 0 ; out < N_output_arrays ; ++out)
+	{
+	switch	(operation_codes[out])
+		{
+	#ifdef HAVE_OP_I
+	  case 0:	want_I = true;		break;
+	#endif
+	#ifdef HAVE_OP_DX
+	  case 1:	want_dx = true;		break;
+	#endif
+	#ifdef HAVE_OP_DY
+	  case 2:	want_dy = true;		break;
+	#endif
+	#ifdef HAVE_OP_DZ
+	  case 3:	want_dz = true;		break;
+	#endif
+	#ifdef HAVE_OP_DXX
+	  case 11:	want_dxx = true;	break;
+	#endif
+	#ifdef HAVE_OP_DXY
+	  case 12:
+	  case 21:	want_dxy = true;	break;
+	#endif
+	#ifdef HAVE_OP_DXZ
+	  case 13:
+	  case 31:	want_dxz = true;	break;
+	#endif
+	#ifdef HAVE_OP_DYY
+	  case 22:	want_dyy = true;	break;
+	#endif
+	#ifdef HAVE_OP_DYZ
+	  case 23:
+	  case 32:	want_dyz = true;	break;
+	#endif
+	#ifdef HAVE_OP_DZZ
+	  case 33:	want_dzz = true;	break;
+	#endif
+	default:
+CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "Generalized interpolation operation_code %d not supported!",
+	   operation_codes[out]);				/*NOTREACHED*/
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		}
+	}
+  }
+
+/*
+ * save origin/delta variables, precompute 1/delta factors
+ * ... in theory we could only compute those factors we're going to use,
+ *     but it's not worth the trouble, so we just compute them all
+ */
+  {
+#if N_DIMS >= 1
+  const fp origin_x = coord_system_origin[X_AXIS];
+  const fp delta_x = grid_spacing[X_AXIS];
+  #if defined(HAVE_OP_DX) || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ)
+    const fp inverse_delta_x = 1.0 / delta_x;
+  #endif
+  #ifdef HAVE_OP_DXX
+    const fp inverse_delta_x2 = 1.0 / (delta_x*delta_x);
+  #endif
+#endif
+#if N_DIMS >= 2
+  const fp origin_y = coord_system_origin[Y_AXIS];
+  const fp delta_y = grid_spacing[Y_AXIS];
+  #if defined(HAVE_OP_DY) || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ)
+    const fp inverse_delta_y = 1.0 / delta_y;
+  #endif
+  #ifdef HAVE_OP_DYY
+    const fp inverse_delta_y2 = 1.0 / (delta_y*delta_y);
+  #endif
+#endif
+#if N_DIMS >= 3
+  const fp origin_z = coord_system_origin[Z_AXIS];
+  const fp delta_z = grid_spacing[Z_AXIS];
+  #if defined(HAVE_OP_DZ) || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ)
+    const fp inverse_delta_z = 1.0 / delta_z;
+  #endif
+  #ifdef HAVE_OP_DZZ
+    const fp inverse_delta_z2 = 1.0 / (delta_z*delta_z);
+  #endif
+#endif
+
+/*
+ * interpolate at each point
+ */
+  {
+int pt;
+	for (pt = 0 ; pt < N_interp_points ; ++pt)
+	{
+	/* declare all the data-values variables */
+	#include DATA_VAR_DCL_FILE_NAME
+
+	/* declare all the interpolation coefficients */
+	#ifdef HAVE_OP_I
+	  #include COEFF_I_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DX
+	  #include COEFF_DX_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DY
+	  #include COEFF_DY_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DZ
+	  #include COEFF_DZ_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DXX
+	  #include COEFF_DXX_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DXY
+	  #include COEFF_DXY_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DXZ
+	  #include COEFF_DXZ_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DYY
+	  #include COEFF_DYY_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DYZ
+	  #include COEFF_DYZ_DCL_FILE_NAME
+	#endif
+	#ifdef HAVE_OP_DZZ
+	  #include COEFF_DZZ_DCL_FILE_NAME
+	#endif
+
+
+	/*
+	 * ***fetch*** interpolation point coordinates
+	 */
+	fp interp_coords_fp[N_DIMS];
+	int axis;
+		for (axis = 0 ; axis < N_DIMS ; ++axis)
+		{
+		/* pointer to array of interp coords for this axis */
+		const void *const interp_coords_ptr = interp_coords[axis];
+
+		switch	(interp_coords_type_code)
+			{
+
+		case CCTK_VARIABLE_REAL:
+			  {
+			const CCTK_REAL *const interp_coords_ptr_real
+				= (const CCTK_REAL *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real[pt];
+			break;
+			  }
+
+		#ifdef HAVE_CCTK_REAL4
+		case CCTK_VARIABLE_REAL4:
+			  {
+			const CCTK_REAL4 *const interp_coords_ptr_real4
+				= (const CCTK_REAL4 *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real4[pt];
+			break;
+			  }
+		#endif
+
+		#ifdef HAVE_CCTK_REAL8
+		case CCTK_VARIABLE_REAL8:
+			  {
+			const CCTK_REAL8 *const interp_coords_ptr_real8
+				= (const CCTK_REAL8 *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real8[pt];
+			break;
+			  }
+		#endif
+
+		#ifdef HAVE_CCTK_REAL16
+		case CCTK_VARIABLE_REAL16:
+			  {
+			/* FIXME: maybe we should warn (once per cactus run) */
+			/*        that we may be doing arithmetic in lower */
+			/*        precision than the interp coords? */
+			const CCTK_REAL16 *const interp_coords_ptr_real16
+				= (const CCTK_REAL16 *) interp_coords_ptr;
+			interp_coords_fp[axis]
+				= interp_coords_ptr_real16[pt];
+			break;
+			  }
+		#endif
+
+		default:
+			CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "interp-coords datatype %d not supported!",
+				   interp_coords_type_code);
+								/*NOTREACHED*/
+			return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
+			/* end of switch (interp_coords_type_code) */
+			}
+
+		/* end of for (axis = ...) loop */
+		}
+
+
+	/*
+	 * locate interpolation molecules with respect to the grid, i.e.
+	 * compute (x,y,z) = coordinates of interpolation point
+	 * relative to molecule center, in units of the grid spacing
+	 *
+	 * n.b. we need the final answers in variables with the magic
+	 *      names (x,y,z) (machine-generated code uses these names),
+	 *      but we use temp variables as intermediates for (likely)
+	 *	better performance: the temp variables have their addresses
+	 *	taken and so may not be register-allocated, whereas the
+	 *	final (x,y,z) are "clean" and thus more likely to be
+	 *      register-allocated
+	 */
+	  {
+	#if (N_DIMS >= 1)
+	  fp x_temp;
+	  const int center_i
+		= LocalInterp_molecule_posn(origin_x, delta_x,
+					    input_array_min_subscripts[X_AXIS],
+					    input_array_max_subscripts[X_AXIS],
+					    out_of_range_tolerance[X_AXIS],
+					    MOLECULE_SIZE,
+					    interp_coords_fp[X_AXIS],
+					    &x_temp,
+					    (int *) NULL, (int *) NULL);
+	  const fp x = x_temp;
+	#endif
+	#if (N_DIMS >= 2)
+	  fp y_temp;
+	  const int center_j
+		= LocalInterp_molecule_posn(origin_y, delta_y,
+					    input_array_min_subscripts[Y_AXIS],
+					    input_array_max_subscripts[Y_AXIS],
+					    out_of_range_tolerance[Y_AXIS],
+					    MOLECULE_SIZE,
+					    interp_coords_fp[Y_AXIS],
+					    &y_temp,
+					    (int *) NULL, (int *) NULL);
+	  const fp y = y_temp;
+	#endif
+	#if (N_DIMS >= 3)
+	  fp z_temp;
+	  const int center_k
+		= LocalInterp_molecule_posn(origin_z, delta_z,
+					    input_array_min_subscripts[Z_AXIS],
+					    input_array_max_subscripts[Z_AXIS],
+					    out_of_range_tolerance[Z_AXIS],
+					    MOLECULE_SIZE,
+					    interp_coords_fp[Z_AXIS],
+					    &z_temp,
+					    (int *) NULL, (int *) NULL);
+	  const fp z = z_temp;
+	#endif
+
+	/* is the interpolation point out-of-range? */
+	#if (N_DIMS >= 1)
+	  if ((center_i == INT_MIN) || (center_i == INT_MAX))
+	     then {
+		  *error_pt   = pt;
+		  *error_axis = X_AXIS;
+		  *error_end  = (center_i == INT_MIN) ? -1 : +1;
+		  return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
+		  }
+	#endif
+	#if (N_DIMS >= 2)
+	  if ((center_j == INT_MIN) || (center_j == INT_MAX))
+	     then {
+		  *error_pt   = pt;
+		  *error_axis = Y_AXIS;
+		  *error_end  = (center_j == INT_MIN) ? -1 : +1;
+		  return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
+		  }
+	#endif
+	#if (N_DIMS >= 3)
+	  if ((center_k == INT_MIN) || (center_k == INT_MAX))
+	     then {
+		  *error_pt   = pt;
+		  *error_axis = Z_AXIS;
+		  *error_end  = (center_k == INT_MIN) ? -1 : +1;
+		  return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
+		  }
+	#endif
+
+	/* 1D subscript of molecule center in input arrays */
+	/* (this is used in DATA data-fetching macros below) */
+	  {
+	#if   (N_DIMS == 1)
+	  const int center_sub = SUB1(center_i);
+	#elif (N_DIMS == 2)
+	  const int center_sub = SUB2(center_i, center_j);
+	#elif (N_DIMS == 3)
+	  const int center_sub = SUB3(center_i, center_j, center_k);
+	#else
+	  #error "N_DIMS must be 1, 2, or 3!"
+	#endif
+
+
+	/*
+	 * compute the coefficients for whichever derivatives are wanted
+	 * using machine-generated coefficient expressions
+	 * (these are polynomials in the variables (x,y,z))
+	 */
+	#ifdef HAVE_OP_I
+	  if (want_I)
+	     then {
+		  #include COEFF_I_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DX
+	  if (want_dx)
+	     then {
+		  #include COEFF_DX_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DY
+	  if (want_dy)
+	     then {
+		  #include COEFF_DY_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DZ
+	  if (want_dz)
+	     then {
+		  #include COEFF_DZ_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DXX
+	  if (want_dxx)
+	     then {
+		  #include COEFF_DXX_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DXY
+	  if (want_dxy)
+	     then {
+		  #include COEFF_DXY_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DXZ
+	  if (want_dxz)
+	     then {
+		  #include COEFF_DXZ_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DYY
+	  if (want_dyy)
+	     then {
+		  #include COEFF_DYY_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DYZ
+	  if (want_dyz)
+	     then {
+		  #include COEFF_DYZ_COMPUTE_FILE_NAME
+		  }
+	#endif
+	#ifdef HAVE_OP_DZZ
+	  if (want_dzz)
+	     then {
+		  #include COEFF_DZZ_COMPUTE_FILE_NAME
+		  }
+	#endif
+
+
+	/*
+	 * compute each output array at this point
+	 */
+	  {
+	int out;
+	const void *input_array_ptr = NULL;
+		for (out = 0 ; out < N_output_arrays ; ++out)
+		{
+		const int in = operand_indices[out];
+		const int input_offset = input_array_offsets[in];
+		const int sub_temp = input_offset + center_sub;
+
+		/*
+ 		 * ***decode*** the input/output array datatypes
+		 * to determine whether they're real or complex,
+		 * and verify that they're both the same in this regard
+		 * ==> define
+		 *	const int N_parts = data is complex ? 2 : 1;
+		 */
+		const int N_input_parts
+		  = LocalInterp_decode_N_parts(input_array_type_codes[in]);
+		const int N_output_parts
+		  = LocalInterp_decode_N_parts(output_array_type_codes[out]);
+		if (N_input_parts != N_output_parts)
+		   then {
+CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "\n"
+	   "   can't do real input --> complex output\n"
+	   "   or complex input --> real output interpolation!\n"
+	   "   (0-origin) input #%d datatype = %d, output #%d datatype = %d\n"
+	   ,
+	   in, (int) input_array_type_codes[in],
+	   out, (int) output_array_type_codes[out]);	/*NOTREACHED*/
+			return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
+			}
+
+		  {
+		const int N_parts = N_input_parts;
+		int part;
+			for (part = 0 ; part < N_parts ; ++part)
+			{
+			if ( (input_arrays[in] != input_array_ptr)
+			     || (part != 0) )
+			   then {
+				/*
+				 * ***fetch*** the input array values
+				 *             into local variables
+				 */
+				input_array_ptr = input_arrays[in];
+				switch	(input_array_type_codes[in])
+					{
+
+case CCTK_VARIABLE_REAL:
+	  {
+	const CCTK_REAL *const input_array_ptr_real
+		= (const CCTK_REAL *) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_real[sub_temp + SUB1(i)]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_real[sub_temp + SUB2(i,j)]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_real[sub_temp + SUB3(i,j,k)]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+
+#ifdef HAVE_CCTK_REAL4
+case CCTK_VARIABLE_REAL4:
+	  {
+	const CCTK_REAL4 *const input_array_ptr_real4
+		= (const CCTK_REAL4 *) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_real4[sub_temp + SUB1(i)]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_real4[sub_temp + SUB2(i,j)]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_real4[sub_temp + SUB3(i,j,k)]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_REAL4 */
+
+#ifdef HAVE_CCTK_REAL8
+case CCTK_VARIABLE_REAL8:
+	  {
+	const CCTK_REAL8 *const input_array_ptr_real8
+		= (const CCTK_REAL8 *) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_real8[sub_temp + SUB1(i)]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_real8[sub_temp + SUB2(i,j)]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_real8[sub_temp + SUB3(i,j,k)]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_REAL8 */
+
+#ifdef HAVE_CCTK_REAL16
+case CCTK_VARIABLE_REAL16:
+	  {
+	/* FIXME: maybe we should warn (once per cactus run) that we may be */
+	/*        doing arithmetic in lower precision than the data type? */
+	const CCTK_REAL16 *const input_array_ptr_real16
+		= (const CCTK_REAL16 *) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_real16[sub_temp + SUB1(i)]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_real16[sub_temp + SUB2(i,j)]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_real16[sub_temp + SUB3(i,j,k)]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_REAL16 */
+
+case CCTK_VARIABLE_COMPLEX:
+	  {
+	const CCTK_REAL (*const input_array_ptr_complex)[COMPLEX_N_PARTS]
+		= (const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_complex[sub_temp + SUB1(i)][part]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_complex[sub_temp + SUB2(i,j)][part]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_complex[sub_temp + SUB3(i,j,k)][part]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+
+#ifdef HAVE_CCTK_COMPLEX8
+case CCTK_VARIABLE_COMPLEX8:
+	  {
+	const CCTK_REAL4 (*const input_array_ptr_complex8)[COMPLEX_N_PARTS]
+		= (const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_complex8[sub_temp + SUB1(i)][part]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_complex8[sub_temp + SUB2(i,j)][part]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_complex8[sub_temp + SUB3(i,j,k)][part]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX8 */
+
+#ifdef HAVE_CCTK_COMPLEX16
+case CCTK_VARIABLE_COMPLEX16:
+	  {
+	const CCTK_REAL8 (*const input_array_ptr_complex16)[COMPLEX_N_PARTS]
+		= (const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_complex16[sub_temp + SUB1(i)][part]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_complex16[sub_temp + SUB2(i,j)][part]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_complex16[sub_temp + SUB3(i,j,k)][part]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX16 */
+
+#ifdef HAVE_CCTK_COMPLEX32
+case CCTK_VARIABLE_COMPLEX32:
+	  {
+	const CCTK_REAL16 (*const input_array_ptr_complex32)[COMPLEX_N_PARTS]
+		= (const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr;
+	#undef DATA
+  #if   (N_DIMS == 1)
+    #define DATA(i)	input_array_ptr_complex32[sub_temp + SUB1(i)][part]
+  #elif (N_DIMS == 2)
+    #define DATA(i,j)	input_array_ptr_complex32[sub_temp + SUB2(i,j)][part]
+  #elif (N_DIMS == 3)
+    #define DATA(i,j,k)	input_array_ptr_complex32[sub_temp + SUB3(i,j,k)][part]
+  #else
+    #error "N_DIMS must be 1, 2, or 3!"
+  #endif
+	#include DATA_VAR_ASSIGN_FILE_NAME
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX32 */
+
+default:
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "input datatype %d not supported!",
+		   input_array_type_codes[in]);		/*NOTREACHED*/
+	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+					}
+				}
+
+			/*
+			 * compute the interpolant itself
+			 * as a linear combination of the data variables
+			 */
+			  {
+			fp result;
+			switch	(operation_codes[out])
+				{
+			#ifdef HAVE_OP_I
+			  case 0:
+				#include INTERP_I_COMPUTE_FILE_NAME
+				break;
+			#endif
+			#ifdef HAVE_OP_DX
+			  case 1:
+				#include INTERP_DX_COMPUTE_FILE_NAME
+				result *= inverse_delta_x;
+				break;
+			#endif
+			#ifdef HAVE_OP_DY
+			  case 2:
+				#include INTERP_DY_COMPUTE_FILE_NAME
+				result *= inverse_delta_y;
+				break;
+			#endif
+			#ifdef HAVE_OP_DZ
+			  case 3:
+				#include INTERP_DZ_COMPUTE_FILE_NAME
+				result *= inverse_delta_z;
+				break;
+			#endif
+			#ifdef HAVE_OP_DXX
+			  case 11:
+				#include INTERP_DXX_COMPUTE_FILE_NAME
+				result *= inverse_delta_x2;
+				break;
+			#endif
+			#ifdef HAVE_OP_DXY
+			  case 12:
+			  case 21:
+				#include INTERP_DXY_COMPUTE_FILE_NAME
+				result *= inverse_delta_x * inverse_delta_y;
+				break;
+			#endif
+			#ifdef HAVE_OP_DXZ
+			  case 13:
+			  case 31:
+				#include INTERP_DXZ_COMPUTE_FILE_NAME
+				result *= inverse_delta_x * inverse_delta_z;
+				break;
+			#endif
+			#ifdef HAVE_OP_DYY
+			  case 22:
+				#include INTERP_DYY_COMPUTE_FILE_NAME
+				result *= inverse_delta_y2;
+				break;
+			#endif
+			#ifdef HAVE_OP_DYZ
+			  case 23:
+			  case 32:
+				#include INTERP_DYZ_COMPUTE_FILE_NAME
+				result *= inverse_delta_y * inverse_delta_z;
+				break;
+			#endif
+			#ifdef HAVE_OP_DZZ
+			  case 33:
+				#include INTERP_DZZ_COMPUTE_FILE_NAME
+				result *= inverse_delta_z2;
+				break;
+			#endif
+			default:
+CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "Generalized interpolation operation_code %d not supported!",
+	   operation_codes[out]);				/*NOTREACHED*/
+				return UTIL_ERROR_BAD_INPUT;
+							/*** ERROR RETURN ***/
+				/* end of switch (operation_codes[out]) */
+				}
+
+
+			/*
+			 * ***store*** the result in the output array
+			 */
+			switch	(output_array_type_codes[out])
+				{
+
+case CCTK_VARIABLE_REAL:
+	  {
+	CCTK_REAL *const output_array_ptr_real
+		= (CCTK_REAL *) output_arrays[out];
+	output_array_ptr_real[pt] = result;
+	break;
+	  }
+
+#ifdef HAVE_CCTK_REAL4
+case CCTK_VARIABLE_REAL4:
+	  {
+	CCTK_REAL4 *const output_array_ptr_real4
+		= (CCTK_REAL4 *) output_arrays[out];
+	output_array_ptr_real4[pt] = result;
+	break;
+	  }
+#endif
+
+#ifdef HAVE_CCTK_REAL8
+case CCTK_VARIABLE_REAL8:
+	  {
+	CCTK_REAL8 *const output_array_ptr_real8
+		= (CCTK_REAL8 *) output_arrays[out];
+	output_array_ptr_real8[pt] = result;
+	break;
+	  }
+#endif
+
+#ifdef HAVE_CCTK_REAL16
+case CCTK_VARIABLE_REAL16:
+	  {
+	CCTK_REAL16 *const output_array_ptr_real16
+		= (CCTK_REAL16 *) output_arrays[out];
+	output_array_ptr_real16[pt] = result;
+	break;
+	  }
+#endif
+
+case CCTK_VARIABLE_COMPLEX:
+	  {
+	CCTK_COMPLEX *const output_array_ptr_complex
+		= (CCTK_COMPLEX *) output_arrays[out];
+	((CCTK_REAL *)  & output_array_ptr_complex[pt]) [part]
+		= result;
+	break;
+	  }
+
+#ifdef HAVE_CCTK_COMPLEX8
+case CCTK_VARIABLE_COMPLEX8:
+	  {
+	CCTK_COMPLEX8 *const output_array_ptr_complex8
+		= (CCTK_COMPLEX8 *) output_arrays[out];
+	((CCTK_REAL4 *)  & output_array_ptr_complex8[pt]) [part]
+		= result;
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX8 */
+
+#ifdef HAVE_CCTK_COMPLEX16
+case CCTK_VARIABLE_COMPLEX16:
+	  {
+	CCTK_COMPLEX16 *const output_array_ptr_complex16
+		= (CCTK_COMPLEX16 *) output_arrays[out];
+	((CCTK_REAL8 *)  & output_array_ptr_complex16[pt]) [part]
+		= result;
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX16 */
+
+#ifdef HAVE_CCTK_COMPLEX32
+case CCTK_VARIABLE_COMPLEX32:
+	  {
+	CCTK_COMPLEX32 *const output_array_ptr_complex32
+		= (CCTK_COMPLEX32 *) output_arrays[out];
+	((CCTK_REAL16 *)  & output_array_ptr_complex32[pt]) [part]
+		= result;
+	break;
+	  }
+#endif	/* HAVE_CCTK_COMPLEX32 */
+
+default:
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "output datatype %d not supported",
+		   output_array_type_codes[out]);		/*NOTREACHED*/
+	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+				/* end of switch (output type code) */
+				}
+			  }
+			/* end of  for (part = ...)  loop */
+			}
+		  }
+		/* end of  for (out = ...)  loop */
+		}
+	  }
+	  }
+	  }
+
+	/* end of  for (pt = ...)  loop */
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+  }
+  }
+}