Forgot to cvs add this important source file.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHSlab/trunk@23 10716dce-81a3-4424-a2c8-48026a0d3035

1 files changed, 930 insertions, 0 deletions

diff --git a/src/Hyperslab.c b/src/Hyperslab.c
new file mode 100644
index 0000000..8b28b66
--- /dev/null
+++ b/src/Hyperslab.c
@@ -0,0 +1,930 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "Hyperslab.h"
+
+
+/* enable debugging output */
+/*#define DEBUG 1*/
+
+/* macros for getting the minimum/maximum value */
+#ifdef MIN
+#undef MIN
+#endif
+#define MIN(x, y)             ((x) < (y) ? (x) : (y))
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(x, y)             ((x) > (y) ? (x) : (y))
+
+/* shortcuts for the local/global start/endpoints on this processor */
+#define MY_LOCAL_SP(extras, istag, dim)    (extras->ownership [istag][0][dim])
+#define MY_LOCAL_EP(extras, istag, dim)    (extras->ownership [istag][1][dim])
+#define MY_GLOBAL_SP(extras, myproc, istag, dim)                              \
+        (extras->lb [myproc][dim] + MY_LOCAL_SP (extras, istag, dim))
+#define MY_GLOBAL_EP(extras, myproc, istag, dim)                              \
+        (extras->lb [myproc][dim] + MY_LOCAL_EP (extras, istag, dim))
+
+
+/* macro for copying all relevant data points
+   (ranging from startpoint[] to endpoint[])
+   into the typed hyperslab data buffer */
+#define PICKUP_HYPERSLAB_DATA(cctk_type, vdim, vdata, hdata, point,           \
+                              startpoint, endpoint, downsample,points_per_dim)\
+    {                                                                         \
+      int i, dim, vindex;                                                     \
+      cctk_type *typed_vdata = (cctk_type *) vdata;                           \
+      cctk_type *typed_hdata = (cctk_type *) hdata;                           \
+                                                                              \
+                                                                              \
+      /* set point to local startpoint */                                     \
+      memcpy (point, startpoint, vdim * sizeof (point [0]));                  \
+                                                                              \
+      /* do the nested loops starting with the innermost */                   \
+      dim = 0;                                                                \
+      while (1)                                                               \
+      {                                                                       \
+        /* check for end of current loop */                                   \
+        if (point [dim] >= endpoint [dim])                                    \
+        {                                                                     \
+          /* increment outermost loopers */                                   \
+          for (dim++; dim < vdim; dim++)                                      \
+          {                                                                   \
+            point [dim] += downsample [dim];                                  \
+            if (point [dim] < endpoint [dim])                                 \
+              break;                                                          \
+          }                                                                   \
+                                                                              \
+          /* done if beyond outermost loop */                                 \
+          if (dim >= vdim)                                                    \
+            break;                                                            \
+                                                                              \
+          /* reset innermost loopers */                                       \
+          for (dim--; dim >= 0; dim--)                                        \
+            point [dim] = startpoint [dim];                                   \
+          dim = 0;                                                            \
+        }                                                                     \
+                                                                              \
+        /* copy the data point */                                             \
+        vindex = point [0];                                                   \
+        for (i = 1; i < vdim; i++)                                            \
+          vindex += point [i] * points_per_dim [i];                           \
+        *typed_hdata++ = typed_vdata [vindex];                                \
+                                                                              \
+        /* increment current looper */                                        \
+        point [dim] += downsample [dim];                                      \
+                                                                              \
+      } /* end of nested loops over all dimensions */                         \
+    }
+
+
+/*** local function prototypes ***/
+/* Gerd's routine to get 1D lines from a 3D array */
+int Hyperslab_CollectData1D (cGH *GH, int vindex, int vtimelvl,
+                             const int *origin,
+                             const int *directions,
+                             int downsampling,
+                             int length,
+                             void **hdata,
+                             int *hsize,
+                             int proc);
+/* routine to check parameters passed to the Hyperslab routines */
+static const char *checkParameters (cGH *GH, int vindex, int vtimelvl,
+                                    int hdim,
+                                    const int global_startpoint [/*vdim*/],
+                                    const int directions [/*vdim*/],
+                                    const int extents [/*hdim*/],
+                                    const int downsample_ [/*hdim*/],
+                                    void **hdata,
+                                    int hsize [/*hdim*/]);
+
+
+/*@@
+  @routine    Hyperslab_GetLocalHyperslab
+  @date       Fri May 12 2000
+  @author     Thomas Radke
+  @desc
+              Extract a hyperslab from the processor-local chunk
+              of a domain-decomposed Cactus array variable.
+
+              This routine delivers the local hyperslab data
+              to be collected into a global hyperslab
+              by Hyperslab_GetHyperslab().
+              IO methods can call this routine as well to collect the
+              local hyperslab data and output it in parallel.
+  @enddesc 
+  @calls
+  @calledby   Hyperslab_GetHyperslab
+  @history 
+  @endhistory 
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        vindex
+  @vdesc      index of variable to get a hyperslab from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        vtimelvl
+  @vdesc      timelvl of variable to get a hyperslab from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        hdim
+  @vdesc      dimensionality of the requested hyperslab
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        global_startpoint
+  @vdesc      global coordinates of the hyperslab origin
+  @vtype      const int [dimensions of vindex]
+  @vio        in
+  @endvar
+  @var        directions
+  @vdesc      directions which span the hyperslab
+              This is the direction vector for 1D hyperslabs,
+              or the normal vector of the 2 lowest dimensions
+              for hyperslabs of higher dimensions.
+              If the hyperslab is of same dimensionality as the variable
+              the directions vector is ignored.
+  @vtype      const int [dimensions of vindex]
+  @vio        in
+  @endvar
+  @var        extents
+  @vdesc      number of grid points to follow in each hyperslab direction
+              starting from origin
+              Negative values are taken as extents up to the grid boundaries.
+  @vtype      const int [hdim]
+  @vio        in
+  @endvar
+  @var        downsample_
+  @vdesc      downsampling values for each hyperslab dimension
+  @vtype      const int [hdim]
+  @vio        in
+  @endvar
+  @var        hdata
+  @vdesc      pointer to store the address of the hyperslab data buffer
+  @vtype      void **
+  @vio        out
+  @endvar
+  @var        hsize
+  @vdesc      sizes of the (local) hyperslab data buffer in each dimension
+  @vtype      int [hdim]
+  @vio        out
+  @endvar
+  @var        hsize_global
+  @vdesc      sizes of the global hyperslab data buffer in each dimension
+  @vtype      int [hdim]
+  @vio        out
+  @endvar
+  @var        hoffset_global
+  @vdesc      if not NULL, array to save the offsets of the local hyperslab
+              into the global one for each dimension
+  @vtype      int [hdim]
+  @vio        out
+  @endvar
+@@*/
+int Hyperslab_GetLocalHyperslab (cGH *GH, int vindex, int vtimelvl,
+                                 int hdim,
+                                 const int global_startpoint [/*vdim*/],
+                                 const int directions [/*vdim*/],
+                                 const int extents [/*hdim*/],
+                                 const int downsample_ [/*hdim*/],
+                                 void **hdata,
+                                 int hsize [/*hdim*/],
+                                 int hsize_global [/*hdim*/],
+                                 int hoffset_global [/*hdim*/])
+{
+  int *point;                    /* looper over hyperslab dimensions */
+  int *startpoint,               /* hyperslab's local start and endpoint */
+      *endpoint;                 /* within the variable's grid dimensions */
+  int *downsample;               /* the downsample_[] vector extended to vdim */
+  int *my_global_startpoint,     /* hyperslab's global start and endpoint */
+      *global_endpoint;
+  int *points_per_dim;           /* points per subvolume */
+  int *do_dir;                   /* directions in which to span the hyperslab */
+  int stagger_index;             /* stagger index in direction i */
+  int myproc;                    /* local processor ID */
+  cGroup vinfo;                  /* variable's group info */
+  int vdim;                      /* looper over all dimensions */
+  int totals;                    /* total number of hyperslab data points */
+  int retval;                    /* the return value (0 for success) */
+  pGExtras *extras;              /* the variable's info structure from PUGH */
+  const char *errormsg;          /* explanation string in case of errors */
+
+
+  /* do some plausibility checks */
+  errormsg = checkParameters (GH, vindex, vtimelvl, hdim,
+                              global_startpoint, directions, extents,
+                              downsample_, hdata, hsize);
+
+  /* immediately return in case of errors */
+  if (errormsg)
+  {
+    CCTK_WARN (1, errormsg);
+    return (-1);
+  }
+
+  /* set the default return code */
+  retval = 0;
+
+  /* get the info on the variable to extract a hyperslab from */
+  CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
+
+  /* allocate the temporary arrays */
+  point = (int *) malloc (8 * vinfo.dim * sizeof (int));
+  startpoint           = point + 1*vinfo.dim;
+  endpoint             = point + 2*vinfo.dim;
+  my_global_startpoint = point + 3*vinfo.dim;
+  global_endpoint      = point + 4*vinfo.dim;
+  downsample           = point + 5*vinfo.dim;
+  do_dir               = point + 6*vinfo.dim;
+  points_per_dim       = point + 7*vinfo.dim;
+
+  /* get the actual directions in which to spawn the hyperslab */
+  /* Note: hdim == 1    the normal vector is the directions vector itself
+           hdim <  vdim the normal vector is the negated directions vector
+           hdim == vdim the normal vector is completely ignored,
+                        and all directions are selected */
+  if (hdim == 1)
+    memcpy (do_dir, directions, vinfo.dim * sizeof (int));
+  else if (hdim < vinfo.dim)
+  {
+    for (vdim = 0; vdim < vinfo.dim; vdim++)
+      do_dir [vdim] = ! directions [vdim];
+  }
+  else
+  {
+    for (vdim = 0; vdim < vinfo.dim; vdim++)
+      do_dir [vdim] = 1;
+  }
+
+  /* get the variable's info structure from PUGH */
+  extras = ((pGA *) PUGH_pGH (GH)->variables [vindex][vtimelvl])->extras;
+
+  /* compute the global endpoint */
+  for (vdim = hdim = 0; vdim < vinfo.dim; vdim++)
+  {
+    if (do_dir [vdim])
+    {
+      global_endpoint [vdim] = extents [hdim] > 0 ?
+                               MIN (global_startpoint [vdim] + extents [hdim],
+                                    extras->nsize [vdim]) :
+                               extras->nsize [vdim];
+      downsample [vdim] = downsample_ [hdim];
+      hdim++;
+    }
+    else
+    {
+      global_endpoint [vdim] = global_startpoint [vdim] + 1;
+      downsample [vdim] = 1;
+    }
+  }
+
+  /* get the local processor ID */
+  myproc = CCTK_MyProc (GH);
+
+  /* compute my global startpoint from the global ranges */
+  for (vdim = 0; vdim < vinfo.dim; vdim++)
+  {
+    stagger_index = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+
+    if (global_startpoint [vdim] < MY_GLOBAL_EP (extras, myproc,
+                                                 stagger_index, vdim))
+    {
+      if (global_startpoint [vdim] < MY_GLOBAL_SP (extras, myproc,
+                                                   stagger_index, vdim))
+      {
+        int npoints;
+
+        npoints = (MY_GLOBAL_SP (extras, myproc, stagger_index, vdim)
+                   - global_startpoint [vdim]) / downsample [vdim];
+        if ((MY_GLOBAL_SP (extras, myproc, stagger_index, vdim)
+             - global_startpoint [vdim]) % downsample [vdim])
+          npoints++;
+        my_global_startpoint [vdim] = global_startpoint [vdim] +
+                                      npoints*downsample [vdim];
+      }
+      else
+      {
+        my_global_startpoint [vdim] = global_startpoint [vdim];     
+      }
+    }
+    else
+      my_global_startpoint [vdim] = -1;
+  }
+
+  /* compute the local start- and endpoint from the global ranges */
+  totals = 1;
+  for (vdim = hdim = 0; vdim < vinfo.dim; vdim++)
+  {
+    stagger_index = CCTK_StaggerDirIndex (vdim, vinfo.stagtype);
+
+    if (my_global_startpoint [vdim] >= 0 &&
+        my_global_startpoint [vdim] <  MY_GLOBAL_EP (extras, myproc,
+                                                     stagger_index, vdim))
+    {
+      startpoint [vdim] = my_global_startpoint [vdim] -
+                          extras->lb [myproc][vdim];
+    }
+    else
+      startpoint [vdim] = -1;
+
+    if (global_endpoint [vdim] > MY_GLOBAL_SP (extras, myproc,
+                                               stagger_index, vdim))
+    {
+      endpoint [vdim] = MIN (MY_LOCAL_EP (extras, stagger_index, vdim),
+                             global_endpoint [vdim] - extras->lb[myproc][vdim]);
+    }
+    else
+      endpoint [vdim] = -1;
+
+#ifdef DEBUG
+    printf ("direction %d: local ranges [%d, %d)\n",
+            vdim, startpoint [vdim], endpoint [vdim]);
+#endif
+
+    if (endpoint [vdim] < 0 || startpoint [vdim] < 0)
+    {
+      totals = 0;
+      endpoint [vdim] = startpoint [vdim];
+    }
+
+    if (do_dir [vdim])
+    {
+      /* compute the global and local size in each hyperslab dimension */
+      hsize_global [hdim] = (global_endpoint[vdim] - global_startpoint[vdim]) /
+                            downsample [vdim];
+      if ((global_endpoint [vdim] - global_startpoint [vdim]) %
+          downsample [vdim])
+      {
+        hsize_global [hdim]++;
+      }
+      hsize [hdim] = (endpoint [vdim] - startpoint [vdim]) / downsample [vdim];
+      if ((endpoint [vdim] - startpoint [vdim]) % downsample [vdim])
+      {
+        hsize [hdim]++;
+      }
+      totals *= hsize [hdim];
+      hdim++;
+    }
+  }
+
+#ifdef DEBUG
+  printf ("total number of hyperslab data points: %d\n", totals);
+#endif
+
+  /* nested loop over vinfo.dim dimensions */
+  /* NOTE: the following code assumes startpoint [vdim] < endpoint [vdim] */
+  if (totals > 0)
+  {
+    void *vdata = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
+
+
+    /* if requested, compute the offsets into the global hyperslab */
+    if (hoffset_global)
+    {
+      for (vdim = hdim = 0; vdim < vinfo.dim; vdim++)
+      {
+        if (do_dir [vdim])
+        {
+          hoffset_global [hdim] = (my_global_startpoint [vdim] -
+                            global_startpoint [vdim]) / downsample [vdim];
+#ifdef DEBUG
+          printf ("hoffset_global, hsize in direction %d: %d, %d\n",
+                  hdim, hoffset_global [hdim], hsize [hdim]);
+#endif
+          hdim++;
+        }
+      }
+    }
+
+    /* allocate the buffer for the hyperslab data */
+    *hdata = malloc (totals * CCTK_VarTypeSize (vinfo.vartype));
+
+    /* compute the points_per_dim[] vector */
+    /* NOTE: this could be computed at startup and kept in a GH extension
+             once we have one for thorn Hyperslab */
+    points_per_dim [0] = 1;
+    for (vdim = 1; vdim < vinfo.dim; vdim++)
+      points_per_dim [vdim] = points_per_dim [vdim-1] * GH->cctk_lsh [vdim-1];
+
+    /* now get the hyperslab data points using that wonderful macro... */
+    switch (vinfo.vartype)
+    {
+      case CCTK_VARIABLE_CHAR:
+        PICKUP_HYPERSLAB_DATA (CCTK_CHAR, vinfo.dim, vdata, *hdata,
+                               point, startpoint, endpoint, downsample,
+                               points_per_dim);
+        break;
+      case CCTK_VARIABLE_INT:
+        PICKUP_HYPERSLAB_DATA (CCTK_INT, vinfo.dim, vdata, *hdata,
+                               point, startpoint, endpoint, downsample,
+                               points_per_dim);
+        break;
+      case CCTK_VARIABLE_REAL:
+        PICKUP_HYPERSLAB_DATA (CCTK_REAL, vinfo.dim, vdata, *hdata,
+                               point, startpoint, endpoint, downsample,
+                               points_per_dim);
+        break;
+      default:
+        CCTK_WARN (1, "Unsupported variable type");
+        retval = -1;
+        break;
+    }
+  }
+  else
+  {
+    /* clear out the return values if there's no hyperslab data */
+    *hdata = NULL;
+    memset (hsize, 0, hdim * sizeof (int));
+  }
+
+  /* free allocated temporary memory */
+  free (point);
+
+  return (retval);
+}
+
+
+/*@@
+  @routine    Hyperslab_GetHyperslab
+  @date       Fri May 12 2000
+  @author     Thomas Radke
+  @desc
+              Extract a hyperslab from a Cactus array variable.
+              All processor-local hyperslab data are collected
+              into a single global buffer.
+  @enddesc 
+  @calls
+  @calledby
+  @history 
+  @endhistory 
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        vindex
+  @vdesc      index of variable to get a hyperslab from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        vtimelvl
+  @vdesc      timelvl of variable to get a hyperslab from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        hdim
+  @vdesc      dimensionality of the requested hyperslab
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        global_startpoint
+  @vdesc      global coordinates of the hyperslab origin
+  @vtype      const int [dimensions of vindex]
+  @vio        in
+  @endvar
+  @var        directions
+  @vdesc      directions which span the hyperslab
+              This is the direction vector for 1D hyperslabs,
+              or the normal vector of the 2 lowest dimensions
+              for hyperslabs of higher dimensions.
+              If the hyperslab is of same dimensionality as the variable
+              the directions vector is ignored.
+  @vtype      const int [dimensions of vindex]
+  @vio        in
+  @endvar
+  @var        extents
+  @vdesc      number of grid points to follow in each hyperslab direction
+              starting from origin
+              Negative values are taken as extents up to the grid boundaries.
+  @vtype      const int [hdim]
+  @vio        in
+  @endvar
+  @var        downsample_
+  @vdesc      downsampling values for each hyperslab dimension
+  @vtype      const int [hdim]
+  @vio        in
+  @endvar
+  @var        hdata
+  @vdesc      pointer to store the address of the hyperslab data buffer
+  @vtype      void **
+  @vio        out
+  @endvar
+  @var        hsize
+  @vdesc      sizes of the hyperslab data buffer in each dimension
+  @vtype      int [hdim]
+  @vio        out
+  @endvar
+@@*/
+int Hyperslab_GetHyperslab (cGH *GH, int target_proc, int vindex, int vtimelvl,
+                            int hdim,
+                            const int global_startpoint [/*vdim*/],
+                            const int directions [/*vdim*/],
+                            const int extents [/*hdim*/],
+                            const int downsample_ [/*hdim*/],
+                            void **hdata,
+                            int hsize [/*hdim*/])
+{
+  int myproc, nprocs;         /* processor identification */
+  int retval;                 /* the return value */
+  cGroup vinfo;               /* variable's group info structure */
+  void *hdata_local;          /* buffer holding local hyperslab data */
+  void **hdata_ptr;           /* address of local hyperslab buffer pointer */
+  int *hsize_local,           /* sizes of the local and global hyperslab */
+      *hsize_global;          /* for each direction */
+  int *hoffset_local;         /* offsets of local into the global hyperslab */
+  const char *errormsg;       /* explanation string in case of errors */
+#ifdef CCTK_MPI
+  int i, proc;                /* loopers */
+  int totals_local,           /* number of local and global */
+      totals_global;          /* hyperslab points */
+  MPI_Comm comm;              /* MPI communicator to use */
+  MPI_Datatype mpi_vtype;     /* MPI datatype to use */
+  int vtypesize;              /* byte-size of a single data point */
+  void *chunked_hdata;        /* receive buffer for all hyperslab chunks */
+  int *displs, *recvcnts;     /* displacements/receive counts for gathering */
+  CCTK_INT *sizes_local,      /* MPI buffers for communicating the */
+           *sizes_global;     /* hyperslab sizes and offsets */
+#endif
+
+
+  /* do some plausibility checks */
+  errormsg = checkParameters (GH, vindex, vtimelvl, hdim,
+                              global_startpoint, directions, extents,
+                              downsample_, hdata, hsize);
+
+  /* FIXME: hack for getting diagonals from 3D variables
+            This is calling Gerd's CollectData1D() routine which can
+            extract non-axis-parallel lines too but is fixed to 3D data. */
+  if (errormsg && ! strcmp (errormsg, "Given normal vector isn't axis-parallel"))
+  {
+    int length;
+
+
+    /* get the info on the variable to extract a hyperslab from */
+    CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
+
+    if (hdim != 1 || vinfo.dim != 3)
+    {
+      CCTK_WARN (1, "Non-axis-parallel hyperslabs are supported as 1D lines "
+                    "from 3D variables only");
+      return (-1);
+    }
+
+    length = extents [0];
+    if (length > 0)
+      length--;
+    return (Hyperslab_CollectData1D (GH, vindex, vtimelvl, global_startpoint,
+                                     directions, downsample_ [0], length,
+                                     hdata, hsize, target_proc));
+  }
+
+  /* immediately return in case of errors */
+  if (errormsg)
+  {
+    CCTK_WARN (1, errormsg);
+    return (-1);
+  }
+
+  /* get processor information */
+  myproc = CCTK_MyProc (GH);
+  nprocs = CCTK_nProcs (GH);
+  if (target_proc >= nprocs)
+  {
+    CCTK_WARN (1, "Invalid target processor ID");
+    return (-1);
+  }
+
+  /* target processors must pass pointers to where to store the results */
+  if (target_proc < 0 || target_proc == myproc)
+  {
+    if (! hdata || ! hsize)
+    {
+      CCTK_WARN (1, "Must pass valid hyperslab data and sizes buffer pointers");
+      return (-1);
+    }
+
+    /* default return values: no hyperslab data */
+    *hdata = NULL;
+    memset (hsize, 0, hdim * sizeof (int));
+  }
+
+  /* get the info on the variable to extract a hyperslab from */
+  CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
+
+  /* set the pointers to pass to Hyperslab_GetLocalHyperslab() */
+  if (nprocs == 1)
+  {
+    hdata_ptr = hdata;
+    hsize_local = hsize_global = hsize;
+    hoffset_local = NULL;
+  }
+  else
+  {
+    hdata_ptr = &hdata_local;
+    hsize_local = (int *) malloc (3 * hdim * sizeof (int));
+    hoffset_local = hsize_local + hdim;
+    hsize_global = hoffset_local + hdim;
+  }
+
+  /* get the processor-local hyperslab */
+  retval = Hyperslab_GetLocalHyperslab (GH, vindex, vtimelvl, hdim,
+                                        global_startpoint, directions, extents,
+                                        downsample_, hdata_ptr, hsize_local,
+                                        hsize_global, hoffset_local);
+
+  /* that's it for the single processor case */
+  if (nprocs == 1)
+    return (retval);
+
+#ifdef CCTK_MPI
+  if (retval)
+  {
+    CCTK_WARN (1, "Hyperslab_GetLocalHyperslab() failed\n");
+    free (hsize_local);
+    return (retval);
+  }
+
+  /* now build the CCTK_INT buffer to communicate the sizes and offsets */
+  sizes_local = (CCTK_INT *) malloc ((1 + nprocs) * (2*hdim + 1)
+                                                  * sizeof (CCTK_INT));
+  sizes_global = sizes_local + 2*hdim + 1;
+  totals_local = 1;
+  for (i = 0; i < hdim; i++)
+  {
+    totals_local *= hsize_local [i];
+    sizes_local [0*hdim + i] = hsize_local [i];
+    sizes_local [1*hdim + i] = hoffset_local [i];
+  }
+  sizes_local [2*hdim + 0] = totals_local;
+
+  comm = PUGH_pGH (GH)->PUGH_COMM_WORLD;
+  /* FIXME: do an Allgather here so that all procs know how many data points
+            we have. This is useful to decide whether we have to gather 
+            something at all or not and thus can return immediately */
+  CACTUS_MPI_ERROR (MPI_Allgather (sizes_local, 2*hdim + 1, PUGH_MPI_INT,
+                    sizes_global, 2*hdim + 1, PUGH_MPI_INT, comm));
+
+  /* sum up the total number of hyperslab data points */
+  totals_global = sizes_global [2*hdim];
+  for (proc = 1; proc < nprocs; proc++) 
+    totals_global += sizes_global [proc*(2*hdim+1) + 2*hdim];
+
+  /* immediately return if there is no data at all */
+  if (totals_global <= 0)
+  {
+    free (sizes_local);
+    free (hsize_local);
+    return (retval);
+  }
+
+  vtypesize = CCTK_VarTypeSize (vinfo.vartype);
+  if (target_proc < 0 || target_proc == myproc)
+  {
+    /* copy dimensions for the returned hyperslab */
+    memcpy (hsize, hsize_global, hdim * sizeof (int));
+printf ("hsize_global = %d ", hsize [0]);
+for (i = 1; i < hdim; i++)
+  printf ("%d ", hsize [i]);
+printf ("\n");
+
+    /* set the displacements/receive counts for the gather operation */
+    displs = (int *) malloc (2 * nprocs * sizeof (int));
+    recvcnts = displs + nprocs;
+
+    /* for hdim == 1 we let MPI_Gatherv() do the sorting directly into the
+       hyperslab buffer using the given hyperslab offsets and sizes */
+    if (hdim == 1)
+    {
+      for (proc = 0; proc < nprocs; proc++)
+      {
+        displs [proc]   = sizes_global [proc*(2*hdim+1) + 2*hdim - 1];
+        recvcnts [proc] = sizes_global [proc*(2*hdim+1) + 2*hdim];
+      }
+    }
+    else
+    {
+      displs [0] = 0;
+      recvcnts [0] = sizes_global [2*hdim];
+      for (proc = 1; proc < nprocs; proc++)
+      {
+        displs [proc]   = displs [proc-1] + sizes_global [(proc-1)*(2*hdim+1) +
+                                                          2*hdim];
+        recvcnts [proc] = sizes_global [proc*(2*hdim+1) + 2*hdim];
+      }
+    }
+
+    chunked_hdata = malloc (totals_global * vtypesize);
+  }
+  else
+  {
+    displs = recvcnts = NULL;
+    chunked_hdata = NULL;
+  }
+
+  /* detect the MPI datatype to use */
+  mpi_vtype = 0;
+  switch (vinfo.vartype)
+  {
+    case CCTK_VARIABLE_CHAR:
+      mpi_vtype = PUGH_MPI_CHAR; break;
+    case CCTK_VARIABLE_INT:
+      mpi_vtype = PUGH_MPI_INT; break;
+    case CCTK_VARIABLE_REAL:
+      mpi_vtype = PUGH_MPI_REAL; break;
+    default:
+      CCTK_WARN (1, "Unsupported variable type"); break;
+  }
+
+  /* collect the hyperslab chunks from all processors */
+  if (target_proc < 0)
+  {
+    CACTUS_MPI_ERROR (MPI_Allgatherv (*hdata_ptr, totals_local, mpi_vtype,
+                      chunked_hdata, recvcnts, displs, mpi_vtype, comm));
+  }
+  else
+  {
+    CACTUS_MPI_ERROR (MPI_Gatherv (*hdata_ptr, totals_local, mpi_vtype,
+                      chunked_hdata, recvcnts, displs, mpi_vtype,
+                      target_proc, comm));
+  }
+
+  /* Now we got the global hyperslab data in a chunked array.
+     The user wants it unchunked, so let's sort it... */
+  if (target_proc < 0 || target_proc == myproc)
+  {
+    /* for hdim == 1 the hyperslab was already sorted by MPI_Gatherv()
+       we just need to return that buffer */
+    if (hdim == 1)
+    {
+      *hdata = chunked_hdata;
+    }
+    else
+    {
+      int j;
+      int *point;
+      int *points_per_hdim;
+      char *copy_from, *copy_to;
+      CCTK_INT *hsize_chunk, *hoffset_chunk;
+      int linear_hoffset;
+
+
+      /* allocate temporary vectors */
+      point = (int *) malloc (2 * hdim * sizeof (int));
+      points_per_hdim = point + hdim;
+
+      points_per_hdim [0] = 1;
+      for (i = 1; i < hdim; i++)
+        points_per_hdim [i] = points_per_hdim [i-1] * hsize [i-1];
+
+      /* allocate buffer for the returned hyperslab */
+      *hdata = malloc (totals_global * vtypesize);
+
+      /* use char pointers for easy incrementing when copying */
+      copy_from = (char *) chunked_hdata;
+      copy_to   = (char *) *hdata;
+
+      /* now copy the chunks from each processor into the global hyperslab */
+      for (proc = 0; proc < nprocs; proc++)
+      {
+        /* skip processors which didn't contribute any data */
+        if (sizes_global [proc * (2*hdim + 1) + 2*hdim] <= 0)
+          continue;
+
+        hsize_chunk = sizes_global + proc * (2*hdim+1);
+        hoffset_chunk = hsize_chunk + hdim;
+
+        /* set startpoint to zero */
+        memset (point, 0, hdim * sizeof (int));
+
+        i = 1;
+        while (1)
+        {
+          /* check for end of current loop */
+          if (point [i] >= hsize_chunk [i])
+          {
+            /* increment outermost loopers */
+            for (i++; i < hdim; i++)
+              if (++point [i] < hsize_chunk [i])
+                break;
+
+            /* done if beyond outermost loop */
+            if (i >= hdim)
+              break;
+
+            /* reset innermost loopers */
+            for (i--; i > 0; i--)
+              point [i] = 0;
+            i = 1;
+          }
+
+          /* copy the line */
+          for (linear_hoffset = 0, j = 1; j < hdim; j++)
+            linear_hoffset += (hoffset_chunk [j] + point [j]) *
+                               points_per_hdim [j];
+          memcpy (copy_to + linear_hoffset * vtypesize,
+                  copy_from, hsize_chunk [0] * vtypesize);
+          copy_from += hsize_chunk [0] * vtypesize;
+
+          /* increment current looper */
+          point [i]++;
+
+        } /* end of nested loops over all dimensions */
+      } /* end of loop over all processors */
+
+      /* free allocated temporary vectors and the chunk buffer */
+      free (chunked_hdata);
+      free (point);
+    }
+
+    /* free allocated vectors for MPI communication */
+    free (displs);
+  }
+
+  /* free allocated arrays for communicating the hyperslab sizes of offsets */
+  free (sizes_local);
+  free (hsize_local);
+
+#endif /* CCTK_MPI */
+
+  return (retval);
+}
+
+
+/********************** local routines ************************************/
+static const char *checkParameters (cGH *GH, int vindex, int vtimelvl,
+                                    int hdim,
+                                    const int global_startpoint [/*vdim*/],
+                                    const int directions [/*vdim*/],
+                                    const int extents [/*hdim*/],
+                                    const int downsample_ [/*hdim*/],
+                                    void **hdata,
+                                    int hsize [/*hdim*/])
+{
+  int dim;                       /* looper */
+  int num_directions;            /* number of non-zero directions */
+  cGroup vinfo;                  /* variable's group info */
+
+
+  /* check the variable index and timelevel */
+  if (vindex < 0 || vindex >= CCTK_NumVars ())
+    return ("Invalid variable index");
+  if (vtimelvl < 0 || vtimelvl >= CCTK_NumTimeLevelsFromVarI (vindex))
+    return ("Invalid timelevel");
+
+  /* check the passed pointers */
+  if (! global_startpoint || ! directions || ! extents || ! downsample_ ||
+      ! hdata || ! hsize)
+    return ("NULL pointer(s) passed as parameters");
+
+  /* check the extent and downsample parameters */
+  for (dim = 0; dim < hdim; dim++)
+  {
+    if (extents [dim] == 0)
+      return ("Invalid hyperslab extent parameters");
+    if (downsample_ [dim] <= 0)
+      return ( "Invalid hyperslab downsample parameters");
+  }
+
+  /* get the info on the variable to extract a hyperslab from */
+  if (CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo) < 0)
+    return ("Couldn't get group info");
+
+  /* check the variable's grouptype */
+  if (vinfo.grouptype != CCTK_GF && vinfo.grouptype != CCTK_ARRAY)
+    return ("Invalid variable group type");
+
+  /* check the hyperslab dimension */
+  if (hdim < 0 || hdim > vinfo.dim)
+    return ("Invalid hyperslab dimension");
+
+  /* check the direction(s) of the hyperslab
+     if it is less-dimensional than the variable */
+  if (hdim != vinfo.dim)
+  {
+    for (dim = 0, num_directions = 0; dim < vinfo.dim; dim++)
+      if (directions [dim])
+        num_directions++;
+    if (num_directions == 0)
+      return ("Given normal vector is a null vector");
+    if (num_directions > 1)
+      return ("Given normal vector isn't axis-parallel");
+  }
+
+  /* check if PUGH is active */
+  if (! PUGH_pGH (GH))
+    return ("No GH extension for PUGH found. Did you activate thorn PUGH ?");
+
+  return (NULL);
+}