Removed the restriction that the slab types have to be CCTK_REAL.

Boy, is juggling Cactus, MPI, and C types a mess.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Slab/trunk@12 2e825fa2-fb71-486d-8b7f-a5ff3f0f6cb8

3 files changed, 265 insertions, 41 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index c1709e0..27a0f0f 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -328,11 +328,11 @@ are not valid, but will correctly fill in the ghost zones in the
 destination slab.
 
 There are currently some restrictions, which can be remedied if the
-need arises: The datatype of the arrays must be CCTK\_REAL.  The
-dimension must be 3 (lower dimensions can easily be padded).  The
-communication schedule is set up for every slab transfer, which is
-expensive.  The number of ghost zones along the lower and upper
-boundary is assumed to be the same.  There is no Fortran interface.
+need arises: The dimension must be 3 (lower dimensions can easily be
+padded).  The communication schedule is set up for every slab
+transfer, which is expensive.  The number of ghost zones along the
+lower and upper boundary is assumed to be the same.  There is no
+Fortran interface.
 \end{Discussion}
 
 \begin{SeeAlsoSection}
diff --git a/schedule.ccl b/schedule.ccl
index 2a19d26..07ea006 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,2 +1,7 @@
 # Schedule definitions for thorn Slab
 # $Header$
+
+SCHEDULE Slab_InitMPIDatatypes AT startup after Driver_Startup
+{
+  LANG: C
+} "Create MPI datatypes for complex variables in C"
diff --git a/src/slab.c b/src/slab.c
index dc6714b..03b7190 100644
--- a/src/slab.c
+++ b/src/slab.c
@@ -2,7 +2,6 @@
 
 /*
   TODO:
-  Provide facilities for types other than CCTK_REAL
   Provide facilities for dim != 3
   Set up the slab exchange information in advance
   Test slabbing without MPI
@@ -72,14 +71,90 @@ CCTK_FILEVERSION(TAT_Slab_slab_c);
 
 
 
+#ifdef CCTK_MPI
+/* Determine MPI type sizes */
+
+#  define MPI_INT1 MPI_CHAR
+
+#  if SIZEOF_SHORT_INT == 2
+#    define MPI_INT2 MPI_SHORT
+#  elif SIZEOF_INT == 2
+#    define MPI_INT2 MPI_INT
+#  elif SIZEOF_LONG_INT == 2
+#    define MPI_INT2 MPI_LONG
+#  elif SIZEOF_LONG_LONG == 2
+#    define MPI_INT2 MPI_LONG_LONG_INT
+#  endif
+
+#  if SIZEOF_SHORT_INT == 4
+#    define MPI_INT4 MPI_SHORT
+#  elif SIZEOF_INT == 4
+#    define MPI_INT4 MPI_INT
+#  elif SIZEOF_LONG_INT == 4
+#    define MPI_INT4 MPI_LONG
+#  elif SIZEOF_LONG_LONG == 4
+#    define MPI_INT4 MPI_LONG_LONG_INT
+#  endif
+
+#  if SIZEOF_SHORT_INT == 8
+#    define MPI_INT8 MPI_SHORT
+#  elif SIZEOF_INT == 8
+#    define MPI_INT8 MPI_INT
+#  elif SIZEOF_LONG_INT == 8
+#    define MPI_INT8 MPI_LONG
+#  elif SIZEOF_LONG_LONG == 8
+#    define MPI_INT8 MPI_LONG_LONG_INT
+#  endif
+
+#  if SIZEOF_FLOAT == 4
+#    define MPI_REAL4 MPI_FLOAT
+#  elif SIZEOF_DOUBLE == 4
+#    define MPI_REAL4 MPI_DOUBLE
+#  elif SIZEOF_LONG_DOUBLE == 4
+#    define MPI_REAL4 MPI_LONG_DOUBLE
+#  endif
+
+#  if SIZEOF_FLOAT == 8
+#    define MPI_REAL8 MPI_FLOAT
+#  elif SIZEOF_DOUBLE == 8
+#    define MPI_REAL8 MPI_DOUBLE
+#  elif SIZEOF_LONG_DOUBLE == 8
+#    define MPI_REAL8 MPI_LONG_DOUBLE
+#  endif
+
+#  if SIZEOF_FLOAT == 16
+#    define MPI_REAL16 MPI_FLOAT
+#  elif SIZEOF_DOUBLE == 16
+#    define MPI_REAL16 MPI_DOUBLE
+#  elif SIZEOF_LONG_DOUBLE == 16
+#    define MPI_REAL16 MPI_LONG_DOUBLE
+#  endif
+
+static MPI_Datatype MPI_COMPLEX8  = -1;
+static MPI_Datatype MPI_COMPLEX16 = -1;
+static MPI_Datatype MPI_COMPLEX32 = -1;
+
+#endif
+
+
+
 /* Replace MPI functions if MPI is disabled */
 #ifndef CCTK_MPI
 
+typedef int MPI_Comm;
 typedef int MPI_Datatype;
-#define MPI_DOUBLE CCTK_VARIABLE_REAL8
-#define MPI_INT CCTK_VARIABLE_INT4
 
-typedef int MPI_Comm;
+#define MPI_BYTE      CCTK_VARIABLE_BYTE
+#define MPI_INT1      CCTK_VARIABLE_INT1
+#define MPI_INT2      CCTK_VARIABLE_INT2
+#define MPI_INT4      CCTK_VARIABLE_INT4
+#define MPI_INT8      CCTK_VARIABLE_INT8
+#define MPI_REAL4     CCTK_VARIABLE_REAL4
+#define MPI_REAL8     CCTK_VARIABLE_REAL8
+#define MPI_REAL16    CCTK_VARIABLE_REAL16
+#define MPI_COMPLEX8  CCTK_VARIABLE_COMPLEX8
+#define MPI_COMPLEX16 CCTK_VARIABLE_COMPLEX16
+#define MPI_COMPLEX32 CCTK_VARIABLE_COMPLEX32
 
 static int
 MPI_Barrier (MPI_Comm comm)
@@ -159,7 +234,9 @@ MPI_Alltoallv (void * sendbuf, int * sendcnt, int * sendoff, int sendtype,
 
 
 
-static MPI_Comm get_mpi_comm (cGH * restrict const cctkGH)
+/* Get the MPI COMM_WOLRD communicator from the driver */
+static MPI_Comm
+get_mpi_comm (cGH * restrict const cctkGH)
 {
 #ifdef CCTK_MPI
 #  if defined CARPET_CARPET
@@ -173,6 +250,7 @@ static MPI_Comm get_mpi_comm (cGH * restrict const cctkGH)
   }
 #  endif
   assert (0);
+  return -1;
 #else
   return 0;
 #endif
@@ -180,6 +258,125 @@ static MPI_Comm get_mpi_comm (cGH * restrict const cctkGH)
 
 
 
+/* Initialise the MPI datatypes for complex variables */
+void
+Slab_InitMPIDatatypes (void)
+{
+#ifdef CCTK_MPI
+#  ifdef CCTK_REAL4
+  MPI_Type_contiguous (2, MPI_REAL4, &MPI_COMPLEX8);
+  MPI_Type_commit (&MPI_COMPLEX8);
+#  endif
+#  ifdef CCTK_REAL8
+  MPI_Type_contiguous (2, MPI_REAL8, &MPI_COMPLEX16);
+  MPI_Type_commit (&MPI_COMPLEX16);
+#  endif
+#  ifdef CCTK_REAL16
+  MPI_Type_contiguous (2, MPI_REAL16, &MPI_COMPLEX32);
+  MPI_Type_commit (&MPI_COMPLEX32);
+#  endif
+#endif
+}
+
+
+
+/* Normalise a Cactus datatype */
+static int
+normal_type (int cactustype)
+{
+  switch (cactustype) {
+  case CCTK_VARIABLE_INT:
+#ifdef CCTK_INTEGER_PRECISION_1
+    return CCTK_VARIABLE_INT1;
+#endif
+#ifdef CCTK_INTEGER_PRECISION_2
+    return CCTK_VARIABLE_INT2;
+#endif
+#ifdef CCTK_INTEGER_PRECISION_4
+    return CCTK_VARIABLE_INT4;
+#endif
+#ifdef CCTK_INTEGER_PRECISION_8
+    return CCTK_VARIABLE_INT8;
+#endif
+    assert (0);
+    return -1;
+  case CCTK_VARIABLE_REAL:
+#ifdef CCTK_REAL_PRECISION_4
+    return CCTK_VARIABLE_REAL4;
+#endif
+#ifdef CCTK_REAL_PRECISION_8
+    return CCTK_VARIABLE_REAL8;
+#endif
+#ifdef CCTK_REAL_PRECISION_16
+    return CCTK_VARIABLE_REAL16;
+#endif
+    assert (0);
+    return -1;
+  case CCTK_VARIABLE_COMPLEX:
+#ifdef CCTK_REAL_PRECISION_4
+    return CCTK_VARIABLE_COMPLEX8;
+#endif
+#ifdef CCTK_REAL_PRECISION_8
+    return CCTK_VARIABLE_COMPLEX16;
+#endif
+#ifdef CCTK_REAL_PRECISION_16
+    return CCTK_VARIABLE_COMPLEX32;
+#endif
+    assert (0);
+    return -1;
+  }
+  return cactustype;
+}
+
+
+
+/* Find the MPI datatype corresponding to a Cactus datatype */
+static int mpi_type (int const cactustype)
+{
+  int const normaltype = normal_type (cactustype);
+  switch (normaltype) {
+  case CCTK_VARIABLE_BYTE: return MPI_BYTE;
+#if defined CCTK_INT1 && defined MPI_INT1
+  case CCTK_VARIABLE_INT1: return MPI_INT1;
+#endif
+#if defined CCTK_INT2 && defined MPI_INT2
+  case CCTK_VARIABLE_INT2: return MPI_INT2;
+#endif
+#if defined CCTK_INT4 && defined MPI_INT4
+  case CCTK_VARIABLE_INT4: return MPI_INT4;
+#endif
+#if defined CCTK_INT8 && defined MPI_INT8
+  case CCTK_VARIABLE_INT8: return MPI_INT8;
+#endif
+#if defined CCTK_REAL4 && defined MPI_REAL4
+  case CCTK_VARIABLE_REAL4: return MPI_REAL4;
+#endif
+#if defined CCTK_REAL8 && defined MPI_REAL8
+  case CCTK_VARIABLE_REAL8: return MPI_REAL8;
+#endif
+#if defined CCTK_REAL16 && defined MPI_REAL16
+  case CCTK_VARIABLE_REAL16: return MPI_REAL16;
+#endif
+#if defined CCTK_REAL4 && defined CCTK_COMPLEX8
+  case CCTK_VARIABLE_COMPLEX8: return MPI_COMPLEX8;
+#endif
+#if defined CCTK_REAL4 && defined CCTK_COMPLEX16
+  case CCTK_VARIABLE_COMPLEX16: return MPI_COMPLEX16;
+#endif
+#if defined CCTK_REAL16 && defined CCTK_COMPLEX32
+  case CCTK_VARIABLE_COMPLEX32: return MPI_COMPLEX32;
+#endif
+  }
+  assert (0);
+  return -1;
+}
+
+
+
+/*****************************************************************************/
+
+
+
 struct bbox {
   int off, len, str;
 };
@@ -373,6 +570,8 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
   int * restrict dstcount;
   int * restrict dstoffset;
   
+  int srctypesize;
+  int dsttypesize;
   void * restrict srcdata;
   void * restrict dstdata;
   
@@ -388,9 +587,9 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
   assert (cctkGH);
   assert (dim >= 0);
   assert (xferinfo);
-  assert (srctype == CCTK_VARIABLE_REAL);
+  assert (srctype >= 0);
   assert (srcptr);
-  assert (dsttype == CCTK_VARIABLE_REAL);
+  assert (dsttype >= 0);
   assert (dstptr);
   
   info = malloc (dim * sizeof *info);
@@ -452,9 +651,15 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
   MPI_Comm_size (comm, &size);
   MPI_Comm_rank (comm, &rank);
   
-  assert (sizeof(CCTK_REAL) == sizeof(double));
-  srcdatatype = MPI_DOUBLE;
-  dstdatatype = MPI_DOUBLE;
+  assert (srctype == dsttype);
+  srctypesize = CCTK_VarTypeSize (srctype);
+  assert (srctypesize > 0);
+  dsttypesize = CCTK_VarTypeSize (dsttype);
+  assert (dsttypesize > 0);
+  srcdatatype = mpi_type (srctype);
+  assert (srcdatatype >= 0);
+  dstdatatype = mpi_type (dsttype);
+  assert (dstdatatype >= 0);
   
   
   
@@ -542,14 +747,16 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
       ("srccnt n=%d offset=%d count=%d\n", n, srcoffset[n], srccount[n]);
     srcoffset[n+1] = srcoffset[n] + srccount[n];
   }
-  srcdata = malloc (srcoffset[size] * sizeof(CCTK_REAL));
+  srcdata = malloc (srcoffset[size] * srctypesize);
   assert (srcoffset[size] == 0 || srcdata);
   ifcheck {
-    CCTK_REAL * restrict const srcdataptr = srcdata;
-    CCTK_REAL marker;
-    memset (&marker, -1, sizeof marker);
-    for (i = 0; i < srcoffset[size]; ++i) {
-      memcpy (&srcdataptr[i], &marker, sizeof marker);
+    if (srctype == CCTK_VARIABLE_REAL) {
+      CCTK_REAL * restrict const srcdataptr = srcdata;
+      CCTK_REAL marker;
+      memset (&marker, -1, sizeof marker);
+      for (i = 0; i < srcoffset[size]; ++i) {
+	memcpy (&srcdataptr[i], &marker, sizeof marker);
+      }
     }
   }
   
@@ -604,14 +811,16 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
       ("dstcnt n=%d offset=%d count=%d\n", n, dstoffset[n], dstcount[n]);
     dstoffset[n+1] = dstoffset[n] + dstcount[n];
   }
-  dstdata = malloc (dstoffset[size] * sizeof(CCTK_REAL));
+  dstdata = malloc (dstoffset[size] * dsttypesize);
   assert (dstoffset[size] == 0 || dstdata);
   ifcheck {
-    CCTK_REAL * restrict const dstdataptr = dstdata;
-    CCTK_REAL marker;
-    memset (&marker, -1, sizeof marker);
-    for (i = 0; i < dstoffset[size]; ++i) {
-      memcpy (&dstdataptr[i], &marker, sizeof marker);
+    if (dsttype == CCTK_VARIABLE_REAL) {
+      CCTK_REAL * restrict const dstdataptr = dstdata;
+      CCTK_REAL marker;
+      memset (&marker, -1, sizeof marker);
+      for (i = 0; i < dstoffset[size]; ++i) {
+	memcpy (&dstdataptr[i], &marker, sizeof marker);
+      }
     }
   }
   
@@ -671,19 +880,24 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
 	  }
 	  assert (srcind < srclentot);
 	  assert (bufind < (size_t)srccount[n]);
-	  ((CCTK_REAL*)srcdata)[srcoffset[n] + bufind]
-	    = ((const CCTK_REAL*)srcptr)[srcind];
+/* 	  ((CCTK_REAL*)srcdata)[srcoffset[n] + bufind] */
+/* 	    = ((const CCTK_REAL*)srcptr)[srcind]; */
+	  memcpy ((char *) srcdata + srctypesize * (srcoffset[n] + bufind),
+		  (const char *) srcptr + srctypesize * srcind,
+		  srctypesize);
 	}
       }
     }
   }
   
   ifcheck {
-    const CCTK_REAL * restrict const srcdataptr = srcdata;
-    CCTK_REAL marker;
-    memset (&marker, -1, sizeof marker);
-    for (i = 0; i < srcoffset[size]; ++i) {
-      assert (memcmp(&srcdataptr[i], &marker, sizeof marker) != 0);
+    if (srctype == CCTK_VARIABLE_REAL) {
+      const CCTK_REAL * restrict const srcdataptr = srcdata;
+      CCTK_REAL marker;
+      memset (&marker, -1, sizeof marker);
+      for (i = 0; i < srcoffset[size]; ++i) {
+	assert (memcmp(&srcdataptr[i], &marker, sizeof marker) != 0);
+      }
     }
   }
   
@@ -693,11 +907,13 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
      dstdata, dstcount, dstoffset, dstdatatype, comm);
   
   ifcheck {
-    const CCTK_REAL * restrict const dstdataptr = dstdata;
-    CCTK_REAL marker;
-    memset (&marker, -1, sizeof marker);
-    for (i = 0; i < dstoffset[size]; ++i) {
-      assert (memcmp(&dstdataptr[i], &marker, sizeof marker) != 0);
+    if (dsttype == CCTK_VARIABLE_REAL) {
+      const CCTK_REAL * restrict const dstdataptr = dstdata;
+      CCTK_REAL marker;
+      memset (&marker, -1, sizeof marker);
+      for (i = 0; i < dstoffset[size]; ++i) {
+	assert (memcmp(&dstdataptr[i], &marker, sizeof marker) != 0);
+      }
     }
   }
   
@@ -736,8 +952,11 @@ int Slab_Transfer (cGH                   * restrict const cctkGH,
 	  }
 	  assert (bufind < (size_t)dstcount[n]);
 	  assert (dstind < dstlentot);
-	  ((CCTK_REAL*)dstptr)[dstind]
-	    = ((const CCTK_REAL*)dstdata)[dstoffset[n] + bufind];
+/* 	  ((CCTK_REAL*)dstptr)[dstind] */
+/* 	    = ((const CCTK_REAL*)dstdata)[dstoffset[n] + bufind]; */
+	  memcpy ((char *) dstptr + dsttypesize * dstind,
+		  (const char *) dstdata + dsttypesize * (dstoffset[n] + bufind),
+		  dsttypesize);
 	}
       }
     }