Fix for the reduction of grid arrays which don't have any point on a specific

processor.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHReduce/trunk@42 d60812e6-3970-4df4-986e-c251b06effeb

12 files changed, 185 insertions, 96 deletions

diff --git a/src/Reduction.c b/src/Reduction.c
index a088414..30ddfd5 100644
--- a/src/Reduction.c
+++ b/src/Reduction.c
@@ -176,27 +176,15 @@ int PUGH_ReductionArrays (const cGH *GH,
   buffer = malloc (num_outvals * sizeof (CCTK_REAL));
 
   /* do the reduction on the input arrays */
-  retval = reduction_fn (GH,
-                         proc,
-                         num_dims,
-                         from,
-                         to,
-                         iterator,
-                         points_per_dim,
-                         num_points,
-                         num_inarrays,
-                         intypes,
-                         inarrays,
-                         num_outvals,
-                         buffer);
+  retval = reduction_fn (GH, proc, num_dims, from, to, iterator, points_per_dim,
+                         num_points, 1, num_inarrays, intypes, inarrays,
+                         num_outvals, buffer);
 
   if (retval == 0 && (proc < 0 || proc == CCTK_MyProc (GH)))
   {
     /* type-cast the result to the requested datatype */
     retval = copy_real_to_outtype (num_inarrays * num_outvals,
-                                   buffer,
-                                   outtype,
-                                   outvals);
+                                   buffer, outtype, outvals);
   }
 
   free (intypes);
@@ -305,12 +293,12 @@ int PUGH_ReductionGVs (const cGH *GH,
       case CCTK_GF:
       case CCTK_ARRAY:
         this_retval = ReductionGA (GH, invars[i], proc, &result,
-                                        reduction_fn);
+                                   reduction_fn);
         break;
 
       case CCTK_SCALAR:
         this_retval = ReductionScalar (GH, invars[i], proc, &result,
-                                            reduction_fn);
+                                       reduction_fn);
         break;
 
       default:
@@ -381,7 +369,6 @@ static int ReductionGA (const cGH *GH, int vindex, int proc, CCTK_REAL *outval,
   int i, stagger_index, num_points, dir_points, retval;
   pGA *GA;
   const void *data;
-  char *fullname;
   int *from, *to, *iterator, *points_per_dim;
 
 
@@ -389,17 +376,6 @@ static int ReductionGA (const cGH *GH, int vindex, int proc, CCTK_REAL *outval,
   GA = ((pGA ***) PUGH_pGH (GH)->variables)[vindex][0];
   data = GA->data;
 
-  /* check for zero-sized arrays */
-  if (GA->extras->npoints <= 0)
-  {
-    fullname = CCTK_FullName (vindex);
-    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
-               "ReductionGA: Cannot reduce zero-sized grid array '%s'",
-               fullname);
-    free (fullname);
-    return (-1);
-  }
-
   /* allocate temporary vectors */
   from           = malloc (4 * GA->connectivity->dim * sizeof (int));
   to             = from + 1*GA->connectivity->dim;
@@ -441,8 +417,8 @@ static int ReductionGA (const cGH *GH, int vindex, int proc, CCTK_REAL *outval,
 
   /* now do the reduction */
   retval = reduction_fn (GH, proc, GA->connectivity->dim, from, to, iterator,
-                         points_per_dim, num_points, 1, &GA->vtype, &data,
-                         1, outval);
+                         points_per_dim, num_points, GA->extras->npoints > 0, 1,
+                         &GA->vtype, &data, 1, outval);
 
   /* free temporary vectors */
   free (from);
@@ -515,7 +491,7 @@ static int ReductionScalar (const cGH *GH,
 
   /* now do the reduction */
   retval = reduction_fn (GH, proc, 1, &from, &to, &iterator,
-                         &points_per_dim, num_points, 1, &type, &data,
+                         &points_per_dim, num_points, 1, 1, &type, &data,
                          1, outval);
 
   return retval;
diff --git a/src/ReductionAvg.c b/src/ReductionAvg.c
index 792202b..a2bc35d 100644
--- a/src/ReductionAvg.c
+++ b/src/ReductionAvg.c
@@ -27,6 +27,7 @@ static int ReductionAvg (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -205,6 +206,7 @@ static int ReductionAvg (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -321,8 +323,16 @@ static int ReductionAvg (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local average */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points) 
+    { 
+      /* outvals[] contains now the local average */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionCount.c b/src/ReductionCount.c
index ee4c627..3ea4b7c 100644
--- a/src/ReductionCount.c
+++ b/src/ReductionCount.c
@@ -27,6 +27,7 @@ static int ReductionCount (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -200,6 +201,7 @@ static int ReductionCount (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -217,6 +219,7 @@ static int ReductionCount (const cGH *GH,
   (void) (to + 0);
   (void) (iterator + 0);
   (void) (points_per_dim + 0);
+  (void) (have_local_points + 0);
   (void) (intypes + 0);
   (void) (inarrays + 0);
   (void) (num_outvals + 0);
diff --git a/src/ReductionMax.c b/src/ReductionMax.c
index bc02aa0..a1fceb8 100644
--- a/src/ReductionMax.c
+++ b/src/ReductionMax.c
@@ -11,6 +11,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <float.h>
 
 #include "pugh_reductions.h"
 
@@ -27,6 +28,7 @@ static int ReductionMaxVal (const cGH *GH,
                             int iterator[/* dim */],
                             const int points_per_dim[/* dim */],
                             int num_points,
+                            int have_local_points,
                             int num_inarrays,
                             const int intypes[/* num_inarrays */],
                             const void *const inarrays[/* num_inarrays */],
@@ -205,6 +207,7 @@ static int ReductionMaxVal (const cGH *GH,
                             int iterator[/* dim */],
                             const int points_per_dim[/* dim */],
                             int num_points,
+                            int have_local_points,
                             int num_inarrays,
                             const int intypes[/* num_inarrays */],
                             const void *const inarrays[/* num_inarrays */],
@@ -321,8 +324,19 @@ static int ReductionMaxVal (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local maximum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local maximum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to be smallest possible fp number */
+      for (i = 0; i < total_outvals; i++)
+      {
+        local_outvals[i] = -DBL_MAX;
+      }
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionMin.c b/src/ReductionMin.c
index d28e24f..ab76dfc 100644
--- a/src/ReductionMin.c
+++ b/src/ReductionMin.c
@@ -11,6 +11,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <float.h>
 
 #include "pugh_reductions.h"
 
@@ -27,6 +28,7 @@ static int ReductionMinVal (const cGH *GH,
                             int iterator[/* dim */],
                             const int points_per_dim[/* dim */],
                             int num_points,
+                            int have_local_points,
                             int num_inarrays,
                             const int intypes[/* num_inarrays */],
                             const void *const inarrays[/* num_inarrays */],
@@ -206,6 +208,7 @@ static int ReductionMinVal (const cGH *GH,
                             int iterator[/* dim */],
                             const int points_per_dim[/* dim */],
                             int num_points,
+                            int have_local_points,
                             int num_inarrays,
                             const int intypes[/* num_inarrays */],
                             const void *const inarrays[/* num_inarrays */],
@@ -322,8 +325,19 @@ static int ReductionMinVal (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local minimum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points) 
+    { 
+      /* outvals[] contains now the local minimum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to be largest possible fp number */
+      for (i = 0; i < total_outvals; i++)
+      {
+        local_outvals[i] = DBL_MAX;
+      }
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
index 4c4209b..1394a4a 100644
--- a/src/ReductionNorm1.c
+++ b/src/ReductionNorm1.c
@@ -28,6 +28,7 @@ static int ReductionNorm1 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -209,6 +210,7 @@ static int ReductionNorm1 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -352,8 +354,16 @@ static int ReductionNorm1 (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionNorm2.c b/src/ReductionNorm2.c
index 7d0fc04..35d65f8 100644
--- a/src/ReductionNorm2.c
+++ b/src/ReductionNorm2.c
@@ -29,6 +29,7 @@ static int ReductionNorm2 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -210,6 +211,7 @@ static int ReductionNorm2 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -346,8 +348,16 @@ static int ReductionNorm2 (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionNorm3.c b/src/ReductionNorm3.c
index 2ccb67b..f5036af 100644
--- a/src/ReductionNorm3.c
+++ b/src/ReductionNorm3.c
@@ -29,6 +29,7 @@ static int ReductionNorm3 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -210,6 +211,7 @@ static int ReductionNorm3 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -345,8 +347,16 @@ static int ReductionNorm3 (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionNorm4.c b/src/ReductionNorm4.c
index 1476eb3..8132aa2 100644
--- a/src/ReductionNorm4.c
+++ b/src/ReductionNorm4.c
@@ -29,6 +29,7 @@ static int ReductionNorm4 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -210,6 +211,7 @@ static int ReductionNorm4 (const cGH *GH,
                            int iterator[/* dim */],
                            const int points_per_dim[/* dim */],
                            int num_points,
+                           int have_local_points,
                            int num_inarrays,
                            const int intypes[/* num_inarrays */],
                            const void *const inarrays[/* num_inarrays */],
@@ -343,8 +345,16 @@ static int ReductionNorm4 (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/ReductionNormInf.c b/src/ReductionNormInf.c
index 363ca5b..fa4ea97 100644
--- a/src/ReductionNormInf.c
+++ b/src/ReductionNormInf.c
@@ -12,6 +12,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <float.h>
 
 #include "pugh_reductions.h"
 
@@ -28,6 +29,7 @@ static int ReductionNormInf (const cGH *GH,
                              int iterator[/* dim */],
                              const int points_per_dim[/* dim */],
                              int num_points,
+                             int have_local_points,
                              int num_inarrays,
                              const int intypes[/* num_inarrays */],
                              const void *const inarrays[/* num_inarrays */],
@@ -209,6 +211,7 @@ static int ReductionNormInf (const cGH *GH,
                              int iterator[/* dim */],
                              const int points_per_dim[/* dim */],
                              int num_points,
+                             int have_local_points,
                              int num_inarrays,
                              const int intypes[/* num_inarrays */],
                              const void *const inarrays[/* num_inarrays */],
@@ -361,8 +364,19 @@ static int ReductionNormInf (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to be smallest possible fp number */
+      for (i = 0; i < total_outvals; i++)
+      {
+        local_outvals[i] = -DBL_MAX;
+      }
+    }
     pughGH = PUGH_pGH (GH);
     if (proc < 0)
     {
diff --git a/src/ReductionSum.c b/src/ReductionSum.c
index 46053ac..1610f35 100644
--- a/src/ReductionSum.c
+++ b/src/ReductionSum.c
@@ -27,6 +27,7 @@ static int ReductionSum (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -205,6 +206,7 @@ static int ReductionSum (const cGH *GH,
                          int iterator[/* dim */],
                          const int points_per_dim[/* dim */],
                          int num_points,
+                         int have_local_points,
                          int num_inarrays,
                          const int intypes[/* num_inarrays */],
                          const void *const inarrays[/* num_inarrays */],
@@ -321,8 +323,16 @@ static int ReductionSum (const cGH *GH,
   {
     local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outvals[] contains now the local sum */
-    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (have_local_points)
+    {
+      /* outvals[] contains now the local sum */
+      memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    }
+    else
+    {
+      /* initialize local values to zero */
+      memset (local_outvals, 0, total_outvals * sizeof (CCTK_REAL));
+    }
     if (proc < 0)
     {
       CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
index b65cdfb..f11d012 100644
--- a/src/pugh_reductions.h
+++ b/src/pugh_reductions.h
@@ -67,65 +67,72 @@
       outvals_type typed_outval;                                              \
                                                                               \
                                                                               \
-      for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                     \
+      if (have_local_points)                                                  \
       {                                                                       \
-        /* get the linear index of the element to start with */               \
-        _vindex = from[0];                                                    \
-        for (_i = 1; _i < vdim; _i++)                                         \
-          _vindex += from [_i] * points_per_dim [_i];                         \
-        typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);        \
-                                                                              \
-        /* set iterator to local startpoint */                                \
-        memcpy (iterator, from, vdim * sizeof (int));                         \
-                                                                              \
-        /* do the nested loops starting with the second-innermost */          \
-        _dim = 1;                                                             \
-        while (1)                                                             \
+        for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                   \
         {                                                                     \
-          /* get the linear index */                                          \
-          _vindex = 0;                                                        \
+          /* get the linear index of the element to start with */             \
+          _vindex = from[0];                                                  \
           for (_i = 1; _i < vdim; _i++)                                       \
-            _vindex += iterator [_i] * points_per_dim [_i];                   \
+            _vindex += from [_i] * points_per_dim [_i];                       \
+          typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);      \
                                                                               \
-          /* do the reduction for the innermost loop (lowest dimension) */    \
-          for (_i = from [0]; _i < to [0]; _i++)                              \
-          {                                                                   \
-            REDUCTION_OPERATION (typed_outval, typed_vdata [_i + _vindex]);   \
-          }                                                                   \
+          /* set iterator to local startpoint */                              \
+          memcpy (iterator, from, vdim * sizeof (int));                       \
                                                                               \
-          if (vdim > 1)                                                       \
+          /* do the nested loops starting with the second-innermost */        \
+          _dim = 1;                                                           \
+          while (1)                                                           \
           {                                                                   \
-            /* increment current looper and check for end */                  \
-            if (++iterator [_dim] >= to [_dim])                               \
+            /* get the linear index */                                        \
+            _vindex = 0;                                                      \
+            for (_i = 1; _i < vdim; _i++)                                     \
+              _vindex += iterator [_i] * points_per_dim [_i];                 \
+                                                                              \
+            /* do the reduction for the innermost loop (lowest dimension) */  \
+            for (_i = from [0]; _i < to [0]; _i++)                            \
             {                                                                 \
-              /* increment outermost loopers */                               \
-              for (_dim++; _dim < vdim; _dim++)                               \
+              REDUCTION_OPERATION (typed_outval, typed_vdata [_i + _vindex]); \
+            }                                                                 \
+                                                                              \
+            if (vdim > 1)                                                     \
+            {                                                                 \
+              /* increment current looper and check for end */                \
+              if (++iterator [_dim] >= to [_dim])                             \
               {                                                               \
-                if (++iterator [_dim] < to [_dim])                            \
-                  break;                                                      \
-              }                                                               \
+                /* increment outermost loopers */                             \
+                for (_dim++; _dim < vdim; _dim++)                             \
+                {                                                             \
+                  if (++iterator [_dim] < to [_dim])                          \
+                    break;                                                    \
+                }                                                             \
                                                                               \
-              /* done if beyond outermost loop */                             \
-              if (_dim >= vdim)                                               \
-                break;                                                        \
+                /* done if beyond outermost loop */                           \
+                if (_dim >= vdim)                                             \
+                  break;                                                      \
                                                                               \
-              /* reset innermost loopers */                                   \
-              for (_dim--; _dim >= 0; _dim--)                                 \
-                iterator [_dim] = from [_dim];                                \
-              _dim = 1;                                                       \
+                /* reset innermost loopers */                                 \
+                for (_dim--; _dim >= 0; _dim--)                               \
+                  iterator [_dim] = from [_dim];                              \
+                _dim = 1;                                                     \
+              }                                                               \
+            }                                                                 \
+            else                                                              \
+            {                                                                 \
+              /* exit loop if array is one-dimensional */                     \
+              break;                                                          \
             }                                                                 \
-          }                                                                   \
-          else                                                                \
-          {                                                                   \
-            /* exit loop if array is one-dimensional */                       \
-            break;                                                            \
-          }                                                                   \
                                                                               \
-        } /* end of nested loops over all dimensions */                       \
+          } /* end of nested loops over all dimensions */                     \
                                                                               \
-        outvals [total_outvals++] = (CCTK_REAL) typed_outval;                 \
+          outvals [total_outvals++] = (CCTK_REAL) typed_outval;               \
                                                                               \
-      } /* end of loop over num_outvals */                                    \
+        } /* end of loop over num_outvals */                                  \
+      }                                                                       \
+      else                                                                    \
+      {                                                                       \
+        total_outvals += num_outvals;                                         \
+      }                                                                       \
     }
 
 
@@ -163,6 +170,7 @@ typedef int (*reduction_fn_t) (const cGH *GH,
                                int iterator[/* dim */],
                                const int points_per_dim[/* dim */],
                                int num_points,
+                               int have_local_points,
                                int num_inarrays,
                                const int intypes[/* num_inarrays */],
                                const void *const inarrays[/* num_inarrays */],