Replace Cactus complex number type with C/C++ complex numbers

Map CCTK_COMPLEX to "double complex" in C, and "complex<double>" 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

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