Usable version of thorn Hyperslab which extracts M-dimensional orthogonal

subspaces from an N-dimensional data space.

The Hyperslab API consists of 2 routines:
  - Hyperslab_GetLocalHyperslab() which returns a hyperslab
    from the processor's local data.
    This will probably be used for IO methods that want to output
    the hyperslab chunks in parallel (as is done currently by
    IOHDF5 and IOFlexIO).
  - Hyperslab_GetHyperslab() which collects hyperslab chunks from all
    processors using MPI and returns the global hyperslab data.
    This is used in IOASCII where the output is done by only one processor.

Open questions:
  - How to specify the directions in which to span the hyperslab.
  - Putting the parameters passed to Hyperslab routines into structures
    rather than passing them individually.


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

11 files changed, 197 insertions, 820 deletions

diff --git a/src/CollectData1D.c b/src/CollectData1D.c
index 3605423..9645fea 100644
--- a/src/CollectData1D.c
+++ b/src/CollectData1D.c
@@ -1,64 +1,156 @@
 
 #include <stdio.h>
 #include <stdlib.h>
-#include <math.h>  
 
 #include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "Hyperslabi.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "Hyperslab.h"
+
+/*#define DEBUGME  1*/
+
+
+#define PICKUP_LINE_DATA(GH, cctk_type, vsize, start, directions, vdata, hdata)\
+          {                                                                   \
+            int i, dim, point [3];                                            \
+            cctk_type *typed_data = (cctk_type *) hdata;                      \
+                                                                              \
+                                                                              \
+            for (i = 0; i < vsize; i++)                                       \
+            {                                                                 \
+              for (dim = 0; dim < 3; dim++)                                   \
+                point [dim] = start [dim] + i*directions [dim];               \
+              *typed_data++ = ((cctk_type *) vdata)                           \
+                [CCTK_GFINDEX3D (GH, point [0], point [1], point [2])];       \
+            }                                                                 \
+          }
+
+/* local function prototypes */
+int GetLocalLine(cGH *GH,
+                 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(int cctktype);
+#endif
 
 
-int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
-                  int *S_l, int *u_l, int *ds, int len_l, 
-                  void **dptr, int *dsize, int proc)
+
+static int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl,
+                              const int *origin, const int *directions,
+                              int downsample, int length,
+                              void **hdata, int *hsize)
 {
+  int idim, ierr;
+  cGroup groupinfo;
+  int locstart [3], locend [3];
+  void *vdata = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
 
-  void *loc_dptr;
-  int   loc_dcount[1];
 
+  CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &groupinfo);
+
+  if (groupinfo.dim != 3)
+  {
+    CCTK_WARN (1, "CollectLocalData1D() processes 3D variables only");
+    return (-1);
+  }
+
+  ierr = GetLocalLine(GH, 3, origin, directions, downsample, length, 
+                          locstart, locend, hsize);
+  if (ierr < 0)
+  {
+    CCTK_WARN (1, "GetLocalLine() failed");
+    return (ierr);
+  }
+
+  /* Transform from global to local index coordinate */
+  for(idim=0;idim<3;idim++)
+  {  
+    locstart[idim]-=GH->cctk_lbnd[idim];
+    locend[idim]  -=GH->cctk_lbnd[idim];
+  }
+
+#ifdef DEBUGME
+  printf("CollectLocal1D: %d npoints   lb: %d %d %d\n",
+         *hsize,GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2]);
+#endif
+
+  *hdata= malloc (*hsize * CCTK_VarTypeSize (groupinfo.vartype));
+
+  /* CCTK_VARIABLE_INT */
+  switch (groupinfo.vartype)
+  {
+    case CCTK_VARIABLE_CHAR:
+      PICKUP_LINE_DATA(GH, CCTK_CHAR, *hsize, locstart, directions, vdata, *hdata);
+      break;
+
+    case CCTK_VARIABLE_INT:
+      PICKUP_LINE_DATA(GH, CCTK_INT, *hsize, locstart, directions, vdata, *hdata);
+      break;
+
+    case CCTK_VARIABLE_REAL:
+      PICKUP_LINE_DATA(GH, CCTK_REAL, *hsize, locstart, directions, vdata, *hdata);
+      break;
+
+    default:
+      CCTK_WARN (1, "Unsupported variable type");
+      ierr = -1;
+      break;
+  }
+
+  return (ierr);
+}
+
+
+int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl,
+                            const int *origin, const int *directions,
+                            int downsample, int length,
+                            void **hdata, int *hsize, int proc)
+{
+
+  void *loc_dptr;
+  int   loc_dcount;
   void *all_dptr;
   int   all_dcount=0;
-
   int vtype, vtypesize;
-  int iproc,myproc,nprocs;
+  int myproc,nprocs;
+  int ierr;
 
 #ifdef CCTK_MPI
+  int iproc;
   pGH *pughGH;   
   CCTK_INT  tmp, *rem_dcount;
-
-  MPI_Datatype vmpi_type;
+  MPI_Datatype vmpitype;
   int *recvcnts, *displs;
 #endif
 
-  int idim, i,j,k;
-  int ierr;
-
 
   vtype    = CCTK_VarTypeI(vindex);
   vtypesize= CCTK_VarTypeSize(vtype);
 
   myproc   = CCTK_MyProc(GH); 
   nprocs   = CCTK_nProcs(GH);
-  vmpi_type= CCTKtoMPItype(vtype);
 
-  printf("Calling CollectLocal ... ");  
 
+  ierr = CollectLocalData1D(GH, vindex, vtimelvl,
+                            origin, directions, downsample, length, 
+                            &loc_dptr, &loc_dcount);
 
-  ierr = CollectLocalData1D(GH, vindex, vtimelvl, vdim, 
-                            S_l, u_l, ds, len_l, 
-                            &loc_dptr, loc_dcount);
+  if (ierr)
+  {
+    CCTK_WARN (1, "CollectLocalData1D() failed");
+    return (ierr);
+  }
 
 #ifdef CCTK_MPI
 
+  rem_dcount = NULL;
+  recvcnts = NULL;
+  displs = NULL;
+  vmpitype= CCTKtoMPItype(vtype);
+
   pughGH   = PUGH_pGH(GH);
 
   /* Send local number of data elements */
-  tmp = (CCTK_INT) loc_dcount [0];
+  tmp = (CCTK_INT) loc_dcount;
 
   if (proc < 0 || proc == myproc)
   {
@@ -89,11 +181,13 @@ int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
       all_dcount   += recvcnts[iproc];
     }
 
-    printf ("  dsize: %d\n", dsize [0]);
+#ifdef DEBUGME
+    printf ("  local hsize: %d\n", loc_dcount);
     for (iproc=0;iproc<nprocs;iproc++)
       printf(" DISPLS[%d] rem: %d   displs: %d   all %d \n",
              iproc, recvcnts[iproc],displs[iproc],all_dcount);
     fflush(stdout);
+#endif
 
     all_dptr = malloc (all_dcount * vtypesize);
   }
@@ -104,12 +198,12 @@ int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
 
   if (proc < 0)
   {
-    CACTUS_MPI_ERROR (MPI_Allgatherv(loc_dptr, loc_dcount [0], vmpi_type, 
-                                     all_dptr, recvcnts, displs, vmpi_type,
+    CACTUS_MPI_ERROR (MPI_Allgatherv(loc_dptr, loc_dcount, vmpitype, 
+                                     all_dptr, recvcnts, displs, vmpitype,
                                      pughGH->PUGH_COMM_WORLD));
   } else {
-    CACTUS_MPI_ERROR (MPI_Gatherv(loc_dptr, loc_dcount [0], vmpi_type, 
-                                  all_dptr, recvcnts, displs, vmpi_type,
+    CACTUS_MPI_ERROR (MPI_Gatherv(loc_dptr, loc_dcount, vmpitype, 
+                                  all_dptr, recvcnts, displs, vmpitype,
                                   proc, pughGH->PUGH_COMM_WORLD));
   }
 
@@ -128,137 +222,41 @@ int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
   /* assign return values */
   if (proc < 0 || proc == myproc)
   {
-    *dptr = all_dptr;
-    dsize[0] = all_dcount;
+    *hdata = all_dptr;
+    *hsize = all_dcount;
   }
 
+#ifdef DEBUGME
   printf(" CollectData1D done \n");
+#endif
   return(0);
 }
 
 
-int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, int vdim, 
-                       int *S_l, int *u_l, int *ds, int len_l,
-                       void **loc_dptr, int *dsize)
+#ifdef CCTK_MPI
+static MPI_Datatype CCTKtoMPItype(int cctktype)
 {
-  int *locstart, *locend, nlocpoints;
-  int ipnt, gridpnt[MAX_DIM];
-  int ierr;
-
-  int vtype,vtypesize;
-  int idim, lincount;
-
-  void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-
-
-  if (vdim > MAX_DIM)
-    return (-1);
+  MPI_Datatype mpitype;
+  switch (cctktype) {
 
-  vtype      = CCTK_VarTypeI(vindex);
-  vtypesize  = CCTK_VarTypeSize(vtype);
+  case CCTK_VARIABLE_CHAR:
+      mpitype  = PUGH_MPI_CHAR;
+      break;
 
-  locstart = (int*)malloc(vdim*sizeof(int));
-  locend   = (int*)malloc(vdim*sizeof(int));
+  case CCTK_VARIABLE_INT:
+    mpitype = PUGH_MPI_INT;
+    break;
 
-  ierr  = GetLocalLine(GH, vdim, S_l, u_l, ds[0], len_l, 
-                          locstart, locend, &nlocpoints);
+  case CCTK_VARIABLE_REAL:
+    mpitype = PUGH_MPI_REAL;
+    break;
 
-  if (ierr<0)
-  {
-    CCTK_WARN (1, "GetLocalLine() failed");
-    return (ierr);
+  default:
+    CCTK_WARN (1, "Unsupported MPI variable type");
+    mpitype = 0;
+    break;
   }
 
-  /* Transform from global to local index coordinate */
-  for(idim=0;idim<vdim;idim++)
-  {  
-    locstart[idim]-=GH->cctk_lbnd[idim];
-    locend[idim]  -=GH->cctk_lbnd[idim];
-  }
-
-  dsize[0]= nlocpoints;
-
-  printf("CollectLocal1D: %d npoints   lb: %d %d %d\n",
-         nlocpoints,GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2]);
-
-  *loc_dptr= malloc(nlocpoints * vtypesize);
-
-  lincount=0; 
-
-  /* zero out index array so that unused dimensions are 0 */
-  memset (gridpnt, 0, sizeof (gridpnt));
-
-  /* CCTK_VARIABLE_INT */
-  if (vtype == CCTK_VARIABLE_INT)
-  {
-    CCTK_INT *int_dptr = (CCTK_INT *) *loc_dptr;
-
-    for (ipnt=0; ipnt<nlocpoints; ipnt++)
-    { 
-      for (idim=0;idim<vdim;idim++)
-        gridpnt[idim] = locstart[idim]+ipnt * u_l[idim];
-      int_dptr[lincount++]=((CCTK_INT *) data)
-        [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-    }
-  }   
-  /* CCTK_VARIABLE_CHAR */
-  else if (vtype == CCTK_VARIABLE_CHAR)
-  {
-    CCTK_CHAR *char_dptr = (CCTK_CHAR *) *loc_dptr;
-
-    for (ipnt=0; ipnt<nlocpoints; ipnt++)
-    { 
-      for (idim=0;idim<vdim;idim++)
-        gridpnt[idim] = locstart[idim]+ipnt * u_l[idim];
-      char_dptr[lincount++]=((CCTK_CHAR *) data)
-        [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-    }
-  }
- 
-  /* CCTK_VARIABLE_REAL */
-  else  if (vtype == CCTK_VARIABLE_REAL)
-  {
-    CCTK_REAL *real_dptr = (CCTK_REAL *) *loc_dptr;
-
-    for (ipnt=0; ipnt<nlocpoints; ipnt++)
-    { 
-      printf(" nlocpoints: %d \n",nlocpoints);
-      for (idim=0;idim<vdim;idim++) 
-        gridpnt[idim] = locstart[idim]+ipnt * u_l[idim];
-      printf(" gridpnts: ( %d %d %d ) ", 
-             gridpnt[0],gridpnt[1],gridpnt[2]);
-        
-      real_dptr[lincount++]=((CCTK_REAL *) data)
-        [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-      printf(" Value:  %f \n",real_dptr[lincount-1]);
-    }
-  }   
-  /* CCTK_VARIABLE_REAL4 */
-  else if (vtype == CCTK_VARIABLE_REAL4)
-  {    
-    CCTK_REAL4 *real4_dptr = (CCTK_REAL4*) *loc_dptr;
-
-    for (ipnt=0; ipnt<nlocpoints; ipnt++)
-    { 
-      for (idim=0;idim<vdim;idim++)
-        gridpnt[idim] = locstart[idim]+ipnt * u_l[idim];
-      
-      real4_dptr[lincount++]=((CCTK_REAL4 *) data)
-        [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-    }
-  }
-  else
-  {
-    char *fullname = CCTK_FullName (vindex);
-
-    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "Unsupported type for variable '%s'", fullname);
-    free (fullname);
-    ierr = -1;
-  }
-
-  free(locstart);
-  free(locend);
-
-  return (ierr);
+  return(mpitype);
 }
+#endif
diff --git a/src/CollectData2D.c b/src/CollectData2D.c
deleted file mode 100644
index 3019a1e..0000000
--- a/src/CollectData2D.c
+++ /dev/null
@@ -1,115 +0,0 @@
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "Hyperslabi.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#define MAX(a,b) ((a)>(b) ? (a) : (b))
-
-int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		  int *S_s, int *dir, int *ds, int *len,
-		  void **dptr, int *dsize)
-{
-  return(0);
-}
-
-
- /*@@
-   @routine    CollectLocalDataND
-   @date       Fri May 12 17:51:30 2000
-   @author     Gerd Lanfermann
-   @desc 
-   S_s       :  starting point of surface
-   linedir[2]:  the two directions that span the surface: 0~x, 1~y, ...
-   ds[2]     :  downsampling in each dir, 
-   len[2]    :  number of (downsampled) point
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-
-int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		       int *S_s, int *linedir, int *ds, int *len,
-		       void **loc_dptr, int *dirsize,int *totsize)
-{
-
-  int *locstart, *locend, nlocpoints;
-  int ipnt[MAX_DIM], gridpnt[MAX_DIM];
-
-  int vtype,vtypesize;
-  int idim, idir, lincount;
-  int ierr;
- 
-  void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-
-  if (vdim > MAX_DIM)
-    return (-1);
-
-  vtype      = CCTK_VarTypeI(vindex);
-  vtypesize  = CCTK_VarTypeSize(vtype);
-
-  locstart = (int*)malloc(vdim*sizeof(int));
-  locend   = (int*)malloc(vdim*sizeof(int));
-
-  ierr = GetLocalSlab(GH, vdim, sdim, 
-		      S_s, linedir, ds, len, 
-		      locstart, locend, dirsize, &totsize);
-  if (ierr<0)
-  {
-    CCTK_WARN (1, "GetLocalSlab() failed");
-    return (ierr);
-  }
-
-  /* Transform from global to local index coordinate */
-  for(idim=0;idim<vdim;idim++)
-  {  
-    locstart[idim]-=GH->cctk_lbnd[idim];
-    locend[idim]  -=GH->cctk_lbnd[idim];
-  }
-
-  *loc_dptr= malloc(nlocpoints * vtypesize);
-  lincount=0; 
-
-  /* zero out index array so that unused dimensions are 0 */
-  memset (gridpnt, 0, sizeof (gridpnt));
-  memset (ipnt   , 0, sizeof (ipnt));
-
-    /* CCTK_VARIABLE_INT */
-  if (vtype == CCTK_VARIABLE_REAL)
-  {
-    CCTK_INT *int_dptr = (CCTK_REAL *) *loc_dptr;
-
-    for (ipnt[0]=0; ipnt[0]<dirsize[0] && (0 < sdim); ipnt[0]++)
-      for (ipnt[1]=0; ipnt[1]<dirsize[1] && (1 < sdim); ipnt[1]++)
-	for (ipnt[2]=0; ipnt[2]<dirsize[2] && (2 < sdim); ipnt[2]++)
-	{ 
-	  gridpnt[idim] = locstart[idim]; 
-	  for (idim=0;idim<vdim;idim++)
-	    for (idir = 0; idir < sdim; idir ++)
-	      gridpnt[idim] += ipnt[idir] * u_l[idir][idim];
-
-	}
-    int_dptr[lincount++]=((CCTK_INT *) data)
-      [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-    }
-  }
-}
-   
-  
-  
-  
diff --git a/src/CollectDataND.c b/src/CollectDataND.c
deleted file mode 100644
index f70e8e3..0000000
--- a/src/CollectDataND.c
+++ /dev/null
@@ -1,135 +0,0 @@
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "Hyperslabi.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-
-/* SMART returns the loop boundary: 
-   1 if dir>sdim (the hardcoded index coordinate is less than dimension of slab) 
-   sz[idir] otherwise (size of the slab in direction <dir>.
-
-   This way we make sure that we enter the loop of that particular direction
-   at least one time, even though there is no size associated along that direction.
-*/
-#define SMART(sz,dir,sd) ((dir)>(sd) ? 1 : (sz)[(dir)])
-
-int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		  int *S_s, int *dir, int *ds, int *len,
-		  void **dptr, int *dsize)
-{
- return(0);
-}
-
-
- /*@@
-   @routine    CollectLocalDataND
-   @date       Fri May 12 17:51:30 2000
-   @author     Gerd Lanfermann
-   @desc 
-   S_s       :  starting point of surface
-   linedir[2]:  the two directions that span the surface: 0~x, 1~y, ...
-   ds[2]     :  downsampling in each dir, 
-   len[2]    :  number of (downsampled) point
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-
-int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		       int *S_s, int *linedir, int *ds, int *len,
-		       void **loc_dptr, int *dirsize, int *totsize)
-{
-
-  int *locstart, *locend, nlocpoints;
-  int ipnt[MAX_DIM], gridpnt[MAX_DIM];
-
-  int vtype,vtypesize;
-  int idim, idir, lincount;
-  int ierr=0;
- 
-  void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-
-  if (vdim > MAX_DIM)
-    return (-1);
-
-  vtype      = CCTK_VarTypeI(vindex);
-  vtypesize  = CCTK_VarTypeSize(vtype);
-
-  locstart = (int*)malloc(vdim*sizeof(int));
-  locend   = (int*)malloc(vdim*sizeof(int));
-  dirsize  = (int*)malloc(sdim*sizeof(int));
-
-  printf("Calling GetLocalSlab \n");
-
-  ierr = GetLocalSlab(GH, vdim, sdim, 
-		      S_s, linedir, ds, len, 
-		      locstart, locend, dirsize, &nlocpoints);
-  if (ierr<0)
-  {
-    CCTK_WARN (1, "GetLocalSlab() failed");
-    return (ierr);
-  }
-
-  /* Transform from global to local index coordinate */
-  for(idim=0;idim<vdim;idim++)
-  {  
-    locstart[idim]-=GH->cctk_lbnd[idim];
-    locend[idim]  -=GH->cctk_lbnd[idim];
-  }
-
-  *loc_dptr= malloc(nlocpoints * vtypesize);
-  lincount=0; 
-
-  /* zero out index array so that unused dimensions are 0 */
-  memset (gridpnt, 0, sizeof (gridpnt));
-  memset (ipnt   , 0, sizeof (ipnt));
-
-  /* CCTK_VARIABLE_INT */
-  gridpnt[0] = locstart[0]; 
-  gridpnt[1] = locstart[1]; 
-  gridpnt[2] = locstart[2]; 
-
-  if (vtype == CCTK_VARIABLE_REAL)
-  {
-    CCTK_REAL *real_dptr = (CCTK_REAL *) *loc_dptr;
-    
-    for (ipnt[2]=0 ; ipnt[2]<SMART(dirsize, 2, sdim); ipnt[2]++) 
-
-      gridpnt[linedir[0]] = locstart[linedir[0]] + ipnt[0] * ds[0];
-
-      for (ipnt[1]=0; ipnt[1]<SMART(dirsize, 1, sdim); ipnt[1]++) {
-	
-	gridpnt[linedir[1]] = locstart[linedir[1]] + ipnt[1] * ds[1];
-	
-	for (ipnt[0]=0; ipnt[0]<SMART(dirsize, 0, sdim); ipnt[0]++) {
-	  
-	  { 
-	    gridpnt[linedir[2]] = locstart[linedir[2]] + ipnt[2] * ds[2];    
-
-	    real_dptr[lincount++]=((CCTK_REAL *) data)
-	      [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
-	  }
-      }
-    }
-
-  }
-}
-   
-  
-  
-  
diff --git a/src/GetLocalLine.c b/src/GetLocalLine.c
index b105628..19c6008 100644
--- a/src/GetLocalLine.c
+++ b/src/GetLocalLine.c
@@ -3,18 +3,25 @@
 #include <math.h>  
 
 #include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
 
-#include "Hyperslab.h"
-#include "Hyperslabi.h"
+#define MAX_DIM 3
+#define MAX_FACE 6
+#define BAD -42
+#define EPS 1e-10
+
+#define ABS(a) ((a)<0 ? -(a) : (a))
+#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0))
+#define LIN_INDEX3D(GH,P) ((P[0])+ (GH)->cctk_gsh[0]*((P[1]) + (GH)->cctk_gsh[1]*(P[2])))
+
+#define MIN(a,b) ((a)<(b) ? (a) : (b))
+#define MAX(a,b) ((a)>(b) ? (a) : (b))
 
 /*$#define LINE_DEBUG$*/
 
-int IsPointOnProc_INT (int dim, int *gownership, int *pnt);
+int IsPointOnProc_INT (int dim, int *gownership, const int *pnt);
 int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt);
 
-int IsPointOnProc_INT(int dim, int *gownership, int *pnt)
+int IsPointOnProc_INT(int dim, int *gownership, const int *pnt)
 {
   int retval=1; /*assume YES */
   int idim;
@@ -43,14 +50,14 @@ int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt)
 }
 
 int GetLocalLine(cGH *GH, 
-		 int vardim, int *S_l, int *u_l, int ds, int len_l, 
+		 int vardim, const int *S_l, const int *u_l, int ds, int len_l, 
 		 int *locstart, int *locend, int *nlocpoints)
 {
 
   const int n_s[MAX_DIM][MAX_FACE] = 
-  { -1, 1, 0, 0, 0, 0,
-     0, 0,-1, 1, 0, 0,
-     0, 0, 0, 0,-1, 1};
+  {{ -1, 1, 0, 0, 0, 0},
+   { 0, 0,-1, 1, 0, 0},
+   { 0, 0, 0, 0,-1, 1}};
 
 
   /* variables for proc. surfaces */
@@ -58,20 +65,17 @@ int GetLocalLine(cGH *GH,
   int P_s[MAX_DIM][MAX_FACE];
   int gownership[MAX_FACE];
 
-  int E_l[MAX_DIM];
-  int notcoll[MAX_FACE], isloc[2];
+  int notcoll[MAX_FACE];
   CCTK_REAL t_pnt1[MAX_DIM];
 
   int prochasline=0;
 
-  int tmp1,tmp2,t;
+  int tmp1,tmp2;
   CCTK_REAL rtmp1;
 
-  int idim, iface, myproc, linsec;
+  int idim, iface, linsec;
 
 
-  myproc = CCTK_MyProc(GH);
-  
 #ifdef LINE_DEBUG
   printf("\n+++++++++ GET LOCAL LINE \n");
   printf("S_l [XYZ] : %d %d %d \n",
diff --git a/src/GetLocalSlab.c b/src/GetLocalSlab.c
deleted file mode 100644
index 67c559c..0000000
--- a/src/GetLocalSlab.c
+++ /dev/null
@@ -1,155 +0,0 @@
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "Hyperslabi.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-#define SLAB_DEBUG
-
-
- /*@@
-   @routine    GetLocalSlab.c
-   @date       Fri May 12 17:51:30 2000
-   @author     Gerd Lanfermann
-   @desc 
-   S_s       :  starting point of surface
-   linedir[2]:  the two directions that span the surface: 0~x, 1~y, ...
-   ds[2]     :  downsampling in each dir, 
-   len[2]    :  number of (downsampled) point
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-
-int GetLocalSlab(cGH *GH, int vdim, int sdim,
-		 int *S_s, int *linedir, int *ds, int *len,
-		 int *locstart, int *locend, int *dirsize, int *totsize)
-{
-
-  
-
-  int tstart[MAX_DIM];     /* temp start points */  
-  int tlocstart[MAX_DIM]; 
-  int tlocend[MAX_DIM];
-
-  int u_l[MAX_DIM][MAX_DIM];
-
-  int vtype,vtypesize;
-  int idim, idir;
-  int ierr;
-  int phasplane=1;
-
-  /* Initialize the local start/end to BAD values and hope for the best... */
-  for (idim=0; idim<vdim; idim++) 
-  {
-    tstart[idim]   = BAD;
-    locstart[idim] = BAD;
-    locend[idim]   = BAD;
-
-    /* line vectors 1..3 are 000, 000 */
-    u_l[0][idim]  = 0;
-    u_l[1][idim]  = 0;
-    u_l[2][idim]  = 0;
-
-    /* make sure that all procs above and below the plane, 
-       know that they don't have a plane */
-    if ((idim!=linedir[0]) && (idim!=linedir[1]))
-      if((GH->cctk_lbnd[idim]<S_s[idim]) && 
-	 (GH->cctk_ubnd[idim]>S_s[idim])) 
-      {
-	locstart[idim] = S_s[idim];
-	locend[idim]   = S_s[idim];
-	phasplane = 1;
-      }
-      else phasplane = 0;
-  }
-
-  /* construct line vectors that span the surface:
-     e.g. linedir[0,1] = 1,2 means YZ surface:  
-     u_l[0][]=(000) -> (010)
-     u_l[1][]=(000) -> (001) */
-  for (idir=0; idir<sdim; idir++)
-  {
-    u_l[idir][linedir[idir]]=1;
-  }
-
-  /* init to 1 (neutral element for multiplication) */
-  totsize[0] = 1;
-  
-  for (idir=0;idir<2;idir++) 
-  {
-#ifdef SLAB_DEBUG
-    printf(" testing arm %d:  S_s (%d)  < ubnd (%d) ",
-	   idir, S_s[linedir[idir]], GH->cctk_ubnd[linedir[idir]]);
-#endif
-    if (S_s[linedir[idir]]<GH->cctk_ubnd[linedir[idir]])
-    {
-#ifdef SLAB_DEBUG
-      printf(" .... YES \n");
-#endif
-      /* Initialize startpoints */
-      for (idim=0; idim<vdim; idim++) 
-      {
-	tstart[idim]   = GH->cctk_lbnd[idim]+GH->cctk_nghostzones[idim];
-      }
-      tstart[linedir[idir]]=S_s[linedir[idir]];
-      ierr = GetLocalLine(GH, vdim, tstart, u_l[idir], ds[idir], len[idir], 
-			  tlocstart, tlocend, &dirsize[idir]);
-      locstart[linedir[idir]]= tlocstart[linedir[idir]];
-      locend  [linedir[idir]]= tlocend[linedir[idir]];
-
-      totsize[0] *=dirsize[idir];
-
-      if (locstart[linedir[idir]]==BAD) phasplane &=0;
-    }
-    else 
-    {
-#ifdef SLAB_DEBUG
-      printf(" .... NO \n");
-#endif
-      phasplane &= 0;
-      totsize[0] = 0;
-    }
-    
-  }
- 
-  if (phasplane)
-  {
-#ifdef SLAB_DEBUG
-    printf("LOC START:  %d %d %d   %d %d %d   -> %d\n",
-	   locstart[0],locstart[1],locstart[2],
-	   locend[0],locend[1],locend[2],
-	   totsize[0]);
-#endif
-  }
-  else 
-  {
-#ifdef SLAB_DEBUG
-    printf("proc has no plane\n");
-#endif
-  }
-}
-   
-  
-  
-  
-  
-  /*$printf("LB/UB [%d , %d, %d]  - [%d, %d, %d] \n",
-	 GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2],
-	 GH->cctk_ubnd[0],GH->cctk_ubnd[1],GH->cctk_ubnd[2]);$*/
-  
- 
diff --git a/src/Hyperslab.h b/src/Hyperslab.h
index 7376341..cfc41f0 100644
--- a/src/Hyperslab.h
+++ b/src/Hyperslab.h
@@ -1,34 +1,29 @@
-
-
-int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
-		  int *S_l, int *u_l, int *ds, int len_l, 
-		  void **dptr, int *dsize, int proc);
-
-int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, int vdim, 
-		       int *S_l, int *u_l, int *ds, int len_l,
-		       void **loc_dptr, int *dsize);
-
-int CollectData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
-		  int *S_s, int dir[2], int ds[2], int len[2],
-		  void **dptr, int *dsize);
-
-int CollectLocalData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
-		       int *S_s, int dir[2], int ds[2], int len[2],
-		       void **dptr, int *idx1_ptr, int *idx2_ptr, int *dsize);
-
-int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		  int *S_s, int *dir, int *ds, int *len,
-		  void **dptr, int *dsize);
-
-int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
-		       int *S_s, int *dir, int *ds, int *len,
-		       void **dptr,  int *dirsize,int *totsize);
-int GetLocalLine(cGH *GH, 
-		 int vardim, int *S_l, int *u_l, int ds, int len_l, 
-		 int *locstart, int *locend, int *nlocpoints);
-
-int GetLocalSlab(cGH *GH, int vdim, int sdim,
-		 int *S_s, int *linedir, int *ds, int *len,
-		 int *locstart, int *locend, int *dirsize, int *totsize);
-
-int IsOnGlobGrid(cGH *GH, int vdim, int *grdpnt);
+ /*@@
+   @header    Hyperslab.h
+   @date      Sun 28 May 2000
+   @author    Thomas Radke
+   @desc 
+              Function declarations of thorn Hyperslab
+   @enddesc
+   @history
+   @endhistory
+ @@*/
+
+
+/* function prototypes */
+int Hyperslab_GetLocalHyperslab (cGH *GH, int vindex, int vtimelvl,
+                                 int hdim,
+                                 const int global_startpoint [/*vdim*/],
+                                 const int directions [/*vdim*/],
+                                 const int lengths [/*hdim*/],
+                                 const int downsample_ [/*hdim*/],
+                                 void **hdata,
+                                 int hsize [/*hdim*/], int ghsize [/*hdim*/],
+                                 int hoffset [/*hdim*/]);
+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 lengths [/*hdim*/],
+                            const int downsample_ [/*hdim*/],
+                            void **hdata, int hsize [/*hdim*/]);
diff --git a/src/Hyperslabi.h b/src/Hyperslabi.h
deleted file mode 100644
index 2adae72..0000000
--- a/src/Hyperslabi.h
+++ /dev/null
@@ -1,15 +0,0 @@
-
-
-#define MAX_DIM 3
-#define MAX_FACE 6
-#define BAD -42
-#define EPS 1e-10
-
-#define ABS(a) ((a)<0 ? -(a) : (a))
-#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0))
-#define LIN_INDEX3D(GH,P) ((P[0])+ (GH)->cctk_gsh[0]*((P[1]) + (GH)->cctk_gsh[1]*(P[2])))
-
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#define MAX(a,b) ((a)>(b) ? (a) : (b))
-
-
diff --git a/src/Misc.c b/src/Misc.c
deleted file mode 100644
index 7c6747a..0000000
--- a/src/Misc.c
+++ /dev/null
@@ -1,44 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-
-
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-
-int IsOnGlobGrid(cGH *GH, int vdim, int *grdpnt)
-{
-  int idim, retval=0;
-  for (idim=0;idim<vdim;idim++)
-    if ((grdpnt[idim]<0) || (grdpnt[idim]>GH->cctk_gsh[idim]))
-      retval-=1;
-  return(retval);
-}
-      
-
-MPI_Datatype CCTKtoMPItype(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;
-
-  default:
-    printf("Unsupported MPI variable type in Hyperslab\n");
-    mpitype = -1;
-  }
-  
-  return(mpitype);
-}
diff --git a/src/TestSlab.c b/src/TestSlab.c
deleted file mode 100644
index 9f48e3f..0000000
--- a/src/TestSlab.c
+++ /dev/null
@@ -1,85 +0,0 @@
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
-
-#define MAX_DIM 3
-#define MAX_FACE 6
-#define EPS 1e-10
-#define BAD -42
-
-#define ABS(a) ((a)<0 ? -(a) : (a))
-
-
-void InitTestGF(CCTK_ARGUMENTS) 
-{  
-  DECLARE_CCTK_ARGUMENTS
-  int i,j,k;
-
-
-  for(i=0; i < cctk_lsh[0]; i++)
-  {
-    for(j=0; j < cctk_lsh[1]; j++)
-    {
-      for(k=0; k < cctk_lsh[2]; k++)
-      {
-        tx[CCTK_GFINDEX3D(cctkGH,i,j,k)] = i;
-        ty[CCTK_GFINDEX3D(cctkGH,i,j,k)] = j;
-        tz[CCTK_GFINDEX3D(cctkGH,i,j,k)] = k;
-      }
-    }
-  }
-}
-
-
-void TestSlab(CCTK_ARGUMENTS) {  
-  DECLARE_CCTK_ARGUMENTS
-
-  int gs[3], di[3], ds[1]; /* global start/direction/downsampling*/
-  int ls[3], le[3];        /* local start/end */
-  int dsize[1];          /* number of points */
-  int nprocs,iproc;        /* number of procs */
-
-  int l;
-
-
-  int vindex, vdim, vtimelvl;
-  CCTK_REAL *dptr;
-
-  /*$vindex   = CCTK_VarIndex("hyperslab::tz");$*/
-  vindex   = CCTK_VarIndex("grid::z");
-  vtimelvl = CCTK_NumTimeLevelsFromVarI (vindex) - 1;
-  vdim     = CCTK_GroupDimFromVarI(vindex);
-
-  gs[0]=1;
-  gs[1]=1;
-  gs[2]=1;
- 
-  di[0]= 0;
-  di[1]= 0;
-  di[2]= 1;
-
-
-  ds[0]=1;
-
-  printf("Calling CollectData \n");
-  CollectData1D(cctkGH, vindex, vtimelvl, vdim, 
-		gs, di, ds, 10, 
-		&dptr, dsize, 0);
-
-  if (CCTK_MyProc (cctkGH) == 0)
-  {
-    printf("TestSlab: %d \n",dsize[0]);
-    for (l=0;l<dsize[0];l++)
-      printf("%f \n",dptr[l]); 
-  }
-
-}
diff --git a/src/TestSlab2D.c b/src/TestSlab2D.c
deleted file mode 100644
index f8cf6fa..0000000
--- a/src/TestSlab2D.c
+++ /dev/null
@@ -1,67 +0,0 @@
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-#include "Hyperslab.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-#include "CactusBase/IOUtil/src/ioGH.h"
-
-#define MAX_DIM 3
-#define MAX_FACE 6
-#define EPS 1e-10
-#define BAD -42
-
-#define ABS(a) ((a)<0 ? -(a) : (a))
-
-
-void TestSlab2D(CCTK_ARGUMENTS) {  
-  DECLARE_CCTK_ARGUMENTS
-
-  int gs[3], di[2];    /* global start/direction (2times) */
-  int ds[2], len[2];           /* downsampling/length */
-  int ls[3], le[3];            /* local start/end */
-  int dirsize[2],  totsize;      /* number of points (dir./tot)*/
-  int locstart[3], locend[3];
-  int nprocs,iproc;            /* number of procs */
-
-  int l;
-  int *index1, *index2;
-
-  int vindex, vdim, vtimelvl;
-  CCTK_REAL *dptr;
-
-  /*$vindex   = CCTK_VarIndex("hyperslab::tz");$*/
-  vindex   = CCTK_VarIndex("grid::z");
-  vtimelvl = CCTK_NumTimeLevelsFromVarI (vindex) - 1;
-  vdim     = CCTK_GroupDimFromVarI(vindex);
-
-  gs[0]=6;
-  gs[1]=6;
-  gs[2]=3;
- 
-  di[0]= 2;
-  di[1]= 1;
-
-  ds[0]=1;
-  ds[1]=1;
-
-  len[0] = -1;
-  len[1] = -1;
-
-  
-
-  printf("Calling CollectLocalDataND \n");
-  CollectLocalDataND(cctkGH, vindex, vtimelvl, vdim, 2,
-		     gs, di, ds, len, 
-		     &dptr, dirsize, &totsize);
-
-}
-  /*$GetLocalSlab(cctkGH, vdim, 2,
-	       gs, di, ds, len,
-	       locstart, locend, dirsize, &totsize);$*/
diff --git a/src/make.code.defn b/src/make.code.defn
index b448166..2597d24 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,8 +2,4 @@
 # $Header$
 
 # Source files in this directory
-SRCS = GetLocalLine.c GetLocalSlab.c TestSlab.c TestSlab2D.c CollectData1D.c Misc.c CollectDataND.c
-
-# Subdirectories containing source files
-SUBDIRS = 
-
+SRCS = Hyperslab.c CollectData1D.c GetLocalLine.c