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.

4 files changed, 35 insertions, 74 deletions

diff --git a/Carpet/CarpetIOBasic/src/iobasic.cc b/Carpet/CarpetIOBasic/src/iobasic.cc
index 28abd6454..f8366a811 100644
--- a/Carpet/CarpetIOBasic/src/iobasic.cc
+++ b/Carpet/CarpetIOBasic/src/iobasic.cc
@@ -482,14 +482,10 @@ namespace CarpetIOBasic {
         
           int const handle = ireduction->handle;
         
-          union {
-#define TYPECASE(N,T) T var_##T;
-#include "typecase.hh"
-#undef TYPECASE
-          } result;
+          char result[100];     // assuming no type is larger than this
         
           int const ierr
-            = CCTK_Reduce (cctkGH, 0, handle, 1, vartype, &result, 1, n);
+            = CCTK_Reduce (cctkGH, 0, handle, 1, vartype, result, 1, n);
           assert (not ierr);
         
           if (CCTK_MyProc(cctkGH) == 0) {
@@ -500,7 +496,7 @@ namespace CarpetIOBasic {
 #define TYPECASE(N,T)                                                   \
               case N:                                                   \
                 {                                                       \
-                  T const val = result.var_##T;                         \
+                  T const& val = *(T const*)result;                     \
                   if (not isint) {                                      \
                     if (UseScientificNotation (val)) {                  \
                       cout << scientific << setprecision(real_prec_sci); \
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
index 927c9740d..feccb8fbe 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
@@ -445,32 +445,32 @@ static void* SetupGH (tFleshConfig* const fleshconfig,
   HDF5_ERROR (myGH->HDF5_COMPLEX =
               H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX)));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "real",
-                         offsetof (CCTK_COMPLEX, Re), HDF5_REAL));
+                         0, HDF5_REAL));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX, "imag",
-                         offsetof (CCTK_COMPLEX, Im), HDF5_REAL));
+                         sizeof (CCTK_REAL), HDF5_REAL));
 #ifdef CCTK_REAL4
   HDF5_ERROR (myGH->HDF5_COMPLEX8 =
               H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX8)));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "real",
-                         offsetof (CCTK_COMPLEX8, Re), H5T_NATIVE_FLOAT));
+                         0, H5T_NATIVE_FLOAT));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX8, "imag",
-                         offsetof (CCTK_COMPLEX8, Im), H5T_NATIVE_FLOAT));
+                         sizeof (CCTK_REAL4), H5T_NATIVE_FLOAT));
 #endif
 #ifdef CCTK_REAL8
   HDF5_ERROR (myGH->HDF5_COMPLEX16 =
               H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX16)));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "real",
-                         offsetof (CCTK_COMPLEX16, Re), H5T_NATIVE_DOUBLE));
+                         0, H5T_NATIVE_DOUBLE));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX16, "imag",
-                         offsetof (CCTK_COMPLEX16, Im), H5T_NATIVE_DOUBLE));
+                         sizeof (CCTK_REAL8), H5T_NATIVE_DOUBLE));
 #endif
 #ifdef CCTK_REAL16
   HDF5_ERROR (myGH->HDF5_COMPLEX32 =
               H5Tcreate (H5T_COMPOUND, sizeof (CCTK_COMPLEX32)));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "real",
-                         offsetof (CCTK_COMPLEX32, Re), H5T_NATIVE_LDOUBLE));
+                         0, H5T_NATIVE_LDOUBLE));
   HDF5_ERROR (H5Tinsert (myGH->HDF5_COMPLEX32, "imag",
-                         offsetof (CCTK_COMPLEX32, Im), H5T_NATIVE_LDOUBLE));
+                         sizeof (CCTK_REAL16), H5T_NATIVE_LDOUBLE));
 #endif
 
   return (myGH);
diff --git a/Carpet/CarpetIOScalar/src/ioscalar.cc b/Carpet/CarpetIOScalar/src/ioscalar.cc
index 6dbb9eaae..1cd0ee102 100644
--- a/Carpet/CarpetIOScalar/src/ioscalar.cc
+++ b/Carpet/CarpetIOScalar/src/ioscalar.cc
@@ -447,11 +447,7 @@ namespace CarpetIOScalar {
         
         int const handle = ireduction->handle;
 
-        union {
-#define TYPECASE(N,T) T var_##T;
-#include "typecase.hh"
-#undef TYPECASE
-        } result;
+        char result[100];       // assuming no type is larger
 
         int const firstvar
           = one_file_per_group ? CCTK_FirstVarIndexI(group) : n;
@@ -461,14 +457,14 @@ namespace CarpetIOScalar {
         for (int n=firstvar; n<firstvar+numvars; ++n) {
           
           int const ierr
-            = CCTK_Reduce (cctkGH, 0, handle, 1, vartype, &result, 1, n);
+            = CCTK_Reduce (cctkGH, 0, handle, 1, vartype, result, 1, n);
           if (ierr) {
             char * const fullname = CCTK_FullName (n);
             CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                         "Error during reduction for variable \"%s\"",
                         fullname);
             free (fullname);
-            memset (&result, 0, sizeof result);
+            memset (result, 0, sizeof result);
           }
           
           if (CCTK_MyProc(cctkGH)==0) {
@@ -477,7 +473,7 @@ namespace CarpetIOScalar {
             switch (specific_cactus_type(vartype)) {
 #define TYPECASE(N,T)                           \
               case N:                           \
-                file << result.var_##T;         \
+                file << *(T const*)result;      \
               break;
 #include "typecase.hh"
 #undef TYPECASE
diff --git a/Carpet/CarpetLib/src/dist.cc b/Carpet/CarpetLib/src/dist.cc
index 49b73a078..8546ca6ca 100644
--- a/Carpet/CarpetLib/src/dist.cc
+++ b/Carpet/CarpetLib/src/dist.cc
@@ -63,16 +63,16 @@ namespace dist {
 #define CARPET_NO_COMPLEX
 #include "typecase.hh"
 #undef TYPECASE
-#define TYPECASE(N,T)                                           \
-    if (datatype == mpi_datatype<T>()) {                        \
-      assert(not done);                                         \
-      done = true;                                              \
-      T *restrict const invec = (T*)invec_;                     \
-      T *restrict const inoutvec = (T*)inoutvec_;               \
-      for (int n=0; n<len; ++n) {                               \
-        inoutvec[n].Re = max(inoutvec[n].Re, invec[n].Re);      \
-        inoutvec[n].Im = max(inoutvec[n].Im, invec[n].Im);      \
-      }                                                         \
+#define TYPECASE(N,T)                                                   \
+    if (datatype == mpi_datatype<T>()) {                                \
+      assert(not done);                                                 \
+      done = true;                                                      \
+      T *restrict const invec = (T*)invec_;                             \
+      T *restrict const inoutvec = (T*)inoutvec_;                       \
+      for (int n=0; n<len; ++n) {                                       \
+        inoutvec[n] = T(max(real(inoutvec[n]), real(invec[n])),         \
+                        max(imag(inoutvec[n]), imag(invec[n])));        \
+      }                                                                 \
     }
 #define CARPET_COMPLEX
 #include "typecase.hh"
@@ -101,16 +101,16 @@ namespace dist {
 #define CARPET_NO_COMPLEX
 #include "typecase.hh"
 #undef TYPECASE
-#define TYPECASE(N,T)                                           \
-    if (datatype == mpi_datatype<T>()) {                        \
-      assert(not done);                                         \
-      done = true;                                              \
-      T *restrict const invec = (T*)invec_;                     \
-      T *restrict const inoutvec = (T*)inoutvec_;               \
-      for (int n=0; n<len; ++n) {                               \
-        inoutvec[n].Re = min(inoutvec[n].Re, invec[n].Re);      \
-        inoutvec[n].Im = min(inoutvec[n].Im, invec[n].Im);      \
-      }                                                         \
+#define TYPECASE(N,T)                                                   \
+    if (datatype == mpi_datatype<T>()) {                                \
+      assert(not done);                                                 \
+      done = true;                                                      \
+      T *restrict const invec = (T*)invec_;                             \
+      T *restrict const inoutvec = (T*)inoutvec_;                       \
+      for (int n=0; n<len; ++n) {                                       \
+        inoutvec[n] = T(min(real(inoutvec[n]), real(invec[n])),         \
+                        min(imag(inoutvec[n]), imag(invec[n])));        \
+      }                                                                 \
     }
 #define CARPET_COMPLEX
 #include "typecase.hh"
@@ -136,22 +136,6 @@ namespace dist {
         inoutvec[n] *= invec[n];                        \
       }                                                 \
     }
-#define CARPET_NO_COMPLEX
-#include "typecase.hh"
-#undef TYPECASE
-#define TYPECASE(N,T)                                   \
-    if (datatype == mpi_datatype<T>()) {                \
-      assert(not done);                                 \
-      done = true;                                      \
-      T *restrict const invec = (T*)invec_;             \
-      T *restrict const inoutvec = (T*)inoutvec_;       \
-      for (int n=0; n<len; ++n) {                       \
-        complex<T>& inout = *(complex<T>*)&inoutvec[n]; \
-        complex<T>& in    = *(complex<T>*)&invec[n];    \
-        inout *= in;                                    \
-      }                                                 \
-    }
-#define CARPET_COMPLEX
 #include "typecase.hh"
 #undef TYPECASE
     assert(done);
@@ -175,21 +159,6 @@ namespace dist {
         inoutvec[n] += invec[n];                        \
       }                                                 \
     }
-#define CARPET_NO_COMPLEX
-#include "typecase.hh"
-#undef TYPECASE
-#define TYPECASE(N,T)                                   \
-    if (datatype == mpi_datatype<T>()) {                \
-      assert(not done);                                 \
-      done = true;                                      \
-      T *restrict const invec = (T*)invec_;             \
-      T *restrict const inoutvec = (T*)inoutvec_;       \
-      for (int n=0; n<len; ++n) {                       \
-        inoutvec[n].Re += invec[n].Re;                  \
-        inoutvec[n].Im += invec[n].Im;                  \
-      }                                                 \
-    }
-#define CARPET_COMPLEX
 #include "typecase.hh"
 #undef TYPECASE
     assert(done);