Added (still incomplete) implementation of the new Hyperslab API.

Not really usable yet, only single-processor case works. More changes to follow.


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

6 files changed, 1120 insertions, 58 deletions

diff --git a/src/CollectData1D.c b/src/CollectData1D.c
index ea5a052..38b5f60 100644
--- a/src/CollectData1D.c
+++ b/src/CollectData1D.c
@@ -29,16 +29,16 @@ CCTK_FILEVERSION(CactusPUGH_PUGHSlab_CollectData1D_c)
           }
 
 /* local function prototypes */
-int GetLocalLine(cGH *GH, pGExtras *extras,
+int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl,
+                            const int *origin, const int *directions,
+                            int downsample, int length,
+                            void **hdata, int *hsize, int proc);
+int GetLocalLine(const cGH *GH, pGExtras *extras,
                  int vardim, const int *S_l, const int *u_l, int ds, int len_l,
                  int *locstart, int *locend, int *nlocpoints);
-#ifdef CCTK_MPI
-static MPI_Datatype CCTKtoMPItype(pGH *pughGH, int cctktype);
-#endif
-
 
 
-static int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl,
+static int CollectLocalData1D(const cGH *GH, int vindex, int vtimelvl,
                               const int *origin, const int *directions,
                               int downsample, int length,
                               void **hdata, int *hsize)
@@ -110,7 +110,7 @@ static int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl,
 }
 
 
-int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl,
+int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl,
                             const int *origin, const int *directions,
                             int downsample, int length,
                             void **hdata, int *hsize, int proc)
@@ -126,7 +126,7 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl,
 
 #ifdef CCTK_MPI
   int iproc;
-  pGH *pughGH;   
+  const pGH *pughGH;   
   CCTK_INT  tmp, *rem_dcount;
   MPI_Datatype vmpitype;
   int *recvcnts, *displs;
@@ -157,7 +157,7 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl,
   rem_dcount = NULL;
   recvcnts = NULL;
   displs = NULL;
-  vmpitype= CCTKtoMPItype(pughGH, vtype);
+  vmpitype = PUGH_MPIDataType (pughGH, vtype);
 
   /* Send local number of data elements */
   tmp = (CCTK_INT) loc_dcount;
@@ -236,41 +236,5 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl,
     *hsize = all_dcount;
   }
 
-#ifdef DEBUGME
-  printf(" CollectData1D done \n");
-#endif
-  return(0);
+  return (0);
 }
-
-
-#ifdef CCTK_MPI
-static MPI_Datatype CCTKtoMPItype(pGH *pughGH, int cctktype)
-{
-  MPI_Datatype mpitype;
-  switch (cctktype) {
-
-  case CCTK_VARIABLE_CHAR:
-      mpitype  = PUGH_MPI_CHAR;
-      break;
-
-  case CCTK_VARIABLE_INT:
-    mpitype = PUGH_MPI_INT;
-    break;
-
-  case CCTK_VARIABLE_REAL:
-    mpitype = PUGH_MPI_REAL;
-    break;
-
-  case CCTK_VARIABLE_COMPLEX:
-    mpitype = pughGH->PUGH_mpi_complex;
-    break;
-
-  default:
-    CCTK_WARN (1, "Unsupported MPI variable type");
-    mpitype = 0;
-    break;
-  }
-
-  return(mpitype);
-}
-#endif
diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c
new file mode 100644
index 0000000..0c883cc
--- /dev/null
+++ b/src/GetHyperslab.c
@@ -0,0 +1,573 @@
+ /*@@
+   @file      GetHyperslab.c
+   @date      Sun 2 Dec 2001
+   @author    Thomas Radke
+   @desc
+              Routines to extract hyperslabs from CCTK array variables
+   @enddesc
+   @version   $Id$
+ @@*/
+
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "PUGHSlab.h"
+#include "PUGHSlabi.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHSlab_GetHyperslab_c)
+
+/********************************************************************
+ ********************    Macro Definitions   ************************
+ ********************************************************************/
+/* define this if you want debugging output */
+/* #define DEBUG 1 */
+
+
+/********************************************************************
+ ********************    Internal Routines   ************************
+ ********************************************************************/
+static int GetLocalHyperslab (const cGH *GH,
+                              const hslab_mapping_t *mapping,
+                              int vindex,
+                              int timelevel,
+                              int hdatatype,
+                              void *hdata);
+static const char *checkParameters (const cGH *GH,
+                                    const hslab_mapping_t *mapping,
+                                    int vindex,
+                                    int timelevel,
+                                    void *hdata);
+/********************************************************************
+ ********************    External Routines   ************************
+ ********************************************************************/
+/* Gerd's routine to get 1D lines from a 3D array */
+int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl,
+                             const int *origin,
+                             const int *directions,
+                             int downsampling,
+                             int length,
+                             void **hdata,
+                             int *hsize,
+                             int proc);
+
+
+
+CCTK_INT Hyperslab_Get (const cGH *GH,
+                        CCTK_INT   mapping_handle,
+                        CCTK_INT   vindex,
+                        CCTK_INT   timelevel,
+                        CCTK_INT   hdatatype,
+                        void      *hdata)
+{
+  int retval;
+  hslab_mapping_t *mapping;
+
+
+  mapping = PUGHSlabi_GetMapping (mapping_handle);
+  if (mapping == NULL)
+  {
+    return (-1);
+  }
+
+  /* check mapping consistency */
+  /*** FIXME ***/
+
+  /* get the processor-local hyperslab */
+  retval = GetLocalHyperslab (GH, mapping, vindex, timelevel, hdatatype, hdata);
+
+  return (retval);
+}
+
+
+CCTK_INT Hyperslab_GetList (const cGH *GH,
+                            CCTK_INT mapping_handle,
+                            CCTK_INT num_arrays,
+                            const CCTK_INT *vindices   /* num_arrays */,
+                            const CCTK_INT *timelevels /* num_arrays */,
+                            const CCTK_INT *hdatatypes /* num_arrays */,
+                            void *const *hdata         /* num_arrays */)
+{
+  int i, retval;
+
+
+  retval = 0;
+  for (i = 0; i < num_arrays; i++)
+  {
+    if (Hyperslab_Get (GH, mapping_handle, vindices[i],
+                       timelevels ? timelevels[i] : 0,
+                       hdatatypes ? hdatatypes[i] : -1,
+                       hdata[i]) == 0)
+    {
+      retval++;
+    }
+  }
+
+  return (retval);
+}
+
+/********************************************************************
+ ********************    Internal Routines   ************************
+ ********************************************************************/
+/*@@
+  @routine    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      PUGHSlab_GetDatatypeConversionFn
+
+  @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        timelevel
+  @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        hdatatype
+  @vdesc      CCTK datatype of the requested hyperslab
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        conversion_fn
+  @vdesc      pointer to a user-supplied data conversion function
+  @vtype      t_hslabConversionFn
+  @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
+  @vtype      const int[hdim times 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        free_data
+  @vdesc      address of flag which decides whether the returned data needs
+              to be freed or not
+  @vtype      int *
+  @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
+@@*/
+static int GetLocalHyperslab (const cGH *GH,
+                              const hslab_mapping_t *mapping,
+                              int vindex,
+                              int timelevel,
+                              int hdatatype,
+                              void *hdata)
+{
+  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 *points_per_dim;           /* points per subvolume */
+  int myproc;                    /* local processor ID */
+  int i;                         /* general looper */
+  int vdim;                      /* looper over all source dimensions */
+  int vdata_size,                /* size of one data point in bytes for */
+      hdata_size;                /* source and hyperslab data */
+  int dim0_points;               /* number of hyperslab points in dim 0 */
+  int dim0_hsize;                /* byte size of hyperslab points in dim 0 */
+  const char *typed_vdata;             /* byte pointers into source and */
+  char *typed_hdata;             /* hyperslab data arrays */
+  const void *vdata;
+  int retval;                    /* the return value (0 for success) */
+  cGroup vinfo;                  /* variable's group info */
+  pGH *pughGH;                   /* pointer to the current pGH */
+  pGA *GA;                       /* the variable's GA structure from PUGH */
+  const char *errormsg;          /* error message string */
+  t_hslabConversionFn conversion_fn;
+
+
+  /* do some plausibility checks */
+  errormsg = checkParameters (GH, mapping, vindex, timelevel, hdata);
+
+  /* immediately return in case of errors */
+  if (errormsg)
+  {
+    CCTK_WARN (1, errormsg);
+    return (-1);
+  }
+
+  /* check if there's any data to extract */
+  if (mapping->totals == 0)
+  {
+    return (0);
+  }
+
+  /* get the info on the variable to extract a hyperslab from */
+  CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
+
+  /* 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 (mapping->is_diagonal_in_3D)
+  {
+    const int origin[] = {0, 0, 0};
+    const int directions[] = {1, 1, 1};
+    int dummy_hsize;
+    void *dummy_hdata;
+
+
+    retval = Hyperslab_CollectData1D (GH, vindex, timelevel, origin, directions,
+                                      1, -1, &dummy_hdata, &dummy_hsize,
+                                      mapping->target_proc);
+    if (! retval)
+    {
+      memcpy (hdata, dummy_hdata,
+              mapping->totals * CCTK_VarTypeSize (vinfo.vartype));
+      free (dummy_hdata);
+    }
+
+    return (retval);
+  }
+
+  /* if datatype conversion was requested
+     get the appropriate predefined datatype conversion routine
+     in case the user didn't supply one by her own */
+  if (hdatatype < 0)
+  {
+    hdatatype = vinfo.vartype;
+  }
+  conversion_fn = mapping->conversion_fn;
+  if (vinfo.vartype != hdatatype)
+  {
+    if (conversion_fn == NULL)
+    {
+      conversion_fn = PUGHSlabi_GetDatatypeConversionFn (vinfo.vartype,
+                                                         hdatatype);
+      if (! conversion_fn)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "No predefined PUGHSlab routine available to convert "
+                    "'%s' into '%s'", CCTK_VarTypeName (vinfo.vartype),
+                    CCTK_VarTypeName (hdatatype));
+        return (-1);
+      }
+    }
+  }
+  else if (conversion_fn)
+  {
+    CCTK_WARN (8, "Datatype conversion routine supplied but no datatype "
+                  "conversion requested. Ignoring conversion routine...");
+    conversion_fn = NULL;
+  }
+
+  /* allocate the temporary arrays */
+  point = (int *) malloc (5 * vinfo.dim * sizeof (int));
+  startpoint           = point + 1*vinfo.dim;
+  endpoint             = point + 2*vinfo.dim;
+  downsample           = point + 3*vinfo.dim;
+  points_per_dim       = point + 4*vinfo.dim;
+
+  memcpy (startpoint, mapping->local_startpoint, vinfo.dim * sizeof (int));
+  memcpy (endpoint, mapping->local_endpoint, vinfo.dim * sizeof (int));
+  memcpy (downsample, mapping->downsample, vinfo.dim * sizeof (int));
+
+  /* get the pGH pointer and the variable's GA structure */
+  pughGH = PUGH_pGH (GH);
+  GA     = (pGA *) pughGH->variables[vindex][timelevel];
+
+  /* get the local processor ID */
+  myproc = CCTK_MyProc (GH);
+
+  /* nested loop over vinfo.dim dimensions */
+  /* NOTE: the following code assumes startpoint[vdim] < endpoint[vdim] */
+  vdata = CCTK_VarDataPtrI (GH, timelevel, vindex);
+
+  if (mapping->full_hyperslab && conversion_fn == NULL)
+  {
+    memcpy (hdata, vdata, mapping->totals * CCTK_VarTypeSize (vinfo.vartype));
+  }
+  else
+  {
+    /* get the byte size of a single data point
+       in the variable and hyperslab data array */
+    vdata_size = CCTK_VarTypeSize (vinfo.vartype);
+    hdata_size = CCTK_VarTypeSize (hdatatype);
+
+    typed_hdata = (char *) hdata;
+
+    /* 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 PUGHSlab */
+    points_per_dim[0] = 1;
+    for (vdim = 1; vdim < vinfo.dim; vdim++)
+    {
+      points_per_dim[vdim] = points_per_dim[vdim-1] *
+                             GA->extras->lnsize[vdim-1];
+    }
+
+    /* get the number of hyperslab points in lowest dimension
+       and their size in bytes */
+    dim0_points = (endpoint[0] - startpoint[0]) / downsample[0];
+    if ((endpoint[0] - startpoint[0]) % downsample[0])
+    {
+      dim0_points++;
+    }
+    dim0_hsize = dim0_points * hdata_size;
+
+    /* transform the ranges into byte ranges */
+    for (i = 0; i < vinfo.dim; i++)
+    {
+      startpoint[i] *= vdata_size;
+      endpoint[i]   *= vdata_size;
+      downsample[i] *= vdata_size;
+    }
+
+    /* initialize the index vector to the local startpoint */
+    memcpy (point, startpoint, vinfo.dim * sizeof (point[0]));
+
+    /* do the nested loops starting with the innermost */
+    vdim = 1;
+    while (1)
+    {
+      /* check for end of current loop */
+      if (vinfo.dim > 1 && point[vdim] >= endpoint[vdim])
+      {
+        /* increment outermost loopers */
+        for (vdim++; vdim < vinfo.dim; vdim++)
+        {
+          point[vdim] += downsample[vdim];
+          if (point[vdim] < endpoint[vdim])
+          {
+            break;
+          }
+        }
+
+        /* done if beyond outermost loop */
+        if (vdim >= vinfo.dim)
+        {
+          break;
+        }
+
+        /* reset innermost loopers */
+        for (vdim--; vdim > 0; vdim--)
+        {
+          point[vdim] = startpoint[vdim];
+        }
+        vdim = 1;
+      }
+
+      /* get the byte pointer into the source array */
+      typed_vdata = (const char *) vdata + point[0];
+#if 0
+fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_FullName (vindex));
+#endif
+      for (i = 1; i < vinfo.dim; i++)
+      {
+        typed_vdata += point[i] * points_per_dim[i];
+      }
+
+      /* copy the data in lowest dimension: if possible copy all data points
+         in a row otherwise do it one by one */
+      if (downsample[0] == vdata_size)
+      {
+        if (mapping->conversion_fn)
+        {
+          mapping->conversion_fn (typed_vdata, typed_hdata, dim0_points, 1,1);
+        }
+        else
+        {
+#if 0
+fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vdata, typed_hdata);
+#endif
+          memcpy (typed_hdata, typed_vdata, dim0_hsize);
+        }
+      }
+      else
+      {
+        if (mapping->conversion_fn)
+        {
+          mapping->conversion_fn (typed_vdata, typed_hdata, dim0_points,
+                                  downsample[0], 1);
+          typed_vdata += downsample[0] * dim0_points;
+        }
+        else
+        {
+          for (i = 0; i < dim0_hsize; i += hdata_size)
+          {
+            memcpy (typed_hdata + i, typed_vdata, vdata_size);
+            typed_vdata += downsample[0];
+          }
+        }
+      }
+      typed_hdata += dim0_hsize;
+
+      if (vinfo.dim > 1)
+      {
+        /* increment current looper */
+        point[vdim] += downsample[vdim];
+      }
+      else
+      {
+        /* exit loop if hyperslab dim is only 1D */
+        break;
+      }
+
+    } /* end of nested loops over all dimensions */
+  } /* end of branch extracting the hyperslab data */
+
+  /* free allocated temporary memory */
+  free (point);
+
+  return (0);
+}
+
+
+static const char *checkParameters (const cGH *GH,
+                                    const hslab_mapping_t *mapping,
+                                    int vindex,
+                                    int timelevel,
+                                    void *hdata)
+{
+  cGroup vinfo;                  /* variable's group info */
+#if 0
+  int i, vdim;                   /* looper */
+  int num_directions;            /* number of non-zero directions */
+#endif
+
+
+  /* check the variable index and timelevel */
+  if (vindex < 0 || vindex >= CCTK_NumVars ())
+  {
+    return ("Invalid variable index");
+  }
+  if (timelevel < 0 || timelevel >= CCTK_NumTimeLevelsFromVarI (vindex))
+  {
+    return ("Invalid timelevel");
+  }
+
+#if 0
+  /* check the passed pointers */
+  if (! mapping->global_origin || ! directions || ! extents || ! mapping->downsample ||
+      ! hdata || ! hsize)
+  {
+    return ("NULL pointer(s) passed as parameters");
+  }
+
+  /* check the extent and downsample parameters */
+  for (vdim = 0; vdim < hdim; vdim++)
+  {
+    if (extents[vdim] == 0)
+    {
+      return ("Invalid hyperslab extent parameters");
+    }
+    if (mapping->downsample[vdim] <= 0)
+    {
+      return ( "Invalid hyperslab downsample parameters");
+    }
+  }
+#endif
+
+  /* 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 (mapping->hdim <= 0 || mapping->hdim > vinfo.dim)
+  {
+    return ("Invalid hyperslab dimension");
+  }
+
+#if 0
+  /* check the direction(s) of the hyperslab */
+  for (i = 0; i < mapping->hdim; i++)
+  {
+    for (vdim = 0, num_directions = 0; vdim < vinfo.dim; vdim++)
+    {
+      if (directions[i * vinfo.dim + vdim])
+      {
+        num_directions++;
+      }
+    }
+    if (num_directions == 0)
+    {
+      return ("Given direction vector is a null vector");
+    }
+    if (num_directions != 1)
+    {
+      return ("Given direction vector isn't orthogonal");
+    }
+  }
+#endif
+
+  /* check if PUGH is active */
+  if (! PUGH_pGH (GH))
+  {
+    return ("No GH extension for PUGH found. Did you activate thorn PUGH ?");
+  }
+
+  return (NULL);
+}
diff --git a/src/Mapping.c b/src/Mapping.c
new file mode 100644
index 0000000..b53b050
--- /dev/null
+++ b/src/Mapping.c
@@ -0,0 +1,519 @@
+ /*@@
+   @file      Mapping.c
+   @date      Sun 2 Dec 2001
+   @author    Thomas Radke
+   @desc
+              Routines to define hyperslab mappings
+   @enddesc
+   @version   $Id$
+ @@*/
+
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "cctk.h"
+
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+
+#include "PUGHSlab.h"
+#include "PUGHSlabi.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHSlab_Mapping_c)
+
+/********************************************************************
+ ********************    Macro Definitions   ************************
+ ********************************************************************/
+/* define this if you want 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))
+
+
+/********************************************************************
+ ********************    Internal Routines   ************************
+ ********************************************************************/
+static int nmapping_list = 0;
+static hslab_mapping_t *mapping_list = NULL;
+
+static int checkFullHyperslab (const pGA *GA,
+                               const CCTK_INT *origin,
+                               const CCTK_INT *extent,
+                               const hslab_mapping_t *mapping);
+
+
+CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
+           const cGH *GH,
+           CCTK_INT vindex,
+           CCTK_INT hdim,
+           const CCTK_INT *direction  /* vdim*hdim */,
+           const CCTK_INT *origin     /* vdim */,
+           const CCTK_INT *extent     /* vdim */,
+           const CCTK_INT *downsample /* hdim */,
+           CCTK_INT table_handle,
+           CCTK_INT target_proc,
+           t_hslabConversionFn conversion_fn,
+           CCTK_INT *hsize            /* hdim */)
+{
+  int i, j, vdim, num_dirs, stagger_type, retval;
+  int stagger_index;
+  int myproc;
+  hslab_mapping_t *mapping;
+  const char *error_msg;
+  pGH *pughGH;                   /* pointer to the current pGH */
+  pGA *GA;                       /* the variable's GA structure from PUGH */
+
+
+  /* PUGHSlab doesn't use table information */
+  if (table_handle >= 0)
+  {
+    CCTK_WARN (1, "Hyperslab_DefineGlobalMappingByIndex: table information is "
+                  "ignored");
+  }
+
+  /* check parameter consistency */
+  retval = 0;
+  error_msg = NULL;
+  vdim = CCTK_GroupDimFromVarI (vindex);
+  if (vdim <= 0)
+  {
+    error_msg = "invalid variable index given";
+    retval = -1;
+  }
+  else if (hdim < 0 || hdim > vdim)
+  {
+    error_msg = "invalid hyperslab dimension given";
+    retval = -2;
+  }
+  else if (hsize == NULL)
+  {
+    error_msg = "NULL pointer passed for hyperslab size return parameter";
+    retval = -3;
+  }
+  else if (target_proc >= CCTK_nProcs (GH))
+  {
+    error_msg = "invalid target procesor ID given";
+    retval = -4;
+  }
+  else
+  {
+    for (i = 0; i < vdim; i++)
+    {
+      retval |= origin[i] < 0 || extent[i] <= 0;
+      if (downsample && i < hdim)
+      {
+        retval |= downsample[i] <= 0;
+      }
+    }
+    if (retval)
+    {
+      error_msg = "invalid hyperslab origin/extent/downsample given";
+      retval = -5;
+    }
+  }
+  if (! retval)
+  {
+    mapping = (hslab_mapping_t *) malloc (sizeof (hslab_mapping_t));
+    if (mapping)
+    {
+      mapping->vectors = (int *) malloc ((6*vdim + (2+vdim)*hdim) * sizeof (int));
+    }
+    if (mapping == NULL || mapping->vectors == NULL)
+    {
+      if (mapping)
+      {
+        free (mapping);
+      }
+      error_msg = "couldn't allocate hyperslab mapping structure";
+      retval = -6;
+    }
+  }
+
+  /* return in case of errors */
+  if (retval)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Hyperslab_DefineGlobalMappingByIndex: %s", error_msg);
+    return (retval);
+  }
+
+  /* assign memory for the other vectors */
+  mapping->local_startpoint  = mapping->vectors + 0*vdim;
+  mapping->local_endpoint    = mapping->vectors + 1*vdim;
+  mapping->global_startpoint = mapping->vectors + 2*vdim;
+  mapping->global_endpoint   = mapping->vectors + 3*vdim;
+  mapping->do_dir            = mapping->vectors + 4*vdim;
+  mapping->downsample        = mapping->vectors + 5*vdim;
+  mapping->local_hsize       = mapping->vectors + 6*vdim + 0*hdim;
+  mapping->global_hsize      = mapping->vectors + 6*vdim + 1*hdim;
+  mapping->directions        = mapping->vectors + 6*vdim + 2*hdim;
+
+  /* check direction vectors */
+  for (j = 0; j < hdim; j++)
+  {
+    num_dirs = 0;
+    for (i = 0; i < vdim; i++)
+    {
+      if (direction[j*vdim + i])
+      {
+        mapping->directions[i] = 1;
+        num_dirs++;
+      }
+    }
+    mapping->is_diagonal_in_3D = num_dirs == 3 && hdim == 1;
+    if (num_dirs != 1 && ! mapping->is_diagonal_in_3D)
+    {
+      free (mapping->vectors);
+      free (mapping);
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Hyperslab_DefineGlobalMappingByIndex: %d-direction vector "
+                  "isn't axis-orthogonal", j);
+      return (-7);
+    }
+  }
+  for (i = 0; i < vdim; i++)
+  {
+    mapping->do_dir[i] = 0;
+    for (j = 0; j < hdim; j++)
+    {
+      if (direction[j*vdim + i])
+      {
+        mapping->do_dir[i]++;
+      }
+    }
+    if (mapping->do_dir[i] > 1)
+    {
+      free (mapping->vectors);
+      free (mapping);
+      CCTK_WARN (1, "Hyperslab_DefineGlobalMappingByIndex: duplicate direction "
+                    "vectors given");
+      return (-8);
+    }
+  }
+
+  /* get the pGH pointer and the variable's GA structure */
+  pughGH = PUGH_pGH (GH);
+  GA     = (pGA *) pughGH->variables[vindex][0];
+  myproc = CCTK_MyProc (GH);
+  stagger_type = CCTK_GroupStaggerIndexGI (CCTK_GroupIndexFromVarI (vindex));
+
+  /* check extent */
+  for (i = j = 0; i < vdim; i++)
+  {
+    if (mapping->do_dir[i])
+    {
+      if (origin[i] + extent[j] > GA->extras->nsize[i])
+      {
+        free (mapping->vectors);
+        free (mapping);
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Hyperslab_DefineGlobalMappingByIndex: extent in "
+                    "%d-direction exceeds grid size", j);
+        return (-8);
+      }
+      j++;
+    }
+  }
+
+  /* now fill out the hyperslab mapping structure */
+  for (i = j = 0; i < vdim; i++)
+  {
+    mapping->downsample[i] = 1;
+
+    if (mapping->do_dir[i])
+    {
+      if (downsample)
+      {
+        mapping->downsample[i] = downsample[j];
+      }
+      mapping->global_hsize[j] = extent[j] / mapping->downsample[i];
+      if (extent[j] % mapping->downsample[i])
+      {
+        mapping->global_hsize[j]++;
+      }
+      /* subtract ghostzones for periodic BC */
+      if (GA->connectivity->perme[i])
+      {
+        mapping->global_hsize[j] -= 2 * GA->extras->nghostzones[i];
+      }
+      j++;
+    }
+  }
+
+  /* check whether the full local data patch was requested as hyperslab */
+  mapping->full_hyperslab = checkFullHyperslab (GA, origin, extent, mapping);
+  if (mapping->full_hyperslab)
+  {
+    memset (mapping->local_startpoint, 0,                  vdim * sizeof (int));
+    memcpy (mapping->local_endpoint,   GA->extras->lnsize, vdim * sizeof (int));
+    mapping->totals = GA->extras->npoints;
+  }
+  else
+  {
+    /* compute the global endpoint */
+    for (i = j = 0; i < vdim; i++)
+    {
+      mapping->global_endpoint[i] = origin[i] +
+                                    mapping->do_dir[i] ? extent[j++] : 1;
+    }
+
+    /* compute this processor's global startpoint from the global ranges */
+    for (i = 0; i < vdim; i++)
+    {
+      stagger_index = CCTK_StaggerDirIndex (i, stagger_type);
+
+      if (origin[i] < MY_GLOBAL_EP (GA->extras, myproc, stagger_index, i))
+      {
+        mapping->global_startpoint[i] = origin[i];
+        if (origin[i] < MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i))
+        {
+          int npoints;
+
+          npoints = (MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i)
+                     - origin[i]) / mapping->downsample[i];
+          if ((MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i)
+               - origin[i]) % mapping->downsample[i])
+          {
+            npoints++;
+          }
+          mapping->global_startpoint[i] += npoints * mapping->downsample[i];
+        }
+      }
+      else
+      {
+        mapping->global_startpoint[i] = -1;
+      }
+    }
+
+    /* compute the local start- and endpoint from the global ranges */
+    mapping->totals = 1;
+    for (i = j = 0; i < vdim; i++)
+    {
+      stagger_index = CCTK_StaggerDirIndex (i, stagger_type);
+
+      if (mapping->global_startpoint[i] >= 0 &&
+          mapping->global_startpoint[i] <  MY_GLOBAL_EP (GA->extras, myproc,
+                                                         stagger_index, i))
+      {
+        mapping->local_startpoint[i] = mapping->global_startpoint[i] -
+                                       GA->extras->lb[myproc][i];
+      }
+      else
+      {
+        mapping->local_startpoint[i] = -1;
+      }
+
+      if (mapping->global_endpoint[i] > MY_GLOBAL_SP (GA->extras, myproc,
+                                                      stagger_index, i))
+      {
+        mapping->local_endpoint[i] = MIN (MY_LOCAL_EP (GA->extras, stagger_index, i),
+                                          mapping->global_endpoint[i] - GA->extras->lb[myproc][i]);
+      }
+      else
+      {
+        mapping->local_endpoint[i] = -1;
+      }
+
+#ifdef DEBUG
+      printf ("direction %d: local ranges[%d, %d)\n",
+              i, mapping->local_startpoint[i], mapping->local_endpoint[i]);
+#endif
+
+      if (mapping->local_endpoint[i] < 0 || mapping->local_startpoint[i] < 0)
+      {
+        mapping->totals = 0;
+        mapping->local_endpoint[i] = mapping->local_startpoint[i];
+      }
+
+      if (mapping->do_dir[i])
+      {
+        /* compute the local size in each hyperslab dimension */
+        mapping->local_hsize[j] = (mapping->local_endpoint[i] -
+                                   mapping->local_startpoint[i]) /
+                                  mapping->downsample[i];
+        if ((mapping->local_endpoint[i] - mapping->local_startpoint[i]) %
+            mapping->downsample[i])
+        {
+          mapping->local_hsize[j]++;
+        }
+        mapping->totals *= mapping->local_hsize[j];
+        j++;
+      }
+    }
+  } /* end of else branch for 'if (mapping->full_hyperslab)' */
+
+  /* diagonals in 3D are special: global_hsize is minimum of grid extents */
+  if (mapping->is_diagonal_in_3D)
+  {
+    mapping->global_hsize[0] = GH->cctk_gsh[0];
+    if (mapping->global_hsize[0] > GH->cctk_gsh[1])
+    {
+      mapping->global_hsize[0] = GH->cctk_gsh[1];
+    }
+    if (mapping->global_hsize[0] > GH->cctk_gsh[2])
+    {
+      mapping->global_hsize[0] = GH->cctk_gsh[2];
+    }
+    mapping->totals = mapping->global_hsize[0];
+  }
+
+#ifdef DEBUG
+  printf ("total number of hyperslab data points: %d\n", mapping->totals);
+#endif
+
+  mapping->hdim = hdim;
+  mapping->vdim = vdim;
+  mapping->target_proc = target_proc;
+  mapping->conversion_fn = conversion_fn;
+
+  /* add this mapping to the mapping list */
+  if (mapping_list)
+  {
+    mapping_list->prev = mapping;
+  }
+  mapping->prev = NULL;
+  mapping->next = mapping_list;
+  mapping_list = mapping;
+
+  mapping->handle = nmapping_list++;
+
+  /* set the global hsize in the return arguments */
+  if (hsize)
+  {
+    for (j = 0; j < hdim; j++)
+    {
+      hsize[j] = mapping->global_hsize[j];
+    }
+  }
+
+#if 0
+  if (mapping->totals > 0)
+  {
+    /* if requested, compute the offsets into the global hyperslab */
+    if (hoffset_global)
+    {
+      for (i = j = 0; i < vdim; i++)
+      {
+        if (mapping->do_dir[i])
+        {
+          if (mapping->full_hyperslab)
+          {
+            hoffset_global[j] = GA->extras->lb[myproc][i];
+          }
+          else
+          {
+            hoffset_global[j] = (mapping->global_startpoint[i] -
+                                    origin[i]) / mapping->downsample[i];
+            if (GA->connectivity->perme[i])
+            {
+              hoffset_global[j] -= GA->extras->nghostzones[i];
+            }
+          }
+#ifdef DEBUG
+          printf ("hoffset_global, hsize in direction %d: %d, %d\n",
+                  j, hoffset_global[j], mapping->local_hsize[j]);
+#endif
+          j++;
+        }
+      }
+    }
+  }
+#endif
+
+  return (mapping->handle);
+}
+
+
+int Hyperslab_FreeMapping (CCTK_INT mapping_handle)
+{
+  hslab_mapping_t *mapping;
+
+
+  mapping = mapping_list;
+  while (mapping)
+  {
+    if (mapping->handle == mapping_handle)
+    {
+      if (mapping == mapping_list)
+      {
+        mapping_list = mapping_list->next;
+      }
+      else
+      {
+        if (mapping->next)
+        {
+          mapping->next->prev = mapping->prev;
+        }
+        if (mapping->prev)
+        {
+          mapping->prev->next = mapping->next;
+        }
+      }
+      free (mapping->vectors);
+      free (mapping);
+      return (0);
+    }
+    mapping = mapping->next;
+  }
+
+  return (-1);
+}
+
+
+hslab_mapping_t *PUGHSlabi_GetMapping (int mapping_handle)
+{
+  hslab_mapping_t *mapping;
+
+
+  mapping = mapping_list;
+  while (mapping && mapping->handle != mapping_handle)
+  {
+    mapping = mapping->next;
+  }
+
+  return (mapping);
+}
+
+
+static int checkFullHyperslab (const pGA *GA,
+                               const CCTK_INT *origin,
+                               const CCTK_INT *extent,
+                               const hslab_mapping_t *mapping)
+{
+  int i;
+  int is_full_hyperslab;
+
+
+  is_full_hyperslab = mapping->hdim == GA->extras->dim;
+
+  if (is_full_hyperslab)
+  {
+    for (i = 0; i < GA->extras->dim; i++)
+    {
+      is_full_hyperslab &= (origin[i] == 0);
+      is_full_hyperslab &= (extent[i] <= 0);
+      is_full_hyperslab &= (mapping->downsample[i] <= 1);
+      is_full_hyperslab &= (GA->connectivity->perme[i] == 0);
+    }
+  }
+
+  return (is_full_hyperslab);
+}
diff --git a/src/PUGHSlab.h b/src/PUGHSlab.h
index c2ff990..718edc6 100644
--- a/src/PUGHSlab.h
+++ b/src/PUGHSlab.h
@@ -25,12 +25,12 @@ typedef void (*t_hslabConversionFn) (const void *src,
                                      CCTK_INT dst_stride);
 
 /* function prototypes */
-CCTK_INT Hyperslab_Get (const cGH     *GH,
-                        CCTK_INT       mapping_handle,
-                        const CCTK_INT vindex,
-                        const CCTK_INT timelevel,
-                        const CCTK_INT hdatatype,
-                        void          *hdata);
+CCTK_INT Hyperslab_Get (const cGH *GH,
+                        CCTK_INT   mapping_handle,
+                        CCTK_INT   vindex,
+                        CCTK_INT   timelevel,
+                        CCTK_INT   hdatatype,
+                        void      *hdata);
 CCTK_INT Hyperslab_GetList (const cGH *GH,
                             CCTK_INT mapping_handle,
                             CCTK_INT num_arrays,
diff --git a/src/PUGHSlabi.h b/src/PUGHSlabi.h
index f251544..47f8907 100644
--- a/src/PUGHSlabi.h
+++ b/src/PUGHSlabi.h
@@ -26,15 +26,21 @@ typedef struct hslab_mapping_t
   int target_proc;
   int hdim;
   int vdim;
-  int *global_origin;                    /* vdim */
-  int *directions;                       /* hdim * vdim */
-  int *extents;                          /* hdim */
-  int *downsample;                       /* hdim */
-  int *global_hsize;                     /* hdim */
+  int *vectors;
+  int *local_startpoint;                 /* vdim */
+  int *local_endpoint;                   /* vdim */
+  int *global_startpoint;                /* vdim */
+  int *global_endpoint;                  /* vdim */
   int *do_dir;                           /* vdim */
+  int *downsample;                       /* vdim */
+  int *global_hsize;                     /* hdim */
+  int *local_hsize;                      /* hdim */
+  int *directions;                       /* hdim * vdim */
   int is_diagonal_in_3D;
+  unsigned int totals;
   t_hslabConversionFn conversion_fn;
   struct hslab_mapping_t *prev, *next;
+  int full_hyperslab;
 } hslab_mapping_t;
 
 
diff --git a/src/make.code.defn b/src/make.code.defn
index 2463bb7..7d322e8 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,4 +2,4 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c
+SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c