diff options
Diffstat (limited to 'src/Interpolate.c')
-rw-r--r-- | src/Interpolate.c | 52 |
1 files changed, 19 insertions, 33 deletions
diff --git a/src/Interpolate.c b/src/Interpolate.c index 748d48a..56df7d0 100644 --- a/src/Interpolate.c +++ b/src/Interpolate.c @@ -133,14 +133,12 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c) * different directions' interpolations so that only 3 scalar temporaries * are needed. */ -#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array, \ +#define INTERPOLATE(cctk_type, in_array, out_array, \ order, point, dims, n, coeff) \ { \ int ii, jj, kk; \ const cctk_type *fi; \ - /* for some reason (probably loop unrolling & pipelining) the compiler \ - produces faster code if arrays are used instead of scalars */ \ - cctk_subtype interp_result, fj[1], fk[1]; \ + cctk_type interp_result, fj, fk; \ \ \ interp_result = 0; \ @@ -148,43 +146,39 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c) /* NOTE-MAXDIM: support >3D arrays by adding more loops */ \ for (kk = 0; kk <= order; kk++) \ { \ - fk[0] = 0; \ + fk = 0; \ for (jj = 0; jj <= order; jj++) \ { \ /* NOTE-MAXDIM: for >3D arrays adapt the index calculation here */ \ fi = (const cctk_type *) in_array + \ point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk)); \ \ - fj[0] = 0; \ + fj = 0; \ for (ii = 0; ii <= order; ii++) \ { \ - fj[0] += fi[ii]subpart * coeff[0][ii]; \ + fj += fi[ii] * coeff[0][ii]; \ } \ /* at this point we have just computed */ \ - /* fj[0] = in_array[*][jj][kk] interpolated to x=offset[0] */ \ + /* fj = in_array[*][jj][kk] interpolated to x=offset[0] */ \ \ - fk[0] += fj[0] * coeff[1][jj]; \ + fk += fj * coeff[1][jj]; \ } \ /* at this point we have just computed */ \ - /* fk[0] = fj[0][*][kk] interpolated to y=offset[1] */ \ + /* fk = fj[*][kk] interpolated to y=offset[1] */ \ /* = in_array[*][*][kk] interpolated to */ \ /* x=offset[0], y=offset[1] */ \ \ - interp_result += fk[0] * coeff[2][kk]; \ + interp_result += fk * coeff[2][kk]; \ } \ /* at this point we have just computed */ \ - /* interp_result = fk[0][*] interpolated to z=offset[2] */ \ + /* interp_result = fk[*] interpolated to z=offset[2] */ \ /* = in_array[*][*][*] interpolated to */ \ /* x=offset[0], y=offset[1], z=offset[2] */ \ \ /* assign the result */ \ - ((cctk_type *) out_array)[n]subpart = interp_result; \ + ((cctk_type *) out_array)[n] = interp_result; \ } /* end of macro */ -/* this is needed by some preprocessors to pass into INTERPOLATE - as a dummy macro */ -#define NOTHING - /******************************************************************************/ /*@@ @@ -606,55 +600,47 @@ if (n == LocalInterp_verbose_debug_n) we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */ if (in_types[a] == CCTK_VARIABLE_REAL) { - INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a], + INTERPOLATE (CCTK_REAL, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } else if (in_types[a] == CCTK_VARIABLE_COMPLEX) { - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a], + INTERPOLATE (CCTK_COMPLEX, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } #ifdef CCTK_REAL4 else if (in_types[a] == CCTK_VARIABLE_REAL4) { - INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a], + INTERPOLATE (CCTK_REAL4, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } else if (in_types[a] == CCTK_VARIABLE_COMPLEX8) { - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a], + INTERPOLATE (CCTK_COMPLEX8, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } #endif #ifdef CCTK_REAL8 else if (in_types[a] == CCTK_VARIABLE_REAL8) { - INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a], + INTERPOLATE (CCTK_REAL8, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } else if (in_types[a] == CCTK_VARIABLE_COMPLEX16) { - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a], + INTERPOLATE (CCTK_COMPLEX16, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } #endif #ifdef CCTK_REAL16 else if (in_types[a] == CCTK_VARIABLE_REAL16) { - INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a], + INTERPOLATE (CCTK_REAL16, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } else if (in_types[a] == CCTK_VARIABLE_COMPLEX32) { - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a], + INTERPOLATE (CCTK_COMPLEX32, in_arrays[a], out_arrays[a], order, point, max_dims, n, coeff); } #endif |