Multiprocessor bugfix for 1D diagonal hyperslab extraction from 3D grid arrays.

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

2 files changed, 32 insertions, 24 deletions

diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c
index aaebbf2..adfd7d3 100644
--- a/src/GetHyperslab.c
+++ b/src/GetHyperslab.c
@@ -592,7 +592,7 @@ static int GetDiagonalFromFrom3D (const cGH *GH,
                                   t_hslabConversionFn conversion_fn,
                                   void *hdata)
 {
-  int i, j, k, npoints, myproc, linear_idx;
+  int i, j, k, myproc, linear_idx;
   int vdatatype, hdata_size, vdata_size;
   const char *vdata;
   const pGExtras *extras;
@@ -608,22 +608,28 @@ static int GetDiagonalFromFrom3D (const cGH *GH,
 
   myproc = CCTK_MyProc (GH);
 
-  i = extras->ownership[0][0][0];
-  j = extras->ownership[0][0][1];
-  k = extras->ownership[0][0][2];
-  for (npoints = 0; npoints < mapping->global_hsize[0]; npoints++)
+  i = mapping->global_startpoint[0] - extras->lb[myproc][0];
+  j = mapping->global_startpoint[1] - extras->lb[myproc][1];
+  k = mapping->global_startpoint[2] - extras->lb[myproc][2];
+  while (i < extras->ownership[0][1][0] &&
+         j < extras->ownership[0][1][1] &&
+         k < extras->ownership[0][1][2])
   {
-    linear_idx = i + j*extras->hyper_volume[1] + k*extras->hyper_volume[2];
-    if (conversion_fn)
+    if (i >= extras->ownership[0][0][0] &&
+        j >= extras->ownership[0][0][1] &&
+        k >= extras->ownership[0][0][2])
     {
-      conversion_fn (vdata + linear_idx*vdata_size, hdata, 1, 1, 1);
-    }
-    else
-    {
-      memcpy (hdata, vdata + linear_idx*vdata_size, hdata_size);
+      linear_idx = i + j*extras->hyper_volume[1] + k*extras->hyper_volume[2];
+      if (conversion_fn)
+      {
+        conversion_fn (vdata + linear_idx*vdata_size, hdata, 1, 1, 1);
+      }
+      else
+      {
+        memcpy (hdata, vdata + linear_idx*vdata_size, hdata_size);
+      }
+      hdata = (char *) hdata + hdata_size;
     }
-
-    hdata = (char *) hdata + hdata_size;
     i += mapping->downsample[0];
     j += mapping->downsample[1];
     k += mapping->downsample[2];
diff --git a/src/Mapping.c b/src/Mapping.c
index 4434c8f..2176494 100644
--- a/src/Mapping.c
+++ b/src/Mapping.c
@@ -208,12 +208,12 @@ static CCTK_INT DefineMapping (const cGH      *GH,
   }
   if (! retval)
   {
-    mapping = (hslab_mapping_t *) malloc (sizeof (hslab_mapping_t));
+    mapping = malloc (sizeof (hslab_mapping_t));
     if (mapping)
     {
-      mapping->do_dir = (int *) malloc ((6*vinfo.dim + 2*nprocs) * sizeof(int));
-      mapping->global_hsize = (CCTK_INT *) malloc (((2*hdim + 1) * (nprocs+1) +
-                                                   hdim) * sizeof (CCTK_INT));
+      mapping->do_dir = malloc ((6*vinfo.dim + 2*nprocs) * sizeof(int));
+      mapping->global_hsize = malloc (((2*hdim + 1) * (nprocs+1) + hdim) *
+                                      sizeof (CCTK_INT));
     }
     if (! mapping || ! mapping->do_dir || ! mapping->global_hsize)
     {
@@ -290,7 +290,7 @@ static CCTK_INT DefineMapping (const cGH      *GH,
     }
   }
 
-  /* diagonals can be extracted from non-staggered 3D variables only */ 
+  /* diagonals can be extracted from non-staggered 3D variables only */
   if (mapping->is_diagonal_in_3D && vinfo.stagtype != 0)
   {
     free (mapping->do_dir); free (mapping->global_hsize); free (mapping);
@@ -408,15 +408,17 @@ static CCTK_INT DefineMapping (const cGH      *GH,
       mapping->global_startpoint[vdim] = origin[vdim];
     }
 
+    mapping->totals = 0;
     i = mapping->global_startpoint[0] - extras->lb[myproc][0];
     j = mapping->global_startpoint[1] - extras->lb[myproc][1];
     k = mapping->global_startpoint[2] - extras->lb[myproc][2];
-    mapping->totals = 0;
-    for (npoints = 0; npoints < mapping->global_hsize[0]; npoints++)
+    while (i < extras->ownership[0][1][0] &&
+           j < extras->ownership[0][1][1] &&
+           k < extras->ownership[0][1][2])
     {
-      if (i >= extras->ownership[0][0][0] && i < extras->ownership[0][1][0] &&
-          j >= extras->ownership[0][0][1] && j < extras->ownership[0][1][1] &&
-          k >= extras->ownership[0][0][2] && k < extras->ownership[0][1][2])
+      if (i >= extras->ownership[0][0][0] &&
+          j >= extras->ownership[0][0][1] &&
+          k >= extras->ownership[0][0][2])
       {
         mapping->totals++;
       }