diff options
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/template.c')
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/template.c | 2033 |
1 files changed, 1098 insertions, 935 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c index 5cb934b..3a09a34 100644 --- a/src/GeneralizedPolynomial-Uniform/template.c +++ b/src/GeneralizedPolynomial-Uniform/template.c @@ -1,15 +1,26 @@ /*@@ @file template.c - @date 23 Oct 2001 + @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). + such inclusion defines a single interpolation function + (plus other supporting functions local to this file), + which does generalized interpolation for a single molecule, + or more accurately for a single molecule for each possible + operation_codes[] value. + @enddesc + + @hdate 29.Aug.2002 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Restructure code and break it up into smaller functions + to reduce excessive compilation cpu/memory with big + 3-D molecules. + @endhdesc + @desc The following header files must be #included before #including this file: <math.h> @@ -29,6 +40,7 @@ "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 @@ -44,6 +56,18 @@ sizing arrays etc. @endvar + @var XYZ + @desc A comma-separated list of the (x,y,z) coordinates in the + interpolation, e.g. + #define XYZ x,y + @endvar + + @var FP_XYZ + @desc A comma-separated list of function-prototype declarations + of the (x,y,z) coordinates in the interpolation, e.g. + #define FP_XYZ fp x, fp y + @endvar + @var MOLECULE_MIN_M @vdesc The minimum m coordinate in the molecule, e.g. #define MOLECULE_MIN_M -1 @@ -87,37 +111,37 @@ fp data_p1_p1; @endvar - @var COEFF_{I,DX,DY,DXX,DXY,DYY,...}_DCL_FILE_NAME + @var COEFFS_{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; + fp coeffs_dx_m1_m1; + fp coeffs_dx_0_m1; + fp coeffs_dx_p1_m1; + fp coeffs_dx_m1_0; + fp coeffs_dx_0_0; + fp coeffs_dx_p1_0; + fp coeffs_dx_m1_p1; + fp coeffs_dx_0_p1; + fp coeffs_dx_p1_p1; @endvar - @var COEFF_{I,DX,DY,DXX,DXY,DYY,...}_VAR_STORE_FILE_NAME + @var COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_VAR_STORE_FILE_NAME @vdesc The name of a file (presumably machine-generated) containing a sequence of C assignment statements to assign the interpolation coefficients to the corresponding COEFF(...) expressions, eg (for 2D size-3 molecules) - COEFF(-1,-1) = coeff_dx_m1_m1; - COEFF(0,-1) = coeff_dx_0_m1; - COEFF(1,-1) = coeff_dx_p1_m1; - COEFF(-1,0) = coeff_dx_m1_0; - COEFF(0,0) = coeff_dx_0_0; - COEFF(1,0) = coeff_dx_p1_0; - COEFF(-1,1) = coeff_dx_m1_p1; - COEFF(0,1) = coeff_dx_0_p1; - COEFF(1,1) = coeff_dx_p1_p1; + COEFF(-1,-1) = coeffs_dx_m1_m1; + COEFF(0,-1) = coeffs_dx_0_m1; + COEFF(1,-1) = coeffs_dx_p1_m1; + COEFF(-1,0) = coeffs_dx_m1_0; + COEFF(0,0) = coeffs_dx_0_0; + COEFF(1,0) = coeffs_dx_p1_0; + COEFF(-1,1) = coeffs_dx_m1_p1; + COEFF(0,1) = coeffs_dx_0_p1; + COEFF(1,1) = coeffs_dx_p1_p1; @endvar @var DATA_VAR_ASSIGN_FILE_NAME @@ -143,31 +167,31 @@ 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), + (cf COEFFS_{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; + coeffs_dx_m1_m1*data_m1_m1 + + coeffs_dx_0_m1*data_0_m1 + + coeffs_dx_p1_m1*data_p1_m1 + + coeffs_dx_p2_m1*data_p2_m1 + + coeffs_dx_m1_0*data_m1_0 + + coeffs_dx_0_0*data_0_0 + + coeffs_dx_p1_0*data_p1_0 + + coeffs_dx_p2_0*data_p2_0 + + coeffs_dx_m1_p1*data_m1_p1 + + coeffs_dx_0_p1*data_0_p1 + + coeffs_dx_p1_p1*data_p1_p1 + + coeffs_dx_p2_p1*data_p2_p1 + + coeffs_dx_m1_p2*data_m1_p2 + + coeffs_dx_0_p2*data_0_p2 + + coeffs_dx_p1_p2*data_p1_p2 + + coeffs_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 + @var COEFFS_{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 @@ -186,15 +210,15 @@ 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; + coeffs_dx_m1_m1 = t40+t41; + coeffs_dx_0_m1 = t36; + coeffs_dx_p1_m1 = t39+t42; + coeffs_dx_m1_0 = t40+t35; + coeffs_dx_0_0 = t36; + coeffs_dx_p1_0 = t35+t39; + coeffs_dx_m1_p1 = t40+t42; + coeffs_dx_0_p1 = t36; + coeffs_dx_p1_p1 = t39+t41; As illustrated, the code may use the macro RATIONAL (defined later in this file) to represent rational-number @@ -208,13 +232,60 @@ @endvar @version $Header$ -@@*/ + @@*/ + +/******************************************************************************/ + +/* + * ***** prototypes for functions local to this file ***** + */ + +#ifdef HAVE_OP_I + static + void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I); +#endif +#ifdef HAVE_OP_DX + static + void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx); +#endif +#ifdef HAVE_OP_DY + static + void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy); +#endif +#ifdef HAVE_OP_DZ + static + void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz); +#endif +#ifdef HAVE_OP_DXX + static + void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx); +#endif +#ifdef HAVE_OP_DXY + static + void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy); +#endif +#ifdef HAVE_OP_DXZ + static + void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz); +#endif +#ifdef HAVE_OP_DYY + static + void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy); +#endif +#ifdef HAVE_OP_DYZ + static + void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz); +#endif +#ifdef HAVE_OP_DZZ + static + void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz); +#endif /******************************************************************************/ /*@@ @routine FUNCTION_NAME - @date 23 Oct 2001 + @date 23.Oct.2001 @author Jonathan Thornburg <jthorn@aei.mpg.de> @desc This function does generalized interpolation of one or more @@ -305,114 +376,168 @@ int FUNCTION_NAME(/***** coordinate system *****/ struct Jacobian_info* Jacobian_info) { /* - * Implementation notes: + * ***** Naming conventions ***** + * + * in, out = 0-origin indices each selecting an input/output array + * pt = 0-origin index selecting an interpolation point + * part = 0-origin index selecting real/imaginary part of a complex number + */ + +/* + * ***** Implementation notes: ***** * * The basic outline of this function is as follows: * - * store molecule structure flags - * if (querying molecule min/max m) - * then store molecule min/max m - * compute "which derivatives are wanted" flags - * precompute 1/dx factors - * for (int pt = 0 ; pt < N_interp_point ; ++pt) + * store molecule structure flags + * if (querying molecule min/max m) + * then store molecule min/max m + * compute "which derivatives are wanted" flags + * precompute 1/dx factors + * set up input array strides + * set up Jacobian strides (or dummy values if no Jacobian info) + * for (int pt = 0 ; pt < N_interp_point ; ++pt) + * { + * struct DATA_STRUCT data; + * struct COEFFS_STRUCT coeffs_I, coeffs_dx, ..., coeffs_dzz; + * for (int axis = 0 ; axis < N_dims ; ++axis) + * { + * switch (interp_coords_type_codes) + * { + * case CCTK_VARIABLE_REAL: + * fetch interp coords for this axis from + * const CCTK_REAL *const interp_coords[] + * break; + * ... + * } + * } + * compute position of interpolation molecules + * with respect to the grid into the XYZ variables + * if (querying molecule positions) + * then store this molecule position + * compute molecule center 1-d position in input arrays + * + * /# compute interpolation coefficients at this point #/ + * if (want_I) + * then compute_coeffs_I(XYZ, &coeffs_I); + * if (want_dx) + * then compute_coeffs_dx(XYZ, &coeffs_dx); + * ... + * if (want_dzz) + * then compute_coeffs_dzz(XYZ, &coeffs_dx); + * + * for (int out = 0 ; out < N_output_arrays ; ++out) + * { + * const int in = operand_indices[out]; + * ***decode*** the input/output array datatypes + * to determine whether they're real or complex + * (they must both be the same in this regard), then + * const int N_parts = data is complex ? 2 : 1; + * for (int part = 0 ; part < N_parts ; ++part) * { - * declare all the coefficients - * declare all the data-values variables - * ***fetch*** interpolation point coordinates - * compute position of interpolation molecules - * with respect to the grid - * if (querying molecule positions) - * then store this molecule position - * compute coefficients for all derivatives which are wanted - * for (int out = 0 ; out < N_output_arrays ; ++out) - * { - * const int in = operand_indices[out]; - * ***decode*** the input/output array datatypes - * to determine whether they're real or complex - * (they must both be the same in this regard), then - * const int N_parts = data is complex ? 2 : 1; - * for (int part = 0 ; part < N_parts ; ++part) + * if ( (input_arrays[in] != NULL) + * && ( (input_arrays[in] != value at last fetch) + * || (part != value at last fetch) ) ) + * then { + * save input_arrays[in] and part for + * "previous value" test above + * switch (input_array_type_codes[in]) + * { + * case CCTK_VARIABLE_REAL: + * FETCH_DATA_REAL(input_array_ptr, + * STRIDE_IJK, + * &data); + * break; + * ... + * case CCTK_VARIABLE_COMPLEX: + * FETCH_DATA_COMPLEX(input_array_ptr, + * STRIDE_IJK, part, + * &data); + * break; + * ... + * } + * } + * + * if (output_arrays[out] != NULL) + * then { + * /# interpolate at this point #/ + * fp result; + * switch (operation_codes[out]) + * { + * case 0: + * result = EVALUATE_MOLECULE(&coeffs_I, + * &data); + * break; + * case 1: + * result = inverse_dx + * * EVALUATE_MOLECULE(&coeffs_dx, + * &data); + * break; + * ... + * case 33: + * result = inverse_dz * inverse_dz + * * EVALUATE_MOLECULE(&coeffs_dzz, + * &data); + * break; + * } + * + * /# store result in output array #/ + * switch (output_array_type_codes[out]) * { - * if ( (input_arrays[in] != NULL) - * && ( (input_arrays[in] != value at - * last fetch) - * || (part != value at last fetch) ) ) - * then { - * save input_arrays[in] and part for - * "previous value" test above - * ***fetch*** the molecule-sized piece - * of input_arrays[in][part] - * at this molecule position - * into local "data" variables - * } - * if (output_arrays[out] != NULL) - * then { - * fp result; - * switch (operation_code) - * { - * case 0: - * result = compute interpolant - * as a linear combination - * of data variables - * and op=0 coeffs - * break; - * case 1: - * result = compute interpolant - * as a linear combination - * of data variables - * and op=1 coeffs - * result *= inverse_dx; - * break; - * case ... - * } - * ***store*** result in output array - * } - * if (querying Jacobian - * && (Jacobian_pointer[out] != NULL)) - * then { - * fp factor; - * switch (operation_code) - * { - * case 0: - * store op=0 Jacobian values - * break; - * case 1: - * factor = inverse_dx; - * store factor - * * op=1 Jacobian values - * break; - * case ... - * } - * } + * case CCTK_VARIABLE_REAL: + * store result in + * CCTK_REAL *const output_arrays[] + * break; + * ... + * case CCTK_VARIABLE_COMPLEX: + * store result in + * CCTK_COMPLEX *const output_arrays[][part] + * break; + * ... + * } + * } + * if (querying Jacobian && (Jacobian_pointer[out] != NULL)) + * then { + * /# store Jacobian coefficients at this point #/ + * CCTK_REAL *const Jacobian_ptr + * = Jacobian_info->Jacobian_pointer[out] + * + Jacobian_info->Jacobian_offset[out]; + * switch (operation_code) + * { + * case 0: + * STORE_COEFFS(1.0, &coeffs_I, + * Jacobian_ptr, JACOBIAN_STRIDE_IJK, + * pt, part, + * break; + * case 1: + * STORE_COEFFS(inverse_dx, &coeffs_dx, + * Jacobian_ptr, JACOBIAN_STRIDE_IJK, + * pt, part, + * break; + * ... + * case 33: + * STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, + * Jacobian_ptr, JACOBIAN_STRIDE_IJK, + * pt, part, + * break; * } * } * } + * } + * } * - * Here "***decode***", "***fetch***", and "***store***" are all - * actually switches on the various array datatypes. For complex - * datatypes "***fetch***" and "***store***" pointer-alias 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. + * For complex datatypes we pointer-alias 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 - * part = 0-origin index selecting real/imaginary part of a complex number + * 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 + * interpolation errors. */ -/* 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 @@ -439,28 +564,6 @@ int FUNCTION_NAME(/***** coordinate system *****/ #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); @@ -472,7 +575,7 @@ int FUNCTION_NAME(/***** coordinate system *****/ /* - * Jacobian structure info + * store molecule structure flags, molecule min/max m (if requested) */ if (molecule_structure_flags != NULL) then { @@ -493,7 +596,7 @@ if (molecule_min_max_m_info != NULL) /* - * compute flags specifying which derivatives are wanted + * compute "which derivatives are wanted" flags */ { #ifdef HAVE_OP_I @@ -590,337 +693,378 @@ int out; #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); + #if defined(HAVE_OP_DX) \ + || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \ + || defined(HAVE_OP_DXX) + const fp inverse_dx = 1.0 / 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); + #if defined(HAVE_OP_DY) \ + || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \ + || defined(HAVE_OP_DYY) + const fp inverse_dy = 1.0 / 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); + #if defined(HAVE_OP_DZ) \ + || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \ + || defined(HAVE_OP_DZZ) + const fp inverse_dz = 1.0 / delta_z; #endif #endif +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif /* + * set up input array strides + */ +#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 > 3) + #error "N_DIMS may not be > 3!" +#endif + + +/* + * set up Jacobian m strides, or dummy values if no Jacobian info + */ +#if N_DIMS >= 1 + const int Jacobian_mi_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[X_AXIS]; +#endif +#if (N_DIMS >= 2) + const int Jacobian_mj_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[Y_AXIS]; +#endif +#if (N_DIMS >= 3) + const int Jacobian_mk_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[Z_AXIS]; +#endif +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + +/* * interpolate at each point */ int pt; - for (pt = 0 ; pt < N_interp_points ; ++pt) + for (pt = 0 ; pt < N_interp_points ; ++pt) + { + /* this struct holds a molecule-sized piece of a single */ + /* real input array, or of a real/complex component of */ + /* a complex input array */ + struct DATA_STRUCT data; + + /* each of these structs holds a molecule-sized set of */ + /* interpolation coefficients -- recall we are representing */ + /* the interpolant as a linear combination of the data values */ + #ifdef HAVE_OP_I + struct COEFFS_STRUCT coeffs_I; + #endif + #ifdef HAVE_OP_DX + struct COEFFS_STRUCT coeffs_dx; + #endif + #ifdef HAVE_OP_DY + struct COEFFS_STRUCT coeffs_dy; + #endif + #ifdef HAVE_OP_DZ + struct COEFFS_STRUCT coeffs_dz; + #endif + #ifdef HAVE_OP_DXX + struct COEFFS_STRUCT coeffs_dxx; + #endif + #ifdef HAVE_OP_DXY + struct COEFFS_STRUCT coeffs_dxy; + #endif + #ifdef HAVE_OP_DXZ + struct COEFFS_STRUCT coeffs_dxz; + #endif + #ifdef HAVE_OP_DYY + struct COEFFS_STRUCT coeffs_dyy; + #endif + #ifdef HAVE_OP_DYZ + struct COEFFS_STRUCT coeffs_dyz; + #endif + #ifdef HAVE_OP_DZZ + struct COEFFS_STRUCT coeffs_dzz; + #endif + + + /* + * fetch the interpolation point coordinates + * from the interp_coords[] arrays + * FIXME: Maybe it would be better (faster) to do this + * with N_DIMS open-coded calls on a function? + * But then we'd have to have a sentinal value + * return for the unknown-type-code error case. + * Yuk! :( :( + */ + fp interp_coords_fp[N_DIMS]; + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) { - /* declare all the data-values variables */ - #include DATA_VAR_DCL_FILE_NAME + /* pointer to array of interp coords for this axis */ + const void *const interp_coords_ptr = interp_coords[axis]; - /* 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) + switch (interp_coords_type_code) { - /* 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, -"\n" -" CCTK_InterpLocalUniform(): interpolation coordinates datatype %d\n" -" is 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 */ - } - - - /* - * compute position of 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], - MOLECULE_SIZE, - out_of_range_tolerance[X_AXIS_MIN], - out_of_range_tolerance[X_AXIS_MAX], - 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], - MOLECULE_SIZE, - out_of_range_tolerance[Y_AXIS_MIN], - out_of_range_tolerance[Y_AXIS_MAX], - 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], - MOLECULE_SIZE, - out_of_range_tolerance[Z_AXIS_MIN], - out_of_range_tolerance[Z_AXIS_MAX], - 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 { - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = X_AXIS; - error_flags->error_end = (center_i > 0) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ + 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; } - #endif - #if (N_DIMS >= 2) - if ((center_j == INT_MIN) || (center_j == INT_MAX)) - then { - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Y_AXIS; - error_flags->error_end = (center_j > 0) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ + + #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 - #if (N_DIMS >= 3) - if ((center_k == INT_MIN) || (center_k == INT_MAX)) - then { - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Z_AXIS; - error_flags->error_end = (center_k > 0) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ + + #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 - /* 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!" + #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 - /* - * molecule position queries - */ - if (molecule_positions != NULL) - then { - #if (N_DIMS >= 1) - molecule_positions[X_AXIS][pt] = center_i; - #endif - #if (N_DIMS >= 2) - molecule_positions[Y_AXIS][pt] = center_j; - #endif - #if (N_DIMS >= 3) - molecule_positions[Z_AXIS][pt] = center_k; - #endif + default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): interpolation coordinates datatype %d\n" +" is 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 */ + } - /* - * 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 position of 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], + MOLECULE_SIZE, + out_of_range_tolerance[X_AXIS_MIN], + out_of_range_tolerance[X_AXIS_MAX], + 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], + MOLECULE_SIZE, + out_of_range_tolerance[Y_AXIS_MIN], + out_of_range_tolerance[Y_AXIS_MAX], + 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], + MOLECULE_SIZE, + out_of_range_tolerance[Z_AXIS_MIN], + out_of_range_tolerance[Z_AXIS_MAX], + interp_coords_fp[Z_AXIS], + &z_temp, + (int *) NULL, (int *) NULL); + const fp z = z_temp; + #endif + #if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" + #endif + + /* is the interpolation point out-of-range? */ + #if (N_DIMS >= 1) + if ((center_i == INT_MIN) || (center_i == INT_MAX)) + then { + if (error_flags != NULL) + then { + error_flags->error_pt = pt; + error_flags->error_axis = X_AXIS; + error_flags->error_end = (center_i > 0) ? +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 { + if (error_flags != NULL) + then { + error_flags->error_pt = pt; + error_flags->error_axis = Y_AXIS; + error_flags->error_end = (center_j > 0) ? +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 { + if (error_flags != NULL) + then { + error_flags->error_pt = pt; + error_flags->error_axis = Z_AXIS; + error_flags->error_end = (center_k > 0) ? +1 : -1; + } + return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ + } + #endif + + + /* + * compute 1-d position of molecule center in input arrays + */ + { + #if (N_DIMS == 1) + const int molecule_center_posn = stride_i*center_i; + #endif + #if (N_DIMS == 2) + const int molecule_center_posn = stride_i*center_i + + stride_j*center_j; + #endif + #if (N_DIMS == 3) + const int molecule_center_posn = stride_i*center_i + + stride_j*center_j + + stride_k*center_k; + #endif + #if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" + #endif + + + /* + * molecule position queries + */ + if (molecule_positions != NULL) + then { + #if (N_DIMS >= 1) + molecule_positions[X_AXIS][pt] = center_i; + #endif + #if (N_DIMS >= 2) + molecule_positions[Y_AXIS][pt] = center_j; + #endif + #if (N_DIMS >= 3) + molecule_positions[Z_AXIS][pt] = center_k; + #endif + } + + + /* + * compute the coefficients at this point for whichever + * operation_codes[] values (derivatatives) are wanted + * (these are polynomials in the variables (x,y,z)) + */ + #ifdef HAVE_OP_I + if (want_I) + then compute_coeffs_I(XYZ, &coeffs_I); + #endif + #ifdef HAVE_OP_DX + if (want_dx) + then compute_coeffs_dx(XYZ, &coeffs_dx); + #endif + #ifdef HAVE_OP_DY + if (want_dy) + then compute_coeffs_dy(XYZ, &coeffs_dy); + #endif + #ifdef HAVE_OP_DZ + if (want_dz) + then compute_coeffs_dz(XYZ, &coeffs_dz); + #endif + #ifdef HAVE_OP_DXX + if (want_dxx) + then compute_coeffs_dxx(XYZ, &coeffs_dxx); + #endif + #ifdef HAVE_OP_DXY + if (want_dxy) + then compute_coeffs_dxy(XYZ, &coeffs_dxy); + #endif + #ifdef HAVE_OP_DXZ + if (want_dxz) + then compute_coeffs_dxz(XYZ, &coeffs_dxz); + #endif + #ifdef HAVE_OP_DYY + if (want_dyy) + then compute_coeffs_dyy(XYZ, &coeffs_dyy); + #endif + #ifdef HAVE_OP_DYZ + if (want_dyz) + then compute_coeffs_dyz(XYZ, &coeffs_dyz); + #endif + #ifdef HAVE_OP_DZZ + if (want_dzz) + then compute_coeffs_dzz(XYZ, &coeffs_dzz); + #endif /* @@ -937,25 +1081,23 @@ int pt; const void* input_array_ptr__last_fetch = NULL; int part__last_fetch = -1; - for (out = 0 ; out < N_output_arrays ; ++out) - { - const int in = operand_indices[out]; - const int input_center_sub - = input_array_offsets[in] + 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 { + for (out = 0 ; out < N_output_arrays ; ++out) + { + const int in = operand_indices[out]; + + /* + * ***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" " CCTK_InterpLocalUniform(): can't do real input --> complex output or\n" @@ -965,51 +1107,41 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, , 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; - const void* const input_array_ptr = input_arrays[in]; - int part; - for (part = 0 ; part < N_parts ; ++part) - { - if ( (input_array_ptr != NULL) - && - ( (input_array_ptr != input_array_ptr__last_fetch) - || (part != part__last_fetch) ) ) - then { - /* remember when we did the following fetch */ - input_array_ptr__last_fetch = input_array_ptr; - part__last_fetch = part; - - /* - * ***fetch*** the molecule-sized piece - * of input_arrays[in][part] - * at this molecule position - * into local "data" variables - */ - switch (input_array_type_codes[in]) - { - + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + + { + const int N_parts = N_input_parts; + const void* const input_array_ptr = input_arrays[in]; + int part; + for (part = 0 ; part < N_parts ; ++part) + { + if ( (input_array_ptr != NULL) + && + ((input_array_ptr != input_array_ptr__last_fetch) + || (part != part__last_fetch)) ) + then { + /* remember when we did the following fetch */ + input_array_ptr__last_fetch = input_array_ptr; + part__last_fetch = part; + + /* + * fetch the molecule-sized piece of + * input_arrays[in][part] at this molecule + * position, into the data struct + */ + { + const int input_posn = molecule_center_posn + + input_array_offsets[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(mi) \ - input_array_ptr_real[input_center_sub + SUB1(mi)] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_real[input_center_sub + SUB2(mi,mj)] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_real[input_center_sub + SUB3(mi,mj,mk)] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL *) input_array_ptr) + input_posn; + FETCH_DATA_REAL(input_array_ptr_real, + STRIDE_IJK, + &data); break; } @@ -1017,47 +1149,25 @@ case CCTK_VARIABLE_REAL: 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(mi) \ - input_array_ptr_real4[input_center_sub + SUB1(mi)] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_real4[input_center_sub + SUB2(mi,mj)] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_real4[input_center_sub + SUB3(mi,mj,mk)] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL4 *) input_array_ptr) + input_posn; + FETCH_DATA_REAL4(input_array_ptr_real4, + STRIDE_IJK, + &data); break; } -#endif /* HAVE_CCTK_REAL4 */ +#endif #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(mi) \ - input_array_ptr_real8[input_center_sub + SUB1(mi)] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_real8[input_center_sub + SUB2(mi,mj)] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_real8[input_center_sub + SUB3(mi,mj,mk)] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL8 *) input_array_ptr) + input_posn; + FETCH_DATA_REAL8(input_array_ptr_real8, + STRIDE_IJK, + &data); break; } -#endif /* HAVE_CCTK_REAL8 */ +#endif #ifdef HAVE_CCTK_REAL16 case CCTK_VARIABLE_REAL16: @@ -1065,43 +1175,22 @@ 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(mi) \ - input_array_ptr_real16[input_center_sub + SUB1(mi)] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_real16[input_center_sub + SUB2(mi,mj)] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_real16[input_center_sub + SUB3(mi,mj,mk)] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL16 *) input_array_ptr) + input_posn; + FETCH_DATA_REAL16(input_array_ptr_real16, + STRIDE_IJK, + &data); break; } -#endif /* HAVE_CCTK_REAL16 */ +#endif 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(mi) \ - input_array_ptr_complex[input_center_sub + SUB1(mi)][part] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_complex[input_center_sub + SUB2(mi,mj)][part] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_complex[input_center_sub + SUB3(mi,mj,mk)][part] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + FETCH_DATA_COMPLEX(input_array_ptr_complex, + STRIDE_IJK, part, + &data); break; } @@ -1109,400 +1198,474 @@ case CCTK_VARIABLE_COMPLEX: 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(mi) \ - input_array_ptr_complex8[input_center_sub + SUB1(mi)][part] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_complex8[input_center_sub + SUB2(mi,mj)][part] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_complex8[input_center_sub + SUB3(mi,mj,mk)][part] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + FETCH_DATA_COMPLEX8(input_array_ptr_complex8, STRIDE_IJK, part, + &data); break; } -#endif /* HAVE_CCTK_COMPLEX8 */ +#endif #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(mi) \ - input_array_ptr_complex16[input_center_sub + SUB1(mi)][part] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_complex16[input_center_sub + SUB2(mi,mj)][part] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_complex16[input_center_sub + SUB3(mi,mj,mk)][part] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + FETCH_DATA_COMPLEX16(input_array_ptr_complex16, + STRIDE_IJK, part, + &data); break; } -#endif /* HAVE_CCTK_COMPLEX16 */ +#endif + #ifdef HAVE_CCTK_COMPLEX32 case CCTK_VARIABLE_COMPLEX32: { + /* 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_complex32)[COMPLEX_N_PARTS] - = (const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr; - #undef DATA - #if (N_DIMS == 1) - #define DATA(mi) \ - input_array_ptr_complex32[input_center_sub + SUB1(mi)][part] - #elif (N_DIMS == 2) - #define DATA(mi,mj) \ - input_array_ptr_complex32[input_center_sub + SUB2(mi,mj)][part] - #elif (N_DIMS == 3) - #define DATA(mi,mj,mk) \ - input_array_ptr_complex32[input_center_sub + SUB3(mi,mj,mk)][part] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - #include DATA_VAR_ASSIGN_FILE_NAME + = ((const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + FETCH_DATA_COMPLEX32(input_array_ptr_complex32, + STRIDE_IJK, part, + &data); break; } -#endif /* HAVE_CCTK_COMPLEX32 */ +#endif default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): input datatype %d not supported!\n" " (0-origin) input #in=%d" - , - (int) input_array_type_codes[in], - in); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - /* end of switch */ - /* (input_array_type_codes[in]) */ - } - /* end of ***fetch*** input array values */ + , + (int) input_array_type_codes[in], + in); /*NOTREACHED*/ +return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + /* end of switch (input_array_type_codes[in]) */ } + } -if (output_arrays[out] != NULL) - then { - /* - * 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, + /* end of fetch input array values */ + } + + + if (output_arrays[out] != NULL) + then { + /* + * interpolate at this point + */ + fp result; + + switch (operation_codes[out]) + { + #ifdef HAVE_OP_I + case 0: + result = EVALUATE_MOLECULE(&coeffs_I, + &data); + break; + #endif + #ifdef HAVE_OP_DX + case 1: + result = inverse_dx + * EVALUATE_MOLECULE(&coeffs_dx, + &data); + break; + #endif + #ifdef HAVE_OP_DY + case 2: + result = inverse_dy + * EVALUATE_MOLECULE(&coeffs_dy, + &data); + break; + #endif + #ifdef HAVE_OP_DZ + case 3: + result = inverse_dz + * EVALUATE_MOLECULE(&coeffs_dz, + &data); + break; + #endif + #ifdef HAVE_OP_DXX + case 11: + result = inverse_dx * inverse_dx + * EVALUATE_MOLECULE(&coeffs_dxx, + &data); + break; + #endif + #ifdef HAVE_OP_DXY + case 12: + case 21: + result = inverse_dx * inverse_dy + * EVALUATE_MOLECULE(&coeffs_dxy, + &data); + break; + #endif + #ifdef HAVE_OP_DXZ + case 13: + case 31: + result = inverse_dx * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dxz, + &data); + break; + #endif + #ifdef HAVE_OP_DYY + case 22: + result = inverse_dy * inverse_dy + * EVALUATE_MOLECULE(&coeffs_dyy, + &data); + break; + #endif + #ifdef HAVE_OP_DYZ + case 23: + case 32: + result = inverse_dy * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dyz, + &data); + break; + #endif + #ifdef HAVE_OP_DZZ + case 33: + result = inverse_dz * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dzz, + &data); + break; + #endif + default: +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" " is not supported!\n" " (0-origin) output #out=%d" - , - (int) operation_codes[out], - out); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - /* end of switch (operation_codes[out]) */ - } + , + (int) operation_codes[out], + out); + 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]) - { + /* + * store result in 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; - } +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_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_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 +#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_REAL (*const output_array_ptr_complex)[COMPLEX_N_PARTS] - = (CCTK_REAL (*)[COMPLEX_N_PARTS]) output_arrays[out]; - output_array_ptr_complex[pt][part] = result; - break; - } +case CCTK_VARIABLE_COMPLEX: + { + CCTK_REAL (*const output_array_ptr_complex)[COMPLEX_N_PARTS] + = (CCTK_REAL (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex[pt][part] = result; + break; + } - #ifdef HAVE_CCTK_COMPLEX8 - case CCTK_VARIABLE_COMPLEX8: - { - CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS] - = (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out]; - output_array_ptr_complex8[pt][part] = result; - break; - } - #endif /* HAVE_CCTK_COMPLEX8 */ +#ifdef HAVE_CCTK_COMPLEX8 +case CCTK_VARIABLE_COMPLEX8: + { + CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS] + = (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex8[pt][part] = result; + break; + } +#endif /* HAVE_CCTK_COMPLEX8 */ - #ifdef HAVE_CCTK_COMPLEX16 - case CCTK_VARIABLE_COMPLEX16: - { - CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS] - = (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out]; - output_array_ptr_complex16[pt][part] = result; - break; - } - #endif /* HAVE_CCTK_COMPLEX16 */ +#ifdef HAVE_CCTK_COMPLEX16 +case CCTK_VARIABLE_COMPLEX16: + { + CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS] + = (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex16[pt][part] = result; + break; + } +#endif /* HAVE_CCTK_COMPLEX16 */ - #ifdef HAVE_CCTK_COMPLEX32 - case CCTK_VARIABLE_COMPLEX32: - { - CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS] - = (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out]; - output_array_ptr_complex32[pt][part] = result; - break; - } - #endif /* HAVE_CCTK_COMPLEX32 */ +#ifdef HAVE_CCTK_COMPLEX32 +case CCTK_VARIABLE_COMPLEX32: + { + CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS] + = (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex32[pt][part] = result; + break; + } +#endif /* HAVE_CCTK_COMPLEX32 */ - default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): output datatype %d not supported!\n" " (0-origin) output #out=%d" - , - (int) output_array_type_codes[out], - out); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - /* end of switch (output type code) */ - } - - /* end of if (output_arrays[out] != NULL) */ - } - + , + (int) output_array_type_codes[out], + out); + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + /* end of switch (output type code) */ + } -/* - * handle querying the Jacobian - */ -if ((Jacobian_info != NULL) && (Jacobian_info->Jacobian_pointer[out] != NULL)) - then { - /* - * set up the Jacobian storage addressing - */ - const int J_offset = Jacobian_info->Jacobian_offset[out]; - const int J_pt_stride = Jacobian_info->Jacobian_interp_point_stride; - const int J_part_stride = Jacobian_info->Jacobian_part_stride; - double *const Jacobian_ptr - = & Jacobian_info->Jacobian_pointer[out][J_offset - + pt*J_pt_stride - + part*J_part_stride]; - - #if (N_DIMS >= 1) - const int J_mx_stride = Jacobian_info->Jacobian_m_strides[X_AXIS]; - #endif - #if (N_DIMS >= 2) - const int J_my_stride = Jacobian_info->Jacobian_m_strides[Y_AXIS]; - #endif - #if (N_DIMS >= 3) - const int J_mz_stride = Jacobian_info->Jacobian_m_strides[Z_AXIS]; - #endif + /* end of if (output_arrays[out] != NULL) */ + } - #if (N_DIMS == 1) - #define COEFF(mi) Jacobian_ptr[mi*J_mx_stride] - #elif (N_DIMS == 2) - #define COEFF(mi,mj) Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride] - #elif (N_DIMS == 3) - #define COEFF(mi,mj,mk) \ - Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride + mk*J_mz_stride] - #else - #error "N_DIMS must be 1, 2, or 3!" - #endif - { - fp factor; - switch (operation_codes[out]) - { - #ifdef HAVE_OP_I - case 0: - #include COEFF_I_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DX - case 1: - factor = inverse_delta_x; - #include COEFF_DX_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DY - case 2: - factor = inverse_delta_y; - #include COEFF_DY_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DZ - case 3: - factor = inverse_delta_z; - #include COEFF_DZ_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DXX - case 11: - factor = inverse_delta_x2; - #include COEFF_DXX_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DXY - case 12: - case 21: - factor = inverse_delta_x * inverse_delta_y; - #include COEFF_DXY_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DXZ - case 13: - case 31: - factor = inverse_delta_x * inverse_delta_z; - #include COEFF_DXZ_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DYY - case 22: - factor = inverse_delta_y2; - #include COEFF_DYY_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DYZ - case 23: - case 32: - factor = inverse_delta_y * inverse_delta_z; - #include COEFF_DYZ_VAR_STORE_FILE_NAME - break; - #endif - #ifdef HAVE_OP_DZZ - case 33: - factor = inverse_delta_z2; - #include COEFF_DZZ_VAR_STORE_FILE_NAME - break; - #endif - default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + /* + * handle querying the Jacobian + */ + if ( (Jacobian_info != NULL) + && (Jacobian_info->Jacobian_pointer[out] != NULL)) + then { + CCTK_REAL *const Jacobian_ptr + = Jacobian_info->Jacobian_pointer[out] + + Jacobian_info->Jacobian_offset[out] + + Jacobian_info->Jacobian_interp_point_stride*pt + + Jacobian_info->Jacobian_part_stride*part; + + switch (operation_codes[out]) + { + #ifdef HAVE_OP_I + case 0: + STORE_COEFFS(1.0, &coeffs_I, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DX + case 1: + STORE_COEFFS(inverse_dx, &coeffs_dx, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DY + case 2: + STORE_COEFFS(inverse_dy, &coeffs_dy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DZ + case 3: + STORE_COEFFS(inverse_dz, &coeffs_dz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXX + case 11: + STORE_COEFFS(inverse_dx*inverse_dx, &coeffs_dxx, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXY + case 12: + case 21: + STORE_COEFFS(inverse_dx*inverse_dy, &coeffs_dxy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXZ + case 13: + case 31: + STORE_COEFFS(inverse_dx*inverse_dz, &coeffs_dxz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DYY + case 22: + STORE_COEFFS(inverse_dy*inverse_dy, &coeffs_dyy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DYZ + case 23: + case 32: + STORE_COEFFS(inverse_dy*inverse_dz, &coeffs_dyz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DZZ + case 33: + STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + default: +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" " is not supported!\n" " (0-origin) output #out=%d" - , - (int) operation_codes[out], - out); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - /* end of switch(operation_codes[out])*/ - } - } - /* end of Jacobian-query code */ - } - /* end of for (part = ...) loop */ + , + (int) operation_codes[out], + out); + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + /* end of switch(operation_codes[out])*/ + } + /* end of Jacobian-query code */ } - } - /* end of for (out = ...) loop */ + + /* end of for (part = ...) loop */ } + } + /* end of for (out = ...) loop */ + } } - } - } + } + } - /* end of for (pt = ...) loop */ - } + /* end of for (pt = ...) loop */ + } return 0; /*** NORMAL RETURN ***/ } } } + +/******************************************************************************/ + +/*@@ + @routine compute_coeffs_I + @routine compute_coeffs_dx + @routine compute_coeffs_dy + @routine compute_coeffs_dz + @routine compute_coeffs_dxx + @routine compute_coeffs_dxy + @routine compute_coeffs_dxz + @routine compute_coeffs_dyy + @routine compute_coeffs_dyz + @routine compute_coeffs_dzz + @date 29.Aug.2002 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc + Each of these functions computes the corresponding set of + interpolation coefficients at the point given by the XYZ + variables, using the machine-generated experessions in the + files + COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME + @enddesc + @@*/ + +#ifdef HAVE_OP_I + static + void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I) + { + #include COEFFS_I_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DX + static + void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx) + { + #include COEFFS_DX_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DY + static + void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy) + { + #include COEFFS_DY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DZ + static + void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz) + { + #include COEFFS_DZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXX + static + void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx) + { + #include COEFFS_DXX_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXY + static + void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy) + { + #include COEFFS_DXY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXZ + static + void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz) + { + #include COEFFS_DXZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DYY + static + void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy) + { + #include COEFFS_DYY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DYZ + static + void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz) + { + #include COEFFS_DYZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DZZ + static + void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz) + { + #include COEFFS_DZZ_COMPUTE_FILE_NAME + } +#endif |