From 474719cb822661afa915df69051886edd65f3de5 Mon Sep 17 00:00:00 2001 From: eschnett Date: Fri, 8 Mar 2013 20:31:02 +0000 Subject: Replace Cactus complex number type with C/C++ complex numbers Map CCTK_COMPLEX to "double complex" in C, and "complex" in C++. (It is already mapped to "double complex" in Fortran.) Update type definitions. Re-implement Cactus complex number math functions by calling the respective C functions. Update thorn that access real and imaginary parts of complex numbers to use standard-conforming methods instead. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/LocalInterp/trunk@216 df1f8a13-aa1d-4dd4-9681-27ded5b42416 --- src/Interpolate.c | 52 +++++++++++++++++++--------------------------------- 1 file changed, 19 insertions(+), 33 deletions(-) (limited to 'src') 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 -- cgit v1.2.3