Switch to the new hyperslabbing API.

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

3 files changed, 99 insertions, 146 deletions

diff --git a/src/DumpUtils.c b/src/DumpUtils.c
index 4ab5688..2045d1a 100644
--- a/src/DumpUtils.c
+++ b/src/DumpUtils.c
@@ -56,6 +56,7 @@ CCTK_FILEVERSION(BetaThorns_IOHDF5Util_DumpUtils_c)
 int IOHDF5Util_DumpGH (const cGH *GH, const int *timers, hid_t file)
 {
   int first_vindex, gindex, timelevels, retval;
+  int *extent_int;
   cGroup gdata;
   ioSlab slab;
   DECLARE_CCTK_PARAMETERS
@@ -65,21 +66,17 @@ int IOHDF5Util_DumpGH (const cGH *GH, const int *timers, hid_t file)
 
   /* set up a hyperslab description for full hyperslabs */
   slab.vdim = CCTK_MaxDim ();
-  slab.vectors       = (CCTK_INT *) calloc (slab.vdim + 6,
-                                            slab.vdim * sizeof(CCTK_INT));
+  slab.vectors      = (CCTK_INT *) calloc ((slab.vdim + 6)*slab.vdim + 1,
+                                           sizeof(CCTK_INT));
   slab.hoffset      = slab.vectors + 0*slab.vdim;
   slab.hsize        = slab.vectors + 1*slab.vdim;
   slab.hsize_chunk  = slab.vectors + 2*slab.vdim;
-  slab.origin       = slab.vectors + 3*slab.vdim;
-  slab.extent       = slab.vectors + 4*slab.vdim;
-  slab.downsample   = slab.vectors + 5*slab.vdim;
-  slab.direction    = slab.vectors + 6*slab.vdim;
+  slab.origin       = slab.vectors + 3*slab.vdim + 1;
+  slab.extent       = slab.vectors + 4*slab.vdim + 1;
+  slab.downsample   = slab.vectors + 5*slab.vdim + 1;
+  slab.direction    = slab.vectors + 6*slab.vdim + 1;
 
-  for (slab.hdim = 0; slab.hdim < slab.vdim; slab.hdim++)
-  {
-    slab.extent[slab.hdim]     = -1;
-    slab.downsample[slab.hdim] = 1;
-  }
+  extent_int = (int *) malloc (slab.vdim * sizeof (int));
 
   /* start CP_PARAMETERS_TIMER timer */
   if (timers)
@@ -121,20 +118,24 @@ int IOHDF5Util_DumpGH (const cGH *GH, const int *timers, hid_t file)
 
     /* dump all timelevels except the oldest (for multi-level groups) */
     CCTK_GroupData (gindex, &gdata);
+    CCTK_GroupgshGI (GH, gdata.dim, extent_int, gindex);
     timelevels = gdata.numtimelevels;
     if (timelevels > 1)
     {
       timelevels--;
     }
-    slab.hdim = slab.vdim = gdata.dim;
+    slab.vdim = gdata.dim;
     slab.check_exist = 0;
     slab.hdatatype = gdata.vartype;
 
-    /* set the hyperslab directions (orthogonal to all axes) */
-    memset (slab.direction, 0, slab.hdim * slab.vdim * sizeof (int));
+    /* set the hyperslab extents, the directions (orthogonal to all axes),
+       and disable downsampling */
+    memset (slab.direction, 0, slab.vdim * slab.vdim * sizeof (int));
     for (slab.hdim = 0; slab.hdim < slab.vdim; slab.hdim++)
     {
+      slab.extent[slab.hdim] = extent_int[slab.hdim];
       slab.direction[slab.hdim * (slab.vdim + 1)] = 1;
+      slab.downsample[slab.hdim] = 1;
     }
 
     /* loop over all variables in this group */
@@ -165,6 +166,7 @@ int IOHDF5Util_DumpGH (const cGH *GH, const int *timers, hid_t file)
 
   /* free temporary resources */
   free (slab.vectors);
+  free (extent_int);
 
   return (retval);
 }
diff --git a/src/DumpVar.c b/src/DumpVar.c
index 96148d9..6492f58 100644
--- a/src/DumpVar.c
+++ b/src/DumpVar.c
@@ -14,8 +14,8 @@
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
+#include "Hyperslab.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
-#include "CactusPUGH/PUGHSlab/src/NewPUGHSlab.h"
 #include "CactusBase/IOUtil/src/ioGH.h"
 #include "ioHDF5UtilGH.h"
 
@@ -41,8 +41,6 @@ static int WriteGS (const cGH *GH, const ioSlab *slab, const char *name,
                     hid_t file);
 static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
                     hid_t file);
-static int GetLocalHyperslab (const cGH *GH, const ioSlab *slab, void **data,
-                              int *free_data);
 static void WriteData (const cGH *GH, const ioSlab *slab, const char *name,
                        const void *data, int proc, hid_t file);
 #if defined(CCTK_MPI) && defined(H5_HAVE_PARALLEL)
@@ -196,6 +194,8 @@ static int WriteGS (const cGH *GH, const ioSlab *slab, const char *name,
   HDF5_ERROR (H5Dwrite (dataset, hdf5type, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                         CCTK_VarDataPtrI (GH, slab->timelevel, slab->vindex)));
 
+  /* scalars have size 0 */
+  slab->hsize[0] = 0;
   IOHDF5Util_DumpCommonAttributes (GH, slab, dataset);
 
   HDF5_ERROR (H5Dclose (dataset));
@@ -212,7 +212,9 @@ static int WriteGS (const cGH *GH, const ioSlab *slab, const char *name,
                Writes a grid array into a HDF5 file.
    @enddesc
 
-   @calls      GetLocalHyperslab
+   @calls      Hyperslab_DefineLocalMappingByIndex
+               Hyperslab_FreeMapping
+               Hyperslab_Get
                WriteData
                WriteDataCollective
 
@@ -239,8 +241,9 @@ static int WriteGS (const cGH *GH, const ioSlab *slab, const char *name,
 
    @returntype int
    @returndesc
-               0 for success
-               or returncode of @seeroutine GetLocalHyperslab
+                0 for success, or<BR>
+               -1 if hyperslab mapping couldn't be defined, or<BR>
+               -2 if hyperslab couldn't be extracted
    @endreturndesc
 @@*/
 static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
@@ -248,12 +251,11 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
 {
   ioGH *ioUtilGH;
   ioHDF5UtilGH *myGH;
-  int myproc, nprocs, retval;
-  void *data;    /* The data pointer to dump ... */
-  int free_data; /* and whether it needs freeing */
+  int i, myproc, nprocs, mapping, hdatasize, retval;
+  void *hdata;
+  char *fullname;
 #ifdef CCTK_MPI
   void *tmpd;
-  int i, j, incoming, outgoing;
   MPI_Comm comm;
   MPI_Status ms;
   MPI_Datatype mpitype;
@@ -264,28 +266,63 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
   myGH = (ioHDF5UtilGH *) CCTK_GHExtension (GH, "IOHDF5Util");
   ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
 
-#ifdef CCTK_MPI
-  mpitype = PUGH_MPIDataType (PUGH_pGH (GH), slab->hdatatype);
-#endif
+  myproc = CCTK_MyProc (GH);
+  nprocs = CCTK_nProcs (GH);
 
-  /* get the pointer to the data we want to output (includes downsampling) */
-  retval = GetLocalHyperslab (GH, slab, &data, &free_data);
-  if (retval < 0)
+  /* define the hyperslab mapping */
+  mapping = Hyperslab_DefineLocalMappingByIndex (GH, slab->vindex, slab->hdim,
+                                                 slab->direction, slab->origin,
+                                                 slab->extent, slab->downsample,
+                                                 -1, NULL, slab->hsize_chunk,
+                                                 slab->hsize, slab->hoffset);
+  if (mapping < 0)
   {
-    return (retval);
+    fullname = CCTK_FullName (slab->vindex);
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Failed to define hyperslab mapping for variable '%s'",
+                fullname);
+    free (fullname);
+    return (-1);
   }
 
-  myproc = CCTK_MyProc (GH);
-  nprocs = CCTK_nProcs (GH);
+  /* calculate the size of the hyperslab */
+  slab->hsize_chunk[slab->hdim] = 1;
+  for (i = 0; i < slab->hdim; i++)
+  {
+    slab->hsize_chunk[slab->hdim] *= slab->hsize_chunk[i];
+  }
+
+  /* get the hyperslab */
+  hdatasize = CCTK_VarTypeSize (slab->hdatatype);
+  hdata = slab->hsize_chunk[slab->hdim] > 0 ?
+          malloc (slab->hsize_chunk[slab->hdim] * hdatasize) : NULL;
+  retval = Hyperslab_Get (GH, mapping, -1, slab->vindex, slab->timelevel,
+                          slab->hdatatype, hdata);
+
+  /* release the mapping structure */
+  Hyperslab_FreeMapping (mapping);
+
+  if (retval)
+  {
+    fullname = CCTK_FullName (slab->vindex);
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Failed to extract hyperslab for variable '%s'", fullname);
+    free (fullname);
+    if (hdata)
+    {
+      free (hdata);
+    }
+    return (-2);
+  }
 
 #ifdef CCTK_MPI
 #ifdef H5_HAVE_PARALLEL
   if (ioUtilGH->unchunked)
   {
-    WriteDataCollective (GH, slab, name, data, file);
-    if (free_data)
+    WriteDataCollective (GH, slab, name, hdata, file);
+    if (hdata)
     {
-      free (data);
+      free (hdata);
     }
     return (0);
   }
@@ -295,10 +332,12 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
   /* dump data held on I/O processor */
   if (myproc == ioUtilGH->ioproc)
   {
-    WriteData (GH, slab, name, data, myproc, file);
+    WriteData (GH, slab, name, hdata, myproc, file);
   }
 #ifdef CCTK_MPI
   comm = PUGH_pGH (GH)->PUGH_COMM_WORLD;
+  mpitype = PUGH_MPIDataType (PUGH_pGH (GH), slab->hdatatype);
+
   if (myproc == ioUtilGH->ioproc)
   {
     /* dump data from all other processors */
@@ -306,22 +345,16 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
     {
       /* receive geometry (this assumes the geometry arrays
          to be contiguous starting at hoffset */
-      CACTUS_MPI_ERROR (MPI_Recv (slab->hoffset, 3*slab->hdim, PUGH_MPI_INT, i,
-                                  2*i + MPITAGBASE + 1, comm, &ms));
-
-      incoming = 1;
-      for (j = 0; j < slab->hdim; j++)
-      {
-        incoming *= slab->hsize_chunk[slab->hdim + j];
-      }
+      CACTUS_MPI_ERROR (MPI_Recv (slab->hoffset, 3*slab->hdim + 1, PUGH_MPI_INT,
+                                  i, 2*i + MPITAGBASE + 1, comm, &ms));
 
       /* receive data */
       tmpd = NULL;
-      if (incoming > 0)
+      if (slab->hsize_chunk[slab->hdim] > 0)
       {
-        tmpd = malloc (incoming * CCTK_VarTypeSize (slab->hdatatype));
-        CACTUS_MPI_ERROR (MPI_Recv (tmpd, incoming, mpitype, i,
-                                    2*i + MPITAGBASE, comm, &ms));
+        tmpd = malloc (slab->hsize_chunk[slab->hdim] * hdatasize);
+        CACTUS_MPI_ERROR (MPI_Recv (tmpd, slab->hsize_chunk[slab->hdim],
+                                    mpitype, i, 2*i + MPITAGBASE, comm, &ms));
       }
 
       /* write data */
@@ -337,28 +370,20 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
   }
   else
   {
-    /*
-     * send the geometry and data from non-I/O processors to the I/O processors
-     */
-    outgoing = 1;
-    for (i = 0; i < slab->hdim; i++)
-    {
-      outgoing *= slab->hsize_chunk[slab->hdim + i];
-    }
-
-    /* send geometry (this assumes the geometry arrays
-       to be contiguous starting at hoffset */
-    CACTUS_MPI_ERROR (MPI_Send (slab->hoffset, 3*slab->hdim, PUGH_MPI_INT,
+    /* send geometry (this assumes the geometry arrays to be contiguous
+       starting at hoffset) */
+    CACTUS_MPI_ERROR (MPI_Send (slab->hoffset, 3*slab->hdim + 1, PUGH_MPI_INT,
                                 ioUtilGH->ioproc, 2*myproc + MPITAGBASE + 1,
                                 comm));
     /* send data */
-    if (outgoing > 0)
+    if (slab->hsize_chunk[slab->hdim] > 0)
     {
-      CACTUS_MPI_ERROR (MPI_Send (data, outgoing, mpitype,
+      CACTUS_MPI_ERROR (MPI_Send (hdata, slab->hsize_chunk[slab->hdim], mpitype,
                                   ioUtilGH->ioproc, 2*myproc + MPITAGBASE,
                                   comm));
 #ifdef IOHDF5UTIL_DEBUG
-      printf ("Processor %d sent %d data points\n", myproc, outgoing);
+      printf ("Processor %d sent %d data points\n",
+              myproc, slab->hsize_chunk[slab->hdim]);
 #endif
     }
   }
@@ -368,89 +393,15 @@ static int WriteGA (const cGH *GH, const ioSlab *slab, const char *name,
 #endif
 
   /* free allocated resources */
-  if (free_data)
+  if (hdata)
   {
-    free (data);
+    free (hdata);
   }
 
   return (retval);
 }
 
 
-static int GetLocalHyperslab (const cGH *GH, const ioSlab *slab, void **data,
-                              int *free_data)
-{
-  int i, retval;
-  char *fullname;
-
-
-#if 0
-  mapping = Hyperslab_DefineLocalMappingByIndex (GH, vindex, slab->hdim,
-                                                 slab->direction, slab->origin,
-                                                 slab->extent, slab->downsample,
-                                                 -1, NULL, slab->hsize_local,
-                                                 slab->hsize_global,
-                                                 slab->hoffset_global);
-           const cGH          *GH,
-           CCTK_INT            vindex,
-           CCTK_INT            hdim,
-           const CCTK_INT     *direction         /* vdim*hdim */,
-           const CCTK_INT     *origin            /* vdim */,
-           const CCTK_INT     *extent            /* hdim */,
-           const CCTK_INT     *downsample        /* hdim */,
-           CCTK_INT            table_handle,
-           t_hslabConversionFn conversion_fn,
-           CCTK_INT           *hsize_local,      /* hdim */
-           CCTK_INT           *hsize_global,     /* hdim */
-           CCTK_INT           *hoffset_global    /* hdim */);
-
-  if (mapping < 0)
-  {
-    fullname = CCTK_FullName (vindex);
-    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "Failed to extract hyperslab for variable '%s'", fullname);
-    free (fullname);
-  }
-
-#else
-{
-  int *hsizes, *hsizes_global, *hsizes_offset; /* geometry information */
-
-  /* allocate array to hold size information of request */
-  hsizes        = (int *) malloc (3 * slab->hdim * sizeof (int));
-  hsizes_global = hsizes + 1*slab->hdim;
-  hsizes_offset = hsizes + 2*slab->hdim;
-
-  retval = NewHyperslab_GetLocalHyperslab (GH, slab->vindex, slab->timelevel, slab->hdim,
-                                      slab->hdatatype, NULL, slab->origin,
-                                      slab->direction, slab->extent,
-                                      slab->downsample, data, free_data,
-                                      hsizes, hsizes_global, hsizes_offset);
-  if (retval == 0)
-  {
-    for (i = 0; i < slab->hdim; i++)
-    {
-      slab->hoffset[i] = hsizes_offset[i];
-      slab->hsize[i] = hsizes_global[i];
-      slab->hsize_chunk[i] = hsizes[i];
-    }
-  }
-  else
-  {
-    fullname = CCTK_FullName (slab->vindex);
-    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "Failed to extract hyperslab for variable '%s'", fullname);
-    free (fullname);
-  }
-
-  free (hsizes);
-}
-#endif
-
-  return (retval);
-}
-
-
 /*@@
    @routine    WriteData
    @author     Thomas Radke
diff --git a/src/ParseVars.c b/src/ParseVars.c
index 4718bbe..2c58703 100644
--- a/src/ParseVars.c
+++ b/src/ParseVars.c
@@ -370,15 +370,15 @@ static ioSlab *DefaultIOHyperslab (const cGH *GH, int vindex)
 
   /* allocate the arrays all in one go
      only initialize those which are mandatory for the Hyperslab API */
-  slab->vectors = (CCTK_INT *) calloc (slab->vdim + 6,
-                                       slab->vdim * sizeof (CCTK_INT));
+  slab->vectors = (CCTK_INT *) calloc ((slab->vdim + 6) * slab->vdim + 1,
+                                       sizeof (CCTK_INT));
   slab->hoffset      = slab->vectors + 0*slab->vdim;
   slab->hsize        = slab->vectors + 1*slab->vdim;
   slab->hsize_chunk  = slab->vectors + 2*slab->vdim;
-  slab->origin       = slab->vectors + 3*slab->vdim;
-  slab->extent       = slab->vectors + 4*slab->vdim;
-  slab->downsample   = slab->vectors + 5*slab->vdim;
-  slab->direction    = slab->vectors + 6*slab->vdim;
+  slab->origin       = slab->vectors + 3*slab->vdim + 1;
+  slab->extent       = slab->vectors + 4*slab->vdim + 1;
+  slab->downsample   = slab->vectors + 5*slab->vdim + 1;
+  slab->direction    = slab->vectors + 6*slab->vdim + 1;
 
   for (slab->hdim = 0; slab->hdim < slab->vdim; slab->hdim++)
   {