Use new PUGHSlab API for passing the directions of a hyperslab.

Also support single-precision output of CCTK_COMPLEX datatypes.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IOHDF5Util/trunk@14 7842ec3a-9562-4be5-9c5b-06ba18f2b668

1 files changed, 88 insertions, 134 deletions

diff --git a/src/DumpVar.c b/src/DumpVar.c
index adfad6b..7ecf7e0 100644
--- a/src/DumpVar.c
+++ b/src/DumpVar.c
@@ -35,31 +35,32 @@ typedef struct
 } dumpInfo;
 
 
-/* local function prototypes */
+/* prototypes of routines defined in this source file */
 static int IOHDF5Util_DumpGS (cGH *GH,
                               int vindex,
                               int timelevel,
                               int iohdf5_type,
+                              ioHDF5Geo_t *request,  
                               int check_exisiting_objects,
                               hid_t file);
 static int IOHDF5Util_DumpGA (cGH *GH,
                               int vindex,
                               int timelevel,
-                              ioHDF5Geo_t *slab,  
+                              ioHDF5Geo_t *request,  
                               dumpInfo *info,
                               int check_exisiting_objects,
                               hid_t file);
 static int IOHDF5Util_getDumpData (cGH *GH,
                                    int vindex,
                                    int timelevel,
-                                   ioHDF5Geo_t *slab, 
+                                   ioHDF5Geo_t *request,  
                                    void **outme,
                                    int *free_outme, 
                                    CCTK_INT *geom);
 static void IOHDF5Util_procDump (cGH *GH,
                                  int vindex,
                                  int timelevel,
-                                 ioHDF5Geo_t *slab, 
+                                 ioHDF5Geo_t *request,  
                                  void *outme,
                                  CCTK_INT *geom,
                                  int proc, 
@@ -71,7 +72,7 @@ static void IOHDF5Util_procDump (cGH *GH,
 static void IOHDF5Util_collectiveDump (cGH *GH,
                                        int vindex,
                                        int timelevel,
-                                       ioHDF5Geo_t *slab,
+                                       ioHDF5Geo_t *request,  
                                        void *outme,
                                        CCTK_INT *geom, 
                                        int hdf5io_type,
@@ -107,8 +108,8 @@ static void IOHDF5Util_collectiveDump (cGH *GH,
    @vtype      int
    @vio        in
    @endvar
-   @var        slab
-   @vdesc      pointer to the hyperslab structure
+   @var        request
+   @vdesc      pointer to the IO request structure
    @vtype      ioHDF5Geo_t *
    @vio        in
    @endvar
@@ -136,7 +137,7 @@ static void IOHDF5Util_collectiveDump (cGH *GH,
 int IOHDF5Util_DumpVar (cGH *GH,
                         int vindex,
                         int timelevel,
-                        ioHDF5Geo_t *slab,
+                        ioHDF5Geo_t *request,
                         hid_t file,
                         int check_exisiting_objects)
 {
@@ -181,10 +182,13 @@ int IOHDF5Util_DumpVar (cGH *GH,
       break;
 
     case CCTK_VARIABLE_COMPLEX:
-      info.iohdf5_type  = myGH->IOHDF5_COMPLEX;
-      info.element_size = sizeof (CCTK_COMPLEX);
+      info.iohdf5_type  = ioUtilGH->out_single ? myGH->IOHDF5_COMPLEX8 :
+                                                 myGH->IOHDF5_COMPLEX;
+      info.element_size = ioUtilGH->out_single ?
+                              sizeof (CCTK_COMPLEX8) : sizeof (CCTK_COMPLEX);
 #ifdef CCTK_MPI
-      info.mpi_type     = pughGH->PUGH_mpi_complex;
+      info.mpi_type     = ioUtilGH->out_single ? pughGH->PUGH_mpi_complex8 :
+                                                 pughGH->PUGH_mpi_complex;
 #endif
       break;
 
@@ -200,11 +204,11 @@ int IOHDF5Util_DumpVar (cGH *GH,
   if (gtype == CCTK_SCALAR)
   {
     retval = IOHDF5Util_DumpGS (GH, vindex, timelevel, info.iohdf5_type,
-                                check_exisiting_objects, file);
+                                request, check_exisiting_objects, file);
   }
   else if (gtype == CCTK_ARRAY || gtype == CCTK_GF)
   {
-    retval = IOHDF5Util_DumpGA (GH, vindex, timelevel, slab, &info,
+    retval = IOHDF5Util_DumpGA (GH, vindex, timelevel, request, &info,
                                 check_exisiting_objects, file);
   }
   else
@@ -245,6 +249,11 @@ int IOHDF5Util_DumpVar (cGH *GH,
    @vtype      int
    @vio        in
    @endvar
+   @var        request
+   @vdesc      pointer to the IO request structure
+   @vtype      ioHDF5Geo_t *
+   @vio        in
+   @endvar
    @var        check_exisiting_objects
    @vdesc      flag indicating if we should check for already existing objects
                before these are created anew
@@ -267,6 +276,7 @@ static int IOHDF5Util_DumpGS (cGH *GH,
                               int vindex,
                               int timelevel,
                               int iohdf5_type,
+                              ioHDF5Geo_t *request,  
                               int check_exisiting_objects,
                               hid_t file)
 {
@@ -307,7 +317,7 @@ static int IOHDF5Util_DumpGS (cGH *GH,
                           CCTK_VarDataPtrI (GH, timelevel, vindex)));
     
   IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, global_shape,
-                                   NULL, dataset);
+                                   request, dataset);
 
   IOHDF5_ERROR (H5Dclose (dataset));
 
@@ -345,8 +355,8 @@ static int IOHDF5Util_DumpGS (cGH *GH,
    @vtype      int
    @vio        in
    @endvar
-   @var        slab
-   @vdesc      pointer to the hyperslab structure
+   @var        request
+   @vdesc      pointer to the IO request structure
    @vtype      ioHDF5Geo_t *
    @vio        in
    @endvar
@@ -380,7 +390,7 @@ static int IOHDF5Util_DumpGS (cGH *GH,
 static int IOHDF5Util_DumpGA (cGH *GH,
                               int vindex,
                               int timelevel,
-                              ioHDF5Geo_t *slab,  
+                              ioHDF5Geo_t *request,  
                               dumpInfo *info,
                               int check_exisiting_objects,
                               hid_t file)
@@ -408,10 +418,10 @@ static int IOHDF5Util_DumpGA (cGH *GH,
   ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
 
   /* allocate the geometry buffer */
-  geom = (CCTK_INT *) malloc (3 * slab->sdim * sizeof (CCTK_INT));
+  geom = (CCTK_INT *) malloc (3 * request->sdim * sizeof (CCTK_INT));
   
   /* Get the pointer to the data we want to output (includes downsampling) */
-  retval = IOHDF5Util_getDumpData (GH, vindex, timelevel, slab, &outme,
+  retval = IOHDF5Util_getDumpData (GH, vindex, timelevel, request, &outme,
                                    &free_outme, geom);
   if (retval < 0)
   {
@@ -425,7 +435,7 @@ static int IOHDF5Util_DumpGA (cGH *GH,
 #ifdef HAVE_PARALLEL
   if (ioUtilGH->unchunked)
   {
-    IOHDF5Util_collectiveDump (GH, vindex, timelevel, slab, outme, geom,
+    IOHDF5Util_collectiveDump (GH, vindex, timelevel, request, outme, geom,
                                info->iohdf5_type, file);
     if (free_outme)
     {
@@ -440,7 +450,7 @@ static int IOHDF5Util_DumpGA (cGH *GH,
   /* Dump data held on IO processor */
   if (myproc == ioUtilGH->ioproc)
   {
-    IOHDF5Util_procDump (GH, vindex, timelevel, slab, outme, geom,
+    IOHDF5Util_procDump (GH, vindex, timelevel, request, outme, geom,
                          myproc, info->iohdf5_type, check_exisiting_objects,
                          file);
   }
@@ -451,14 +461,14 @@ static int IOHDF5Util_DumpGA (cGH *GH,
     for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++)
     {
       /* receive geometry */
-      CACTUS_MPI_ERROR (MPI_Recv (geom, 3*slab->sdim, PUGH_MPI_INT, i,
+      CACTUS_MPI_ERROR (MPI_Recv (geom, 3*request->sdim, PUGH_MPI_INT, i,
                                   2*i + IOTAGBASE + 1, pughGH->PUGH_COMM_WORLD,
                                   &ms));
 
-      incoming = geom[slab->sdim];
-      for (j = 1; j < slab->sdim; j++)
+      incoming = geom[request->sdim];
+      for (j = 1; j < request->sdim; j++)
       {
-        incoming *= geom[slab->sdim + j];
+        incoming *= geom[request->sdim + j];
       }
 
       /* receive data */
@@ -469,7 +479,7 @@ static int IOHDF5Util_DumpGA (cGH *GH,
                                     2*i + IOTAGBASE, pughGH->PUGH_COMM_WORLD,
                                     &ms));
 
-        IOHDF5Util_procDump (GH, vindex, timelevel, slab, tmpd, geom, i,
+        IOHDF5Util_procDump (GH, vindex, timelevel, request, tmpd, geom, i,
                              info->iohdf5_type, check_exisiting_objects, file);
         free (tmpd);
       }
@@ -483,14 +493,14 @@ static int IOHDF5Util_DumpGA (cGH *GH,
     /* Send the geometry and data from the non-IO processors
      * to the IO processors
      */
-    outgoing = geom[slab->sdim];
-    for (i = 1; i < slab->sdim; i++)
+    outgoing = geom[request->sdim];
+    for (i = 1; i < request->sdim; i++)
     {
-      outgoing *= geom[slab->sdim + i];
+      outgoing *= geom[request->sdim + i];
     }
 
     /* send geometry */
-    CACTUS_MPI_ERROR (MPI_Send (geom, 3*slab->sdim, PUGH_MPI_INT,
+    CACTUS_MPI_ERROR (MPI_Send (geom, 3*request->sdim, PUGH_MPI_INT,
                                 ioUtilGH->ioproc, 2*myproc + IOTAGBASE + 1,
                                 pughGH->PUGH_COMM_WORLD));
     if (outgoing > 0)
@@ -526,71 +536,48 @@ static int IOHDF5Util_DumpGA (cGH *GH,
 static int IOHDF5Util_getDumpData (cGH *GH,
                                    int vindex,
                                    int timelevel,
-                                   ioHDF5Geo_t *slab, 
+                                   ioHDF5Geo_t *request, 
                                    void **outme,
                                    int *free_outme, 
                                    CCTK_INT *geom)
 {
+  int i;
+  int cctk_output_type;
   ioGH *ioUtilGH;
-  int myproc;
-  int i, sdim, vdim;             /*iteration, slab dim, variable dim */
   int *hsizes,*hsizes_global,*hsizes_offset; /* geometry information */
-  int sdir[3];                      /* slab direction FIXME */
-  int locdowns[3];                  /* Holds the downsampling in the order 
-                                       specified in slab->direction[] */
-  int slabstart[3];                 /* global start index for slab to extract */
-  void *data;
   char *fullname;
   
 
-  myproc = CCTK_MyProc (GH);
-
   /* get GH extension for IOUtil */
   ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
 
-  /* Shortcuts */
-  vdim = slab->vdim;
-  sdim = slab->sdim;
-  
-  /* FIXME */
-  if (vdim>3)
+  /* determine the datatype for the hyperslab */
+  cctk_output_type = CCTK_VarTypeI (vindex);
+  if (ioUtilGH->out_single)
   {
-    CCTK_WARN(1,"Cannot treat GFs with dim>3, check with cactus support\n");
-    return (-1);
+    if (cctk_output_type == CCTK_VARIABLE_INT)
+    {
+      cctk_output_type = CCTK_VARIABLE_INT4;
+    }
+    else if (cctk_output_type == CCTK_VARIABLE_REAL)
+    {
+      cctk_output_type = CCTK_VARIABLE_REAL4;
+    }
+    else if (cctk_output_type == CCTK_VARIABLE_COMPLEX)
+    {
+      cctk_output_type = CCTK_VARIABLE_COMPLEX8;
+    }
   }
 
-  for (i=0;i<sdim;i++)
-    locdowns[i]=slab->downsample[slab->direction[i]];
-
-  /* allocate array to hold size information of slab */
-  hsizes        = (int *) malloc (3 * sdim * sizeof (int));
-  hsizes_global = hsizes + 1*sdim;
-  hsizes_offset = hsizes + 2*sdim;
-
-  /* TEMPORARY FIX DUE TO INCONSISTENT DIR SPECIFICATION */
-  /* direction vector of 1d line in 3d volume */
-  if (sdim == 1)
-  {
-    sdir[0] = slab->direction[0] == 0;
-    sdir[1] = slab->direction[0] == 1;
-    sdir[2] = slab->direction[0] == 2;
-  }
-  /* norm vector of 2d surface in 3d volume */
-  else if (sdim==2)
-  {
-    sdir[0] = slab->direction[0] == 1 && slab->direction[1] == 2;
-    sdir[1] = slab->direction[0] == 0 && slab->direction[1] == 2;
-    sdir[2] = slab->direction[0] == 0 && slab->direction[1] == 1;
-  }
-  /* spanning directions for 3d */
-  else if (sdim==3)
-  {
-    sdir[0] = slab->direction[0];
-    sdir[1] = slab->direction[1];
-    sdir[2] = slab->direction[2];
-  }
-  
-  /* Get the start indeces for the slab from origin. */
+  /* FIXME */
+  /* allocate array to hold size information of request */
+  hsizes        = (int *) malloc (3 * request->sdim * sizeof (int));
+  hsizes_global = hsizes + 1*request->sdim;
+  hsizes_offset = hsizes + 2*request->sdim;
+
+#if 0
+  FIXME: what is this for ???
+  /* Get the start indeces for the request from origin. */
   /* Slab_start from the geo structure is not a true start index,
      it is currently used as the intersection point of three surfaces,
      that's why two entries have to be set to zero in the case of a
@@ -598,16 +585,19 @@ static int IOHDF5Util_getDumpData (cGH *GH,
      FIXME */
   for (i = 0; i < vdim; i++)
   {
-    slabstart[i] = slab->origin[i];
+    slabstart[i] = request->origin[i];
   }
-  for (i = 0; i < sdim; i++)
+  for (i = 0; i < request->sdim; i++)
   {
-    slabstart[slab->direction[i]] = 0;
+    slabstart[request->direction[i]] = 0;
   }     
+#endif
 
-  if (NewHyperslab_GetLocalHyperslab (GH, vindex, timelevel, sdim, CCTK_VARIABLE_REAL4, NULL, slab->origin,
-                                   sdir, slab->length, locdowns, outme,
-                                   hsizes, hsizes_global, hsizes_offset) < 0)
+  if (NewHyperslab_GetLocalHyperslab (GH, vindex, timelevel, request->sdim,
+                                      cctk_output_type, NULL, request->origin,
+                                      request->direction, request->length,
+                                      request->downsample, outme,
+                                      hsizes, hsizes_global, hsizes_offset) < 0)
   {
     fullname = CCTK_FullName (vindex);
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -619,52 +609,16 @@ static int IOHDF5Util_getDumpData (cGH *GH,
 
   *free_outme = 1;
 
-  for (i = 0; i < sdim; i++)
+  for (i = 0; i < request->sdim; i++)
   {
-    geom[i + 0*sdim] = hsizes_offset[i];
-    geom[i + 1*sdim] = hsizes[i];        
-    geom[i + 2*sdim] = hsizes_global[i];
-    slab->actlen[i]  = hsizes_global[i];
+    geom[i + 0*request->sdim] = hsizes_offset[i];
+    geom[i + 1*request->sdim] = hsizes[i];        
+    geom[i + 2*request->sdim] = hsizes_global[i];
+    request->actlen[i]  = hsizes_global[i];
   }
 
   free (hsizes);
 
-#ifdef DEBUG_ME
-  printf ("***** %dD \n",sdim);
-  printf ("Lower bound: %d", (int) geom[0]);
-  for (i = 1; i < sdim; i++)
-    printf (" %d", (int) geom[i]);
-  printf ("\n");
-  printf ("Chunk size : %d", (int) geom[1*sdim + 0]);
-  for (i = 1; i < sdim; i++)
-    printf (" %d", (int) geom[1*sdim + i]);
-  printf ("\n");
-  printf ("Global size: %d", (int) geom[2*sdim + 0]);
-  for (i = 1; i < sdim; i++) 
-    printf (" %d", (int) geom[2*sdim + i]);
-  printf ("\n");
-  printf ("Downsampling: ");
-  for (i = 0; i < sdim; i++)
-    printf (" %d", locdowns[i]);
-  printf ("\n"); 
-  printf ("Direction: "); 
-  for (i = 0; i < sdim; i++)
-    printf (" %d",slab->direction[i]);
-  printf ("\n");  
-  printf ("Slab Start: %d", (int) slabstart[0]);
-  for (i = 1; i < vdim; i++)
-    printf (" %d", (int) slabstart[i]);
-  printf ("\n"); 
-  printf ("Slab intersect cntr: %d", (int) slab->origin[0]);
-  for (i = 1; i < vdim; i++)
-    printf (" %d", (int) slab->origin[i]);
-  printf ("\n");  
-  printf ("SlabDelta: "); 
-  for (i = 0; i < sdim; i++)
-    printf (" %f",GH->cctk_delta_space[slab->direction[i]]*locdowns[i]);
-  printf ("\n");
-#endif
-
   return (0);
 }
 
@@ -728,7 +682,7 @@ static int IOHDF5Util_getDumpData (cGH *GH,
 static void IOHDF5Util_procDump (cGH *GH,
                                  int vindex,
                                  int timelevel,
-                                 ioHDF5Geo_t *slab,
+                                 ioHDF5Geo_t *request,
                                  void *outme,
                                  CCTK_INT *geom,
                                  int proc, 
@@ -752,7 +706,7 @@ static void IOHDF5Util_procDump (cGH *GH,
   ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
   myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util");
 
-  sdim = slab->sdim;
+  sdim = request->sdim;
   myproc = CCTK_MyProc (GH);
 
   /* Check if there are points to output */
@@ -811,7 +765,7 @@ static void IOHDF5Util_procDump (cGH *GH,
       IOHDF5_ERROR (dataset = H5Dcreate (file, datasetname,
                                          iohdf5_type, filespace, H5P_DEFAULT));
       IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
-                                       slab, dataset);
+                                       request, dataset);
     }
     else
     {
@@ -841,7 +795,7 @@ static void IOHDF5Util_procDump (cGH *GH,
     {
       IOHDF5_ERROR (group = H5Gcreate (file, datasetname, 0));
       IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
-                                       slab, group);
+                                       request, group);
       IOHDF5_ERROR (H5Gclose (group));
     }
 
@@ -927,7 +881,7 @@ static void IOHDF5Util_procDump (cGH *GH,
 static void IOHDF5Util_collectiveDump (cGH *GH,
                                        int vindex,
                                        int timelevel,
-                                       ioHDF5Geo_t *slab,  
+                                       ioHDF5Geo_t *request,  
                                        void *outme,
                                        CCTK_INT *geom, 
                                        int hdf5io_type,
@@ -946,7 +900,7 @@ static void IOHDF5Util_collectiveDump (cGH *GH,
   myGH = (ioHDF5GH *) CCTK_GHExtension (GH, "IOHDF5Util");
 
   /* get the dimension of the variable */
-  sdim = slab->sdim;
+  sdim = request->sdim;
 
   chunk_origin = (hssize_t *) malloc (sdim * sizeof (hssize_t));
   chunk_dims   = (hsize_t *)  malloc (2*sdim * sizeof (hsize_t));
@@ -989,7 +943,7 @@ static void IOHDF5Util_collectiveDump (cGH *GH,
   if (CCTK_MyProc (GH) == 0)
   {
     IOHDF5Util_DumpCommonAttributes (GH, vindex, timelevel, &geom[2*sdim],
-                                     slab, dataset);
+                                     request, dataset);
   }
 
   /* increase the buffer size if the default isn't sufficient */