All reduction operators are stagger-aware now and can be applied to

arrays of arbitrary dimensions.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGH/trunk@256 b61c5cb5-eaca-4651-9a7a-d64986f99364

7 files changed, 1325 insertions, 1851 deletions

diff --git a/src/Reduction.c b/src/Reduction.c
index c0315b8..89e7cab 100644
--- a/src/Reduction.c
+++ b/src/Reduction.c
@@ -8,387 +8,336 @@
    @version $Header$
  @@*/
 
-#include <stdio.h>
 #include <stdlib.h>
-#include <math.h>
 
-#include "cctk.h"
-#include "pugh.h"
+#include "pugh_reductions.h"
 
 static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_Reduction_c)
 
 /* local function prototypes */
-int PUGH_Sum_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                    void *inArray, int inType, CCTK_REAL outVals []);
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
-int getSizeofElement (int varType);
-
-
-/*@@
-  @routine PUGH_ReduceArraySum
-  @author  Thomas Radke
-  @date    19 Aug 1999
-  @desc
-           Registered PUGH reduction routine for computing the "sum"
-           of a set of arrays.
-           The arrays are described by their dimensions and variable type.
-           For the number of output values only 1 is accepted.
-           Type casting of the result is provided by specifying the
-           requested output datatype. The intermediate reduction value
-           is always computed as a CCTK_REAL value internally.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   int *
-  @vio     in
-  @endvar 
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     inArrays
-  @vdesc   field of input arrays
-  @vtype   void **
-  @vio     in
-  @endvar 
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @vdesc   number of values per output array
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     outVals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
-  @vio     in
-  @endvar 
-  @var     outType
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-@@*/
-int PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
-                          int nArrays, void *inArrays [], int inType,
-                          int nOutVals, void *outVals, int outType)
+static int PUGH_ReductionGA (cGH *GH, int index, int proc, CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn);
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray[/* num_elems */],
+                                 int outtype,
+                                 void *outarray /* [num_elems] */);
+
+
+int PUGH_ReductionArrays (cGH *GH,
+                          int proc,
+                          int num_dims,
+                          int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          void *inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn)
 {
-  int retval=0;
-  int i, elemSize, nPoints;
+  int retval;
+  int i, num_points;
+  int from[1], to[1], iterator[1], points_per_dim[1];
+  int *intypes;
   CCTK_REAL *buffer;
 
 
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
+  points_per_dim[0] = 1;
+  from[0] = 0;
+  to[0] = dims[0];
+  for (i = 1; i < num_dims; i++)
   {
-    nPoints *= dims [i];
+    to[0] *= dims[i];
   }
 
-  if (nOutVals != 1) 
+  if (num_outvals != 1)
   {
-    if (nOutVals != nPoints) 
+    if (num_outvals != to[0])
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArraySum: Don't know how to reduce a "
-		  "%d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return -1;
+                  "PUGH_ReductionArrays: Don't know how to reduce "
+                  "a %d-dimensional array with %d elements "
+                  "to an output array of %d elements",
+                  num_dims, to[0], num_outvals);
+      return (-1);
     }
-    nPoints = 1;
+    to[0] = 1;
   }
+  num_points = to[0] * CCTK_nProcs (GH);
 
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArraySum: Unsupported return type");
-    return -1;
-  }
+  /* set the array types to intype */
+  /* FIXME: could allow to pass in arrays of different types now !!! */
+  intypes = (int *) malloc (num_inarrays * sizeof (int));
+  for (i = 0; i < num_inarrays; i++)
+    intypes[i] = intype;
 
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+  buffer = (CCTK_REAL *) malloc (num_outvals * sizeof (CCTK_REAL));
 
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_Sum_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                        buffer) == 0) 
-    {
+  /* 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);
 
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-      }
-      else
-      {
-	retval = -1;
-      }
-    }
+  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);
   }
 
+  free (intypes);
   free (buffer);
-  return retval;
+
+  return (retval);
 }
 
 
-/*********************************************************************/
-/*                        local functions                            */
-/*********************************************************************/
+int PUGH_ReductionGVs (cGH *GH,
+                       int proc,
+                       int num_invars,
+                       int invars[/* num_invars */],
+                       int outtype,
+                       int num_outvals,
+                       void *outvals /* [num_outvals] */,
+                       reduction_fn_t reduction_fn)
+{
+  int i, myproc, outtypesize, this_retval, retval;
+  CCTK_REAL result;
 
 
-/*@@
-   @routine    PUGH_Sum_Array
-   @date       Sep 14 1999
-   @author     Thomas Radke
-   @desc       Returns the sum of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-   @enddesc 
-@@*/
-int PUGH_Sum_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                    void *inArray, int inType, CCTK_REAL outVals [])
-{
-  int i, j;
-#ifdef CCTK_MPI
-  pGH *pughGH = PUGH_pGH (GH);
-#endif
+  if (num_outvals != 1)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Only one output value per "
+                  "input array allowed");
+    return (-1);
+  }
+
+  outtypesize = CCTK_VarTypeSize (outtype);
+  if (outtypesize <= 0)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Invalid output variable type");
+    return (-1);
+  }
 
-  for (i = 0; i < nOutVals; i++) 
+  retval = 0;
+  myproc = CCTK_MyProc (GH);
+
+  for (i = 0; i < num_invars; i++)
   {
-    outVals [i] = 0.0;
-    for (j = 0; j < nPoints; j++) 
+    switch (CCTK_GroupTypeFromVarI (invars[i]))
     {
-      if (inType == CCTK_VARIABLE_CHAR)
-      {
-        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
-      }
-      else if (inType == CCTK_VARIABLE_INT)
-      {
-        outVals [i] += ((CCTK_INT *) inArray) [i*nPoints + j];
-      }
-#ifdef CCTK_INT2
-      else if (inType == CCTK_VARIABLE_INT2)
-      {
-        outVals [i] += ((CCTK_INT2 *) inArray) [i*nPoints + j];
-      }
-#endif
-#ifdef CCTK_INT4
-      else if (inType == CCTK_VARIABLE_INT4)
-      {
-        outVals [i] += ((CCTK_INT4 *) inArray) [i*nPoints + j];
-      }
-#endif
-#ifdef CCTK_INT8
-      else if (inType == CCTK_VARIABLE_INT8)
-      {
-        outVals [i] += ((CCTK_INT8 *) inArray) [i*nPoints + j];
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_REAL)
-      {
-        outVals [i] += ((CCTK_REAL *) inArray) [i*nPoints + j];
-      }
-#ifdef CCTK_REAL4
-      else if (inType == CCTK_VARIABLE_REAL4)
-      {
-        outVals [i] += ((CCTK_REAL4 *) inArray) [i*nPoints + j];
-      }
-#endif
-#ifdef CCTK_REAL8
-      else if (inType == CCTK_VARIABLE_REAL8)
-      {
-        outVals [i] += ((CCTK_REAL8 *) inArray) [i*nPoints + j];
-      }
-#endif
-#ifdef CCTK_REAL16
-      else if (inType == CCTK_VARIABLE_REAL16)
-      {
-        outVals [i] += ((CCTK_REAL16 *) inArray) [i*nPoints + j];
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_COMPLEX) 
-      {
-        CCTK_WARN (1, "Don't know how to reduce complex variables !!!");
-        return (1);
-      } else 
-      {
-        CCTK_WARN (1, "Unknown variable type in PUGH_Sum_Array");
-        return (1);
-      }
+      case GROUP_GF:
+      case GROUP_ARRAY:
+        this_retval = PUGH_ReductionGA (GH, invars[i], proc, &result,
+                                        reduction_fn);
+        break;
+
+      case GROUP_SCALAR:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Reducing scalars doesn't make sense");
+        this_retval = -1;
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Invalid variable index");
+        this_retval = -1;
+        break;
     }
+
+    if (this_retval == 0 && (proc < 0 || proc == myproc))
+    {
+      /* type-cast the result to the requested datatype */
+      this_retval = copy_real_to_outtype (1, &result, outtype,
+                                          (char *) outvals + i*outtypesize);
+    }
+
+    retval |= this_retval;
   }
 
-#ifdef CCTK_MPI
+  return (retval);
+}
+
+/*****************************************************************************/
+/*                            local functions                                */
+/*****************************************************************************/
+static int PUGH_ReductionGA (cGH *GH, int index, int proc, CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn)
+{
+  int i, retval, stagger_index, num_points;
+  pGA *GA;
+  int *from, *to, *iterator, *points_per_dim;
 
-  if (pughGH->nprocs > 1) 
+
+  /* get the GA structure for the current timelevel */
+  i = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (i > 0)
+  {
+    i--;
+  }
+  GA = ((pGA ***) PUGH_pGH (GH)->variables)[index][i];
+
+  /* allocate temporary vectors */
+  from           = (int *) malloc (4 * GA->connectivity->dim * sizeof (int));
+  to             = from + 1*GA->connectivity->dim;
+  iterator       = from + 2*GA->connectivity->dim;
+  points_per_dim = from + 3*GA->connectivity->dim;
+
+  /* get the start and endpoint to iterate over, the total number of points,
+     and the points_per_dim[] vector */
+  num_points = 1;
+  points_per_dim[0] = 1;
+  for (i = 0; i < GA->connectivity->dim; i++)
   {
-    CCTK_REAL *localOutVals;
+    stagger_index = CCTK_StaggerDirIndex (i, GA->stagger);
 
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+    from[i] = GA->extras->ownership[stagger_index][0][i];
+    to[i]   = GA->extras->ownership[stagger_index][1][i];
 
-    /* outVals[] contains now the local sum */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
-    if (proc < 0)
+    if (stagger_index == 0)
     {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+      num_points *= GA->extras->nsize[i];
     }
     else
     {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+      num_points *= GA->extras->nsize[i] - 1;
     }
 
-    free (localOutVals);
+    if (i > 0)
+    {
+      points_per_dim[i] = points_per_dim[i-1] * GA->extras->lnsize[i-1];
+    }
   }
 
-#endif
+  /* now do the reduction */
+  retval = reduction_fn (GH, proc, GA->connectivity->dim, from, to, iterator,
+                         points_per_dim, num_points, 1, &GA->vtype, &GA->data,
+                         1, outval);
 
-  return (0);
-}
+  /* free temporary vectors */
+  free (from);
 
+  return (retval);
+}
 
 
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n)
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray[/* num_elems */],
+                                 int outtype,
+                                 void *outarray /* [num_elems] */)
 {
-  int j;
+  int i, retval;
+
 
-  for (j = 0; j < n; j++) 
+  retval = 0;
+
+  if (outtype == CCTK_VARIABLE_CHAR)
   {
-    if (to_type == CCTK_VARIABLE_CHAR)
+    CCTK_CHAR *_outarray = (CCTK_CHAR *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
     {
-      ((CCTK_CHAR *) to) [j] = (CCTK_CHAR) (from [j]);
+      _outarray[i] = (CCTK_CHAR) inarray[i];
     }
-    else if (to_type == CCTK_VARIABLE_INT)
+  }
+  else if (outtype == CCTK_VARIABLE_INT)
+  {
+    CCTK_INT *_outarray = (CCTK_INT *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
     {
-      ((CCTK_INT *) to) [j] = (CCTK_INT) (from [j]);
+      _outarray[i] = (CCTK_INT) inarray[i];
     }
+  }
 #ifdef CCTK_INT2
-    else if (to_type == CCTK_VARIABLE_INT2)
+  else if (outtype == CCTK_VARIABLE_INT2)
+  {
+    CCTK_INT2 *_outarray = (CCTK_INT2 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
     {
-      ((CCTK_INT2 *) to) [j] = (CCTK_INT2) (from [j]);
+      _outarray[i] = (CCTK_INT2) inarray[i];
     }
+  }
 #endif
 #ifdef CCTK_INT4
-    else if (to_type == CCTK_VARIABLE_INT4)
+  else if (outtype == CCTK_VARIABLE_INT4)
+  {
+    CCTK_INT4 *_outarray = (CCTK_INT4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
     {
-      ((CCTK_INT4 *) to) [j] = (CCTK_INT4) (from [j]);
+      _outarray[i] = (CCTK_INT4) inarray[i];
     }
+  }
 #endif
 #ifdef CCTK_INT8
-    else if (to_type == CCTK_VARIABLE_INT8)
-    {
-      ((CCTK_INT8 *) to) [j] = (CCTK_INT8) (from [j]);
-    }
-#endif
-    else if (to_type == CCTK_VARIABLE_REAL)
-    {
-      ((CCTK_REAL *) to) [j] = from [j];
-    }
-#ifdef CCTK_REAL4
-    else if (to_type == CCTK_VARIABLE_REAL4)
-    {
-      ((CCTK_REAL4 *) to) [j] = from [j];
-    }
-#endif
-#ifdef CCTK_REAL8
-    else if (to_type == CCTK_VARIABLE_REAL8)
-    {
-      ((CCTK_REAL8 *) to) [j] = from [j];
-    }
-#endif
-#ifdef CCTK_REAL16
-    else if (to_type == CCTK_VARIABLE_REAL16)
-    {
-      ((CCTK_REAL16 *) to) [j] = from [j];
-    }
-#endif
-    else
+  else if (outtype == CCTK_VARIABLE_INT8)
+  {
+    CCTK_INT8 *_outarray = (CCTK_INT8 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
     {
-      CCTK_WARN (1, "copy_real_to_outtype: Unsupported output type");
+      _outarray[i] = (CCTK_INT8) inarray[i];
     }
   }
-}
+#endif
+  else if (outtype == CCTK_VARIABLE_REAL)
+  {
+    CCTK_REAL *_outarray = (CCTK_REAL *) outarray;
 
 
-int getSizeofElement (int varType)
-{
-  if (varType == CCTK_VARIABLE_REAL)
-  {
-    return (sizeof (CCTK_REAL));
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL) inarray[i];
+    }
   }
 #ifdef CCTK_REAL4
-  else if (varType == CCTK_VARIABLE_REAL4)
+  else if (outtype == CCTK_VARIABLE_REAL4)
   {
-    return (sizeof (CCTK_REAL4));
+    CCTK_REAL4 *_outarray = (CCTK_REAL4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL4) inarray[i];
+    }
   }
 #endif
 #ifdef CCTK_REAL8
-  else if (varType == CCTK_VARIABLE_REAL8)
+  else if (outtype == CCTK_VARIABLE_REAL8)
   {
-    return (sizeof (CCTK_REAL8));
+	CCTK_REAL8 *_outarray = (CCTK_REAL8 *) outarray;
+
+
+	for (i = 0; i < num_elems; i++)
+	{
+	  _outarray[i] = (CCTK_REAL8) inarray[i];
+	}
   }
 #endif
 #ifdef CCTK_REAL16
-  else if (varType == CCTK_VARIABLE_REAL16)
-  {
-    return (sizeof (CCTK_REAL16));
-  }
-#endif
-  else if (varType == CCTK_VARIABLE_INT)
-  {
-    return (sizeof (CCTK_INT));
-  }
-#ifdef CCTK_INT2
-  else if (varType == CCTK_VARIABLE_INT2)
+  else if (outtype == CCTK_VARIABLE_REAL16)
   {
-    return (sizeof (CCTK_INT2));
-  }
-#endif
-#ifdef CCTK_INT4
-  else if (varType == CCTK_VARIABLE_INT4)
-  {
-    return (sizeof (CCTK_INT4));
-  }
-#endif
-#ifdef CCTK_INT8
-  else if (varType == CCTK_VARIABLE_INT8)
-  {
-    return (sizeof (CCTK_INT8));
+    CCTK_REAL16 *_outarray = (CCTK_REAL16 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL16) inarray[i];
+    }
   }
 #endif
-  else if (varType == CCTK_VARIABLE_CHAR)
-  {
-    return (sizeof (CCTK_CHAR));
-  }
-  else if (varType == CCTK_VARIABLE_COMPLEX)
+  else
   {
-    return (sizeof (CCTK_COMPLEX));
+    CCTK_WARN (1, "copy_real_to_outtype: Unsupported output type");
+    retval = -1;
   }
 
-  return (-1);
+  return (retval);
 }
diff --git a/src/ReductionMax.c b/src/ReductionMax.c
index 7ce59fa..ba36e31 100644
--- a/src/ReductionMax.c
+++ b/src/ReductionMax.c
@@ -1,43 +1,45 @@
  /*@@
-   @file      Reduction.c
+   @file      ReductionMax.c
    @date      Thu Apr  3 11:54:53 1997
    @author    Thomas Radke, Paul Walker
    @desc 
-   Various MPI reduction operators. 
+              Defines the PUGH reduction operator to get the maximum
+              for a set of arbitrary arrays.
    @enddesc 
    @version $Header$
  @@*/
 
-#include <stdio.h>
 #include <stdlib.h>
-#include <math.h>
 
-#include "cctk.h"
-#include "pugh.h"
+#include "pugh_reductions.h"
 
 static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionMax_c)
 
-/* some useful macros */
-#ifdef MAX
-#undef MAX
-#endif
-#define MAX(x, y)   ((x) > (y) ? (x) : (y))
-
 /* local function prototypes */
-int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals []);
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
-int getSizeofElement (int varType);
+static int PUGH_ReductionMaxVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
 
 /*@@
-  @routine PUGH_ReduceArrayMaxVal
+  @routine PUGH_ReductionMaxValArrays
   @author  Thomas Radke
   @date    19 Aug 1999
   @desc
            Registered PUGH reduction routine for computing the maxima
-           of a set of arrays.
+           for a set of arrays.
            The arrays are described by their dimensions and variable type.
            For the number of output values only 1 is accepted.
            Type casting of the result is provided by specifying the
@@ -57,7 +59,7 @@ int getSizeofElement (int varType);
   @vtype   int
   @vio     in
   @endvar 
-  @var     nDims
+  @var     num_dims
   @vdesc   number of dimensions of input arrays
   @vtype   int
   @vio     in
@@ -67,113 +69,66 @@ int getSizeofElement (int varType);
   @vtype   int *
   @vio     in
   @endvar 
-  @var     nArrays
+  @var     num_inarrays
   @vdesc   number of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     inArrays
+  @var     inarrays
   @vdesc   field of input arrays
   @vtype   void **
   @vio     in
   @endvar 
-  @var     inType
+  @var     intype
   @vdesc   (common) variable type of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of values per output array
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   pointer to buffer holding the output values
   @vtype   void *
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   (common) variable type of output arrays
   @vtype   int
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceArrayMaxVal (cGH *GH, int proc, int nDims, int dims [],
-                             int nArrays, void *inArrays [], int inType,
-                             int nOutVals, void *outVals, int outType)
+int PUGH_ReductionMaxValArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
 {
-  int retval=1;
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArrayMaxVal: Don't know how to reduce "
-		  "a %d-dimensional array "
-                  "with %d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return -1;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayMaxVal: Unsupported return type");
-    return -1;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-				    buffer) == 0) 
-    {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-	retval = 0;
-      }
-    }
-    else 
-    {
-      retval = -1;
-    }
-  }
-
-  free (buffer);
-  
-  return retval;
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionMaxVal));
 }
 
 
-
 /*@@
-  @routine PUGH_ReduceMaxVal
+  @routine PUGH_ReductionMaxValGVs
   @author  Thomas Radke
   @date    23 Apr 1999
   @desc
            PUGH reduction routine to be registered for calculating the maximum.
-           It just goes through the list of variables and calls the appropriate
-           grouptype reduction routine.
+           It just calls the general routine PUGH_ReductionGVs() passing it
+           the callback to compute the maxima for a set of arrays.
            If the number of requested output values equals the size of the
            variable, global reduction is done element-wise.
-           Otherwise nOutVals must be 1, and global reduction is done on the
+           Otherwise num_outvals must be 1, and global reduction is done on the
            results of the local reductions.
   @enddesc
   @history
@@ -189,309 +144,182 @@ int PUGH_ReduceArrayMaxVal (cGH *GH, int proc, int nDims, int dims [],
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of requested output elements
   @vtype   int
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   type of return value
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   array for the result values to be stored
   @vtype   void *
   @vio     in
   @endvar 
-  @var     numinvars
+  @var     num_invars
   @vdesc   number of variables in the list
   @vtype   int
   @vio     in
   @endvar 
-  @var     varlist
+  @var     invars
   @vdesc   list of variables (given as their global indices)
   @vtype   int *
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceMaxVal (cGH *GH, int proc,  int nOutVals,
-                       int outType, void *outVals, int numinvars, int *varlist)
+int PUGH_ReductionMaxValGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
 {
-  int i, elemSize, nPoints;
-  int retval=1;
-  pGH *pughGH;
-  CCTK_REAL *buffer;
-  pGExtras  *extras;
-
-  pughGH = PUGH_pGH (GH);
-
-  /* Get the extras pointer from first variable index */
-  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
-
-  nPoints = extras->npoints;
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceMaxVal: Don't know how to reduce a GF to "
-		  "an output array "
-                  "of %d elements", nOutVals);
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceMaxVal: Unsupported return type");
-    return -2;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < numinvars; i++) 
-  {
-    int result;
-    int timelevel;
-    void *inArray;
-    int inType;
-
-    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
-    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
-    inType = CCTK_VarTypeI (varlist [i]);
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
-    {
-      case GROUP_GF:
-        result = PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArray,
-                                    inType, buffer);
-        break;
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing arrays not yet supported");
-        break;
-      case GROUP_SCALAR:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing a scalar doesn't make sense");
-        break;
-      default:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                            outType, nOutVals);
-      retval = 0;
-    }
-    else if (result == 1)
-    {
-      retval = -1;
-    }
-  }
-
-  free (buffer);
-
-  return retval;
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionMaxVal));
 }
 
 
-/*********************************************************************/
-/*                        local functions                            */
-/*********************************************************************/
-
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
 /*@@
-   @routine    PUGH_MaxVal_Array
+   @routine    PUGH_ReductionMaxVal
    @date       Aug 19 1999
    @author     Thomas Radke
-   @desc 
-               Returns the maximum value of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-   @enddesc 
+   @desc       Returns the maximum of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+   @enddesc
 @@*/
-int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals [])
+static int PUGH_ReductionMaxVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
 {
-  int i, j;
+  int i, j, total_outvals;
 #ifdef CCTK_MPI
-  pGH *pughGH = PUGH_pGH (GH);
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
 #endif
-  
 
-  switch (inType) 
-  {
-    case CCTK_VARIABLE_CHAR:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
 
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_CHAR) (outVals [j]));
-        }
-      }
-      break;
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     ((array)[0])
+#define REDUCTION_OPERATION(max, scalar)   if (max < scalar) max = scalar
 
-    case CCTK_VARIABLE_INT:
-      for (j = 0; j < nOutVals; j++) 
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    for (j = 0; j < num_outvals; j++, total_outvals++)
+    {
+      switch (intypes[i])
       {
-        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
+        case CCTK_VARIABLE_CHAR:
+          ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_CHAR, outvals[total_outvals]);
+          break;
 
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_INT:
+          ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT, outvals[total_outvals]);
+          break;
 
 #ifdef CCTK_INT2
-    case CCTK_VARIABLE_INT2:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT2) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_INT2:
+          ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT2, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_INT4
-    case CCTK_VARIABLE_INT4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT4) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_INT4:
+          ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT4, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_INT8
-    case CCTK_VARIABLE_INT8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT8) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_INT8:
+          ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT8, outvals[total_outvals]);
+          break;
 #endif
 
-    case CCTK_VARIABLE_REAL:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, outVals [j]);
-        }
-      }
-      break;
+        case CCTK_VARIABLE_REAL:
+          ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
 
 #ifdef CCTK_REAL4
-    case CCTK_VARIABLE_REAL4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL4) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_REAL4:
+          ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL4, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_REAL8
-    case CCTK_VARIABLE_REAL8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL8) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_REAL8:
+          ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL8, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_REAL16
-    case CCTK_VARIABLE_REAL16:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL16) (outVals [j]));
-        }
-      }
-      break;
+        case CCTK_VARIABLE_REAL16:
+          ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i], from, to,iterator,
+                         points_per_dim, CCTK_REAL16, outvals[total_outvals]);
+          break;
 #endif
 
-    case CCTK_VARIABLE_COMPLEX:
-      CCTK_WARN(1,"Don't know how to compute maximum of complex variables !!!");
-      return (1);
+        case CCTK_VARIABLE_COMPLEX:
+          CCTK_WARN (1, "PUGH_ReductionMaxVal: Don't know how to compute "
+                        "the maximum of complex variables !!!");
+          return (-1);
 
-    default:
-      CCTK_WARN (1, "Unsupported variable type in PUGH_MaxVal_Array");
-      return (1);
+        default:
+          CCTK_WARN (1, "PUGH_ReductionMaxVal: Unknown variable type");
+          return (-1);
+      }
+    }
   }
 
 #ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
 
-  if (pughGH->nprocs > 1) 
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
   {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outVals[] contains now the local max */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    /* outvals[] contains now the local maximum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
     if (proc < 0)
     {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_MAX, pughGH->PUGH_COMM_WORLD));
     }
     else
     {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_MAX, proc, pughGH->PUGH_COMM_WORLD));
     }
-
-    free (localOutVals);
+    free (local_outvals);
   }
 
 #endif
 
   return (0);
 }
-
-
diff --git a/src/ReductionMin.c b/src/ReductionMin.c
index c16c37a..12bf319 100644
--- a/src/ReductionMin.c
+++ b/src/ReductionMin.c
@@ -3,42 +3,43 @@
    @date      Thu Apr  3 11:54:53 1997
    @author    Thomas Radke, Paul Walker
    @desc 
-   MPI reduction operator to calculate minimum value
+              Defines the PUGH reduction operator to get the minimum
+              for a set of arbitrary arrays.
    @enddesc 
    @version $Header$
  @@*/
 
-#include <stdio.h>
 #include <stdlib.h>
-#include <math.h>
 
-#include "cctk.h"
-#include "pugh.h"
+#include "pugh_reductions.h"
 
 static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionMin_c)
 
-/* some useful macros */
-#ifdef MIN
-#undef MIN
-#endif
-#define MIN(x, y)   ((x) < (y) ? (x) : (y))
-
 /* local function prototypes */
-int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals []);
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
-int getSizeofElement (int varType);
+static int PUGH_ReductionMinVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
 
 
 /*@@
-  @routine PUGH_ReduceArrayMinVal
+  @routine PUGH_ReductionMinValArrays
   @author  Thomas Radke
   @date    19 Aug 1999
   @desc
            Registered PUGH reduction routine for computing the minima
-           of a set of arrays.
+           for a set of arrays.
            The arrays are described by their dimensions and variable type.
            For the number of output values only 1 is accepted.
            Type casting of the result is provided by specifying the
@@ -58,7 +59,7 @@ int getSizeofElement (int varType);
   @vtype   int
   @vio     in
   @endvar 
-  @var     nDims
+  @var     num_dims
   @vdesc   number of dimensions of input arrays
   @vtype   int
   @vio     in
@@ -68,105 +69,57 @@ int getSizeofElement (int varType);
   @vtype   int *
   @vio     in
   @endvar 
-  @var     nArrays
+  @var     num_inarrays
   @vdesc   number of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     inArrays
+  @var     inarrays
   @vdesc   field of input arrays
   @vtype   void **
   @vio     in
   @endvar 
-  @var     inType
+  @var     intype
   @vdesc   (common) variable type of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of values per output array
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   pointer to buffer holding the output values
   @vtype   void *
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   (common) variable type of output arrays
   @vtype   int
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceArrayMinVal (cGH *GH, int proc, int nDims, int dims [],
-                             int nArrays, void *inArrays [], int inType,
-                             int nOutVals, void *outVals, int outType)
+int PUGH_ReductionMinValArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
 {
-  int retval =1;
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArrayMinVal: Don't know how to reduce a "
-		  "%d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return -1;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-      CCTK_WARN (1, "PUGH_ReduceArrayMinVal: Unsupported return type");
-      return -1;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArrays [i],
-                           inType, buffer) == 0) 
-    {
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-	retval = 0;
-      }
-    }
-    else
-    {
-      retval = -1;
-    }
-  }
-  
-  free (buffer);
-
-  return retval;
-
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionMinVal));
 }
 
 
-
-
 /*@@
-  @routine PUGH_ReduceMinVal
+  @routine PUGH_ReductionMinValGVs
   @author  Thomas Radke
   @date    23 Apr 1999
   @desc
@@ -175,7 +128,7 @@ int PUGH_ReduceArrayMinVal (cGH *GH, int proc, int nDims, int dims [],
            grouptype reduction routine.
            If the number of requested output values equals the size of the
            variable, global reduction is done element-wise.
-           Otherwise nOutVals must be 1, and global reduction is done on the
+           Otherwise num_outvals must be 1, and global reduction is done on the
            results of the local reductions.
   @enddesc
   @history
@@ -191,302 +144,182 @@ int PUGH_ReduceArrayMinVal (cGH *GH, int proc, int nDims, int dims [],
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of requested output elements
   @vtype   int
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   type of return value
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   array for the result values to be stored
   @vtype   void *
   @vio     in
   @endvar 
-  @var     numinvars
+  @var     num_invars
   @vdesc   number of variables in the list
   @vtype   int
   @vio     in
   @endvar 
-  @var     varlist
+  @var     invars
   @vdesc   list of variables (given as their global indices)
   @vtype   int *
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceMinVal (cGH *GH, 
-		       int proc, 
-		       int nOutVals,
-		       int outType, 
-		       void *outVals, 
-		       int numinvars, 
-		       int *varlist)
+int PUGH_ReductionMinValGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
 {
-
-  int        i;
-  int        retval=1;
-  int        elemSize;
-  int        nPoints;
-  pGH       *pughGH;
-  CCTK_REAL *buffer;
-  pGExtras  *extras;
-
-  pughGH = PUGH_pGH (GH);
-
-  /* Get the extras pointer from first variable index */
-  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
-
-  nPoints = extras->npoints;
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceMinVal: Don't know how to reduce a GF "
-		  "to an output array "
-                  "of %d elements", nOutVals);
-      return -2;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-    CCTK_WARN (1, "PUGH_ReduceMinVal: Unsupported output datatype");
-    return -3;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < numinvars; i++)
-  {
-    int result;
-    int timelevel;
-    void *inArray;
-    int inType;
-
-    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
-    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
-    inType = CCTK_VarTypeI (varlist [i]);
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
-    {
-      case GROUP_GF:
-        result = PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArray,
-                                    inType, buffer);
-        break;
-
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceMinVal: Reducing arrays not yet supported");
-        break;
-
-      case GROUP_SCALAR:
-        CCTK_WARN(1, "PUGH_ReduceMinVal: Reducing a scalar doesn't make sense");
-        break;
-
-      default:
-        CCTK_WARN (1, "PUGH_ReduceMinVal: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                            outType, nOutVals);
-      retval = 0;
-    }
-    else if (result == 1)
-    {
-      retval = -1;
-    }
-  }
-
-  free (buffer);
-
-  return retval;
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionMinVal));
 }
 
 
-/*********************************************************************/
-/*                        local functions                            */
-/*********************************************************************/
-
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
 /*@@
-   @routine    PUGH_MinVal_Array
+   @routine    PUGH_ReductionMinVal
    @date       Aug 19 1999
    @author     Thomas Radke
-   @desc 
-               Returns the minimum value of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-   @enddesc 
+   @desc       Returns the minimum for a set of distributed array with
+               'num_points' elements each.
+               Global reduction is done element-wise (num_outvals == 1)
+               or on the results of the local reductions.
+   @enddesc
 @@*/
-int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals [])
+static int PUGH_ReductionMinVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
 {
-  int i, j;
+  int i, j, total_outvals;
 #ifdef CCTK_MPI
-  pGH *pughGH = PUGH_pGH (GH);
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
 #endif
-  
-  switch (inType) 
-  {
-    case CCTK_VARIABLE_CHAR:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
 
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr)
-          outVals [j] = MIN (*inPtr, (CCTK_CHAR) (outVals [j]));
-      }
-      break;
 
-    case CCTK_VARIABLE_INT:
-      for (j = 0; j < nOutVals; j++) 
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     ((array)[0])
+#define REDUCTION_OPERATION(min, scalar)   if (min > scalar) min = scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    for (j = 0; j < num_outvals; j++, total_outvals++)
+    {
+      switch (intypes[i])
       {
-        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
+        case CCTK_VARIABLE_CHAR:
+          ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_CHAR, outvals[total_outvals]);
+          break;
 
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT) (outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_INT:
+          ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT, outvals[total_outvals]);
+          break;
 
 #ifdef CCTK_INT2
-    case CCTK_VARIABLE_INT2:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT2) (outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_INT2:
+          ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT2, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_INT4
-    case CCTK_VARIABLE_INT4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT4) (outVals [j]));
-      }
-      break;
-#endif 
+        case CCTK_VARIABLE_INT4:
+          ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT4, outvals[total_outvals]);
+          break;
+#endif
 
 #ifdef CCTK_INT8
-    case CCTK_VARIABLE_INT8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT8) (outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_INT8:
+          ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT8, outvals[total_outvals]);
+          break;
 #endif
 
-    case CCTK_VARIABLE_REAL:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, outVals [j]);
-      }
-      break;
+        case CCTK_VARIABLE_REAL:
+          ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
 
 #ifdef CCTK_REAL4
-    case CCTK_VARIABLE_REAL4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL4)( outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_REAL4:
+          ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL4, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_REAL8
-    case CCTK_VARIABLE_REAL8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL8)( outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_REAL8:
+          ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL8, outvals[total_outvals]);
+          break;
 #endif
 
 #ifdef CCTK_REAL16
-    case CCTK_VARIABLE_REAL16:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL16)( outVals [j]));
-      }
-      break;
+        case CCTK_VARIABLE_REAL16:
+          ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i], from, to,iterator,
+                         points_per_dim, CCTK_REAL16, outvals[total_outvals]);
+          break;
 #endif
 
-    case CCTK_VARIABLE_COMPLEX:
-      CCTK_WARN(1,"Don't know how to compute minimum of complex variables !!!");
-      return (1);
+        case CCTK_VARIABLE_COMPLEX:
+          CCTK_WARN (1, "PUGH_ReductionMinVal: Don't know how to compute "
+                        "the minimum of complex variables !!!");
+          return (-1);
 
-    default:
-      CCTK_WARN (1, "Unsupported variable type in PUGH_MinVal_Array");
-      return (1);
+        default:
+          CCTK_WARN (1, "PUGH_ReductionMinVal: Unknown variable type");
+          return (-1);
+      }
+    }
   }
 
-
 #ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
 
-  if (pughGH->nprocs > 1) 
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
   {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outVals[] contains now the local min */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    /* outvals[] contains now the local minimum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
     if (proc < 0)
     {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_MIN, pughGH->PUGH_COMM_WORLD));
     }
     else
     {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_MIN, proc, pughGH->PUGH_COMM_WORLD));
     }
-    free (localOutVals);
+    free (local_outvals);
   }
-
 #endif
 
   return (0);
 }
-
-
diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
index 94a9d3f..32de553 100644
--- a/src/ReductionNorm1.c
+++ b/src/ReductionNorm1.c
@@ -1,49 +1,41 @@
  /*@@
-   @file      Reduction_nm1.c
+   @file      ReductionNorm1.c
    @date      Thu Apr  3 11:54:53 1997
    @author    Thomas Radke, Paul Walker
    @desc 
-   Various MPI reduction operators. 
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm1 = $\Sigma |a_i| / np$
    @enddesc 
    @version $Header$
  @@*/
 
-#include <stdio.h>
 #include <stdlib.h>
-#include <math.h>
 
-#include "cctk.h"
-#include "pugh.h"
+#include "pugh_reductions.h"
 
 static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm1_c)
 
-/* some useful macros */
-#ifdef MIN
-#undef MIN
-#endif
-#ifdef ABS
-#undef ABS
-#endif
-#ifdef SQR
-#undef SQR
-#endif
-#define MIN(x, y)   ((x) < (y) ? (x) : (y))
-#define ABS(x)      ((x) < 0   ?-(x) : (x))
-#define SQR(x)      ((x) * (x))
-
-
 /* local function prototypes */
-int PUGH_Norm1_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
-int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals []);
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
-int getSizeofElement (int varType);
+static int PUGH_ReductionNorm1 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
 
 
 /*@@
-  @routine PUGH_ReduceArrayNorm1
+  @routine PUGH_ReductionNorm1Arrays
   @author  Thomas Radke
   @date    19 Aug 1999
   @desc
@@ -52,7 +44,7 @@ int getSizeofElement (int varType);
            The arrays are described by their dimensions and variable type.
            If the number of requested output values equals the number of
            elements in an input array, global reduction is done element-wise.
-           Otherwise nOutVals must be 1, and global reduction is done on the
+           Otherwise num_outvals must be 1, and global reduction is done on the
            results of the local reductions.
            Type casting of the result is provided by specifying the
            requested output datatype. The intermediate reduction value
@@ -71,7 +63,7 @@ int getSizeofElement (int varType);
   @vtype   int
   @vio     in
   @endvar 
-  @var     nDims
+  @var     num_dims
   @vdesc   number of dimensions of input arrays
   @vtype   int
   @vio     in
@@ -81,105 +73,57 @@ int getSizeofElement (int varType);
   @vtype   int *
   @vio     in
   @endvar 
-  @var     nArrays
+  @var     num_arrays
   @vdesc   number of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     inArrays
+  @var     inarrays
   @vdesc   field of input arrays
   @vtype   void **
   @vio     in
   @endvar 
-  @var     inType
+  @var     intype
   @vdesc   (common) variable type of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of values per output array
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   pointer to buffer holding the output values
   @vtype   void *
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   (common) variable type of output arrays
   @vtype   int
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
-                            int nArrays, void *inArrays [], int inType,
-                            int nOutVals, void *outVals, int outType)
+int PUGH_ReductionNorm1Arrays (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int dims[/* num_dims */],
+                               int num_inarrays,
+                               void *inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
 {
-  int retval = 1;
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArrayNorm1: Don't know how to reduce "
-                  "a %d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return -1;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type");
-    return -1;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-
-    /* do the reduction on input array [i] */
-    if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                           buffer) == 0) 
-    {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-	retval = 0;
-      }
-    }
-    else
-    {
-      retval = -1;
-    }
-  }
-
-  free (buffer);
-  return retval;
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionNorm1));
 }
 
 
-
-
 /*@@
-  @routine PUGH_ReduceNorm1
+  @routine PUGH_ReductionNorm1GVs
   @author  Thomas Radke
   @date    23 Apr 1999
   @desc
@@ -187,7 +131,7 @@ int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
            It just goes through the list of variables and calls the appropriate
            grouptype reduction routine.
            Global reduction is always done on the results of the local
-           reductions, so nOutVals must be 1.
+           reductions, so num_outvals must be 1.
   @enddesc
   @history
   @endhistory 
@@ -207,444 +151,197 @@ int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
   @vtype   int
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   type of return value
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   array for the result values to be stored
   @vtype   void *
   @vio     in
   @endvar 
-  @var     numinvars
+  @var     num_invars
   @vdesc   number of variables in the list
   @vtype   int
   @vio     in
   @endvar 
-  @var     varlist
+  @var     invars
   @vdesc   list of variables (given as their global indices)
   @vtype   int *
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceNorm1 (cGH *GH, int proc,  
-		       int nOutVals,
-                       int outType, 
-		       void *outVals, 
-		       int numinvars, 
-		       int *varlist)
+int PUGH_ReductionNorm1GVs (cGH *GH,
+                            int proc,  
+                            int num_outvals,
+                            int outtype, 
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            int invars[/* num_invars */])
 {
-  int retval=1;
-  int i, elemSize;
-  pGH *pughGH;
-  CCTK_REAL intermediate_result;
-
-  if (nOutVals != 1) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceNorm1: Only one output value per "
-                  "input array allowed");
-    return -2;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-      CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type");
-      return -3;
-  }
-
-  pughGH = PUGH_pGH (GH);
-
-  for (i = 0; i < numinvars; i++) 
-  {
-    int result;
-
-    result = 1;
-
-    if (CCTK_GroupDimFromVarI(varlist[i])<3)
-    {
-      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
-		 "PUGH_ReduceNorm1: Reducing dimension %d not yet supported",
-		 CCTK_GroupDimFromVarI(varlist[i]));
-    }
-    else
-    {
-		 
-      switch (CCTK_GroupTypeFromVarI (varlist [i]))
-      {
-        case GROUP_GF:
-	  result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result);
-	  break;
-        case GROUP_ARRAY:
-	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported");
-	  break;
-        case GROUP_SCALAR:
-	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense");
-	  break;
-        default:
-	  CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type");
-	  break;
-      }
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (&intermediate_result, (char *) outVals +
-                            i*nOutVals*elemSize, outType, nOutVals);
-      retval = 0;
-    }
-    else if (result == 0)
-    {
-      retval = -1;
-    }
-  }
-
-  return retval;
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionNorm1));
 }
 
 
-
-/*********************************************************************/
-/*                        local functions                            */
-/*********************************************************************/
-
-
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
 /*@@
-   @routine    PUGH_Norm1_Array
+   @routine    PUGH_ReductionNorm1
    @date       Aug 19 1999
    @author     Thomas Radke
    @desc       Returns the "norm1" of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
 
                "norm1" is defined as $\Sigma |a_i| / np$.
    @enddesc 
 @@*/
-int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals [])
+static int PUGH_ReductionNorm1 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
 {
-  int i, j;
-  pGH *pughGH = PUGH_pGH (GH);
+  int i, j, total_outvals;
+  pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(data)       0
+#ifdef ABS
+#undef ABS
+#endif
+#define ABS(x)      ((x) < 0 ? -(x) : (x))
 
 
-  for (i = 0; i < nOutVals; i++) 
+  for (i = total_outvals = 0; i < num_inarrays; i++)
   {
-    outVals [i] = 0.0;
-    for (j = 0; j < nPoints; j++) 
+    for (j = 0; j < num_outvals; j++, total_outvals++)
     {
-      if (inType == CCTK_VARIABLE_CHAR)           /* CCTK_CHAR is unsigned */
+      switch (intypes[i])
       {
-        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
-      }
-      else if (inType == CCTK_VARIABLE_INT)
-      {
-        outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += scalar
+          ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_CHAR, outvals[total_outvals]);
+          break;
+
+        case CCTK_VARIABLE_INT:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += ABS (scalar)
+          ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT, outvals[total_outvals]);
+          break;
+
 #ifdef CCTK_INT2
-      else if (inType == CCTK_VARIABLE_INT2)
-      {
-        outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT2:
+          ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT2, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_INT4
-      else if (inType == CCTK_VARIABLE_INT4)
-      {
-        outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT4:
+          ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT4, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_INT8
-      else if (inType == CCTK_VARIABLE_INT8)
-      {
-        outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT8:
+          ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT8, outvals[total_outvals]);
+          break;
 #endif
-      else if (inType == CCTK_VARIABLE_REAL)
-      {
-        outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]);
-      }
+
+        case CCTK_VARIABLE_REAL:
+          ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
+
 #ifdef CCTK_REAL4
-      else if (inType == CCTK_VARIABLE_REAL4)
-      {
-        outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL4:
+          ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL4, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_REAL8
-      else if (inType == CCTK_VARIABLE_REAL8)
-      {
-        outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL8:
+          ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL8, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_REAL16
-      else if (inType == CCTK_VARIABLE_REAL16)
-      {
-        outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL16:
+          ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL16, outvals[total_outvals]);
+          break;
 #endif
-      else if (inType == CCTK_VARIABLE_COMPLEX) 
-      {
-        outVals [i] += 
-          sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
-                 SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) );
-      } 
-      else 
-      {
-        CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type");
-        return (1);
+
+        case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += CCTK_CmplxAbs (scalar)
+          ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i], from, to,iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
+
+        default:
+          CCTK_WARN (1, "PUGH_ReductionNorm1: Unknown variable type");
+          return (-1);
       }
     }
   }
 
-#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
 
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
   if (pughGH->nprocs > 1) 
   {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outVals[] contains now the local sum */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
     if (proc < 0)
     {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
     }
     else
     {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
     }
-    free (localOutVals);
-
-    nPoints *= pughGH->nprocs;
+    free (local_outvals);
   }
 
 #endif
 
+  /* finally assign the return value */
   if (proc < 0 || proc == pughGH->myproc)
   {
-    for (i = 0; i < nOutVals; i++)
-    {
-      outVals [i] /= nPoints;
-    }
-  }
-
-  return (0);
-}
-
-
-
-/*@@
-   @routine    PUGH_Norm1_GF
-   @date       Thu Apr  3 11:56:05 1997
-   @author     Paul Walker
-   @desc 
-   Returns the average of a distributed grid function.
-   @enddesc 
-@@*/
-int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
-{
-  int timelevel;
-  int ii;
-  CCTK_INT tnp;
-  pGH *pughGH;
-  pGA *GA;
-#ifdef CCTK_MPI
-  int jj, kk;
-  CCTK_INT np;
-  CCTK_REAL local_result;
-#endif
-
-  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
-
-  pughGH = PUGH_pGH (GH);
-  GA = ((pGA ***) pughGH->variables) [index][timelevel];
-
-  *outVal = 0.0;
-
-#ifndef CCTK_MPI
-  for (ii = 0; ii < GA->extras->npoints; ii++)
-  {
-    if (GA->vtype == CCTK_VARIABLE_CHAR)
-    {
-      *outVal += ((CCTK_CHAR *) GA->data) [ii];     /* CCTK_CHAR is unsigned */
-    }
-    else if (GA->vtype == CCTK_VARIABLE_INT)
-    {
-      *outVal += ABS (((CCTK_INT *) GA->data) [ii]);
-    }
-#ifdef CCTK_INT2
-    else if (GA->vtype == CCTK_VARIABLE_INT2)
-    {
-      *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT4
-    else if (GA->vtype == CCTK_VARIABLE_INT4)
-    {
-      *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT8
-    else if (GA->vtype == CCTK_VARIABLE_INT8)
-    {
-      *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_REAL)
-    {
-      *outVal += ABS (((CCTK_REAL *) GA->data) [ii]);
-    }
-#ifdef CCTK_REAL4
-    else if (GA->vtype == CCTK_VARIABLE_REAL4)
-    {
-      *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL8
-    else if (GA->vtype == CCTK_VARIABLE_REAL8)
-    {
-      *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL16
-    else if (GA->vtype == CCTK_VARIABLE_REAL16)
-    {
-      *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-    {
-      *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
-                       SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) );
-    } 
-    else 
-    {
-      CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type");
-      return (1);
-    }
-  }
-
-  tnp = GA->extras->npoints;
-
-#else
-
-  /* MPI version depends heavily on the ownership array */
-  np  = 0;
-  local_result = 0.0;
-
-  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
-       kk < GA->extras->ownership[GA->stagger][1][2]; 
-       kk ++)
-  {
-    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
-         jj < GA->extras->ownership[GA->stagger][1][1]; 
-         jj ++)
+    for (i = 0; i < total_outvals; i++)
     {
-      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
-           ii < GA->extras->ownership[GA->stagger][1][0]; 
-           ii++, np++)
-      {
-        if (GA->vtype == CCTK_VARIABLE_CHAR)
-        {
-          local_result +=              /* CCTK_CHAR is unsigned */
-	    ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)];
-        }
-        else if (GA->vtype == CCTK_VARIABLE_INT)
-        {
-          local_result +=
-	    ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_INT2
-        else if (GA->vtype == CCTK_VARIABLE_INT2)
-        {
-          local_result +=
-	    ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT4
-        else if (GA->vtype == CCTK_VARIABLE_INT4)
-        {
-          local_result +=
-	    ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT8
-        else if (GA->vtype == CCTK_VARIABLE_INT8)
-        {
-          local_result +=
-	    ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_REAL)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_REAL4
-        else if (GA->vtype == CCTK_VARIABLE_REAL4)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL8
-        else if (GA->vtype == CCTK_VARIABLE_REAL8)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL16
-        else if (GA->vtype == CCTK_VARIABLE_REAL16)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-        {
-          local_result +=
-            sqrt( SQR((((CCTK_COMPLEX *)GA->data)
-		       [DATINDEX (GA->extras,ii,jj,kk)]).Re)+
-		  SQR((((CCTK_COMPLEX *)GA->data)
-		       [DATINDEX (GA->extras,ii,jj,kk)]).Im) );  
-        } 
-        else 
-        {
-          CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type");
-          return (1);
-        }
-      }
+      outvals[i] /= num_points;
     }
   }
 
-  if (proc < 0) 
-  {
-    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      pughGH->PUGH_COMM_WORLD));
-  } 
-  else 
-  {
-    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      proc, pughGH->PUGH_COMM_WORLD));
-  }
-#endif
- 
-  if (proc < 0 || proc == pughGH->myproc)
-  {
-    *outVal /= tnp;
-  }
-
   return (0);
 }
-
-
diff --git a/src/ReductionNorm2.c b/src/ReductionNorm2.c
index 797f8f5..0fa6d1d 100644
--- a/src/ReductionNorm2.c
+++ b/src/ReductionNorm2.c
@@ -1,41 +1,42 @@
- /*@@
+/*@@
    @file      ReductionNorm2.c
    @date      Thu Apr  3 11:54:53 1997
    @author    Thomas Radke, Paul Walker
-   @desc 
-   Various MPI reduction operators. 
-   @enddesc 
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm2 = $\sqrt{\Sigma (a_i * a_i) / np}$
+   @enddesc
    @version $Header$
  @@*/
 
-#include <stdio.h>
-#include <stdlib.h>
 #include <math.h>
+#include <stdlib.h>
 
-#include "cctk.h"
-#include "pugh.h"
+#include "pugh_reductions.h"
 
 static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm2_c)
 
-/* some useful macros */
-#ifdef SQR
-#undef SQR
-#endif
-#define SQR(x)      ((x) * (x))
-
-
 /* local function prototypes */
-int PUGH_Norm2_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
-int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals []);
-void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
-int getSizeofElement (int varType);
+static int PUGH_ReductionNorm2 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
 
 
 /*@@
-  @routine PUGH_ReduceArrayNorm2
+  @routine PUGH_ReductionNorm2Arrays
   @author  Thomas Radke
   @date    19 Aug 1999
   @desc
@@ -44,7 +45,7 @@ int getSizeofElement (int varType);
            The arrays are described by their dimensions and variable type.
            If the number of requested output values equals the number of
            elements in an input array, global reduction is done element-wise.
-           Otherwise nOutVals must be 1, and global reduction is done on the
+           Otherwise num_outvals must be 1, and global reduction is done on the
            results of the local reductions.
            Type casting of the result is provided by specifying the
            requested output datatype. The intermediate reduction value
@@ -63,7 +64,7 @@ int getSizeofElement (int varType);
   @vtype   int
   @vio     in
   @endvar 
-  @var     nDims
+  @var     num_dims
   @vdesc   number of dimensions of input arrays
   @vtype   int
   @vio     in
@@ -73,103 +74,57 @@ int getSizeofElement (int varType);
   @vtype   int *
   @vio     in
   @endvar 
-  @var     nArrays
+  @var     num_arrays
   @vdesc   number of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     inArrays
+  @var     inarrays
   @vdesc   field of input arrays
   @vtype   void **
   @vio     in
   @endvar 
-  @var     inType
+  @var     intype
   @vdesc   (common) variable type of input arrays
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     num_outvals
   @vdesc   number of values per output array
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   pointer to buffer holding the output values
   @vtype   void *
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   (common) variable type of output arrays
   @vtype   int
   @vio     in
   @endvar 
 @@*/
-int PUGH_ReduceArrayNorm2 (cGH *GH, int proc, int nDims, int dims [],
-                            int nArrays, void *inArrays [], int inType,
-                            int nOutVals, void *outVals, int outType)
+int PUGH_ReductionNorm2Arrays (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int dims[/* num_dims */],
+                               int num_inarrays,
+                               void *inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
 {
-  int retval = 1;
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArrayNorm2: Don't know how to reduce a "
-		  "%d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return -1;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayNorm2: Unsupported return type");
-    return -1;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_Norm2_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                          buffer) == 0) 
-    {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-	retval = 0;
-      }
-    }
-    else
-    {
-      retval = -1;
-    }
-  }
-
-  free (buffer);
-
-  return retval;
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionNorm2));
 }
 
 
 /*@@
-  @routine PUGH_ReduceNorm2
+  @routine PUGH_ReductionNorm2GVs
   @author  Thomas Radke
   @date    23 Apr 1999
   @desc
@@ -177,7 +132,7 @@ int PUGH_ReduceArrayNorm2 (cGH *GH, int proc, int nDims, int dims [],
            It just goes through the list of variables and calls the appropriate
            grouptype reduction routine.
            Global reduction is always done on the results of the local
-           reductions, so nOutVals must be 1.
+           reductions, so num_outvals must be 1.
   @enddesc
   @history
   @endhistory 
@@ -192,445 +147,201 @@ int PUGH_ReduceArrayNorm2 (cGH *GH, int proc, int nDims, int dims [],
   @vtype   int
   @vio     in
   @endvar 
-  @var     nOutVals
+  @var     nOutVars
   @vdesc   number of requested output elements
   @vtype   int
   @vio     in
   @endvar 
-  @var     outType
+  @var     outtype
   @vdesc   type of return value
   @vtype   int
   @vio     in
   @endvar 
-  @var     outVals
+  @var     outvals
   @vdesc   array for the result values to be stored
   @vtype   void *
   @vio     in
   @endvar 
-  @var     numinvars
+  @var     num_invars
   @vdesc   number of variables in the list
   @vtype   int
   @vio     in
   @endvar 
-  @var     varlist
+  @var     invars
   @vdesc   list of variables (given as their global indices)
   @vtype   int *
   @vio     in
   @endvar 
 @@*/
-
-int PUGH_ReduceNorm2 (cGH *GH, int proc, int nOutVals,
-                       int outType, void *outVals, int numinvars, int *varlist)
+int PUGH_ReductionNorm2GVs (cGH *GH,
+                            int proc,  
+                            int num_outvals,
+                            int outtype, 
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            int invars[/* num_invars */])
 {
-  int retval = 1;
-  int i, elemSize;
-  pGH *pughGH;
-  CCTK_REAL intermediate_result;
-
-
-  if (nOutVals != 1)
-  {
-    CCTK_WARN (1, "PUGH_ReduceNorm2 accepts only one output value per "
-                  "input array");
-    return -2;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-      CCTK_WARN (1, "PUGH_ReduceNorm2: Unsupported return type");
-      return -3;
-  }
-
-  pughGH = PUGH_pGH (GH);
-
-  for (i = 0; i < numinvars; i++)
-  {
-    int result;
-
-    result = 1;
-
-    if (CCTK_GroupDimFromVarI(varlist[i])<3)
-    {
-      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
-		 "PUGH_ReduceNorm2: Reducing dimension %d not yet supported",
-		 CCTK_GroupDimFromVarI(varlist[i]));
-    }
-    else
-    {
-      switch (CCTK_GroupTypeFromVarI (varlist [i])) {
-        case GROUP_GF:
-	  result = PUGH_Norm2_GF (GH, varlist [i], proc, &intermediate_result);
-	  break;
-        case GROUP_ARRAY:
-	  CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing arrays not yet supported");
-	  break;
-        case GROUP_SCALAR:
-	  CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing a scalar doesn't make sense");
-	  break;
-        default:
-	  CCTK_WARN (1, "PUGH_ReduceNorm2: Invalid group type");
-	  break;
-      }
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (&intermediate_result, (char *) outVals +
-                            i*nOutVals*elemSize, outType, nOutVals);
-      retval = 0;
-    }
-    else if (result == 1)
-    {
-      retval = -1;
-    }
-  }
-
-  return retval;
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionNorm2));
 }
 
 
-
-/*********************************************************************/
-/*                        local functions                            */
-/*********************************************************************/
-
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
 /*@@
-   @routine    PUGH_Norm2_Array
+   @routine    PUGH_ReductionNorm2
    @date       Aug 19 1999
    @author     Thomas Radke
-   @desc 
-               Returns the "norm2" of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
+   @desc       Returns the "norm2" of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
 
                "norm2" is defined as $\sqrt{\Sigma (a_i * a_i) / np}$.
    @enddesc 
 @@*/
-int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals [])
+static int PUGH_ReductionNorm2 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
 {
-  int i, j;
+  int i, j, total_outvals;
   pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
 
 
-  pughGH = PUGH_pGH (GH);
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)       0
+#ifdef SQR
+#undef SQR
+#endif
+#define SQR(x)      ((x) * (x))
+
 
-  for (i = 0; i < nOutVals; i++) 
+  for (i = total_outvals = 0; i < num_inarrays; i++)
   {
-    outVals [i] = 0.0;
-    for (j = 0; j < nPoints; j++) 
+    for (j = 0; j < num_outvals; j++, total_outvals++)
     {
-      if (inType == CCTK_VARIABLE_CHAR)
+      switch (intypes[i])
       {
-        outVals [i] += SQR (((CCTK_CHAR *) inArray) [i*nPoints + j]);
-      }
-      else if (inType == CCTK_VARIABLE_INT)
-      {
-        outVals [i] += SQR (((CCTK_INT *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm2, scalar)   norm2 += SQR (scalar)
+          ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_CHAR, outvals[total_outvals]);
+          break;
+
+        case CCTK_VARIABLE_INT:
+          ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT, outvals[total_outvals]);
+          break;
+
 #ifdef CCTK_INT2
-      else if (inType == CCTK_VARIABLE_INT2)
-      {
-        outVals [i] += SQR (((CCTK_INT2 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT2:
+          ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT2, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_INT4
-      else if (inType == CCTK_VARIABLE_INT4)
-      {
-        outVals [i] += SQR (((CCTK_INT4 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT4:
+          ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT4, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_INT8
-      else if (inType == CCTK_VARIABLE_INT8)
-      {
-        outVals [i] += SQR (((CCTK_INT8 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_INT8:
+          ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT8, outvals[total_outvals]);
+          break;
 #endif
-      else if (inType == CCTK_VARIABLE_REAL)
-      {
-        outVals [i] += SQR (((CCTK_REAL *) inArray) [i*nPoints + j]);
-      }
+
+        case CCTK_VARIABLE_REAL:
+          ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
+
 #ifdef CCTK_REAL4
-      else if (inType == CCTK_VARIABLE_REAL4)
-      {
-        outVals [i] += SQR (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL4:
+          ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL4, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_REAL8
-      else if (inType == CCTK_VARIABLE_REAL8)
-      {
-        outVals [i] += SQR (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL8:
+          ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL8, outvals[total_outvals]);
+          break;
 #endif
+
 #ifdef CCTK_REAL16
-      else if (inType == CCTK_VARIABLE_REAL16)
-      {
-        outVals [i] += SQR (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
-      }
+        case CCTK_VARIABLE_REAL16:
+          ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL16, outvals[total_outvals]);
+          break;
 #endif
-      else if (inType == CCTK_VARIABLE_COMPLEX) 
-      {
-        outVals [i] += 
-          ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
-            SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im)                
-            );
-      } 
-      else 
-      {
-        CCTK_WARN (1, "PUGH_Norm2_Array: Unknown variable type");
-        return (1);
+
+        case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm2, scalar)   norm2 += SQR ((scalar).Re) +     \
+                                                      SQR ((scalar).Im)
+          ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i], from, to,iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
+
+        default:
+          CCTK_WARN (1, "PUGH_ReductionNorm2: Unknown variable type");
+          return (-1);
       }
     }
   }
 
-#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
 
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
   if (pughGH->nprocs > 1) 
   {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
 
-    /* outVals[] contains now the local sum */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
     if (proc < 0)
     {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
     }
     else
     {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
                         PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
     }
-
-    free (localOutVals);
-
-    nPoints *= pughGH->nprocs;
+    free (local_outvals);
   }
 
 #endif
 
+  /* finally assign the return value */
   if (proc < 0 || proc == pughGH->myproc)
-    for (i = 0; i < nOutVals; i++)
-      outVals [i] = sqrt (outVals [i] / nPoints);
-
-  return (0);
-}
-
-
-/*@@
-   @routine    PUGH_Norm2_GF
-   @date       Thu Apr  3 11:56:36 1997
-   @author     Paul Walker
-   @desc 
-   Returns the "norm2" of a distributed grid function. That
-   is, returns $\sqrt{\Sigma (a_i * a_i) / np}$.
-   @enddesc 
-@@*/
-int PUGH_Norm2_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
-{
-  int timelevel;
-  pGH *pughGH;
-  pGA *GA;
-  int ii;
-  CCTK_INT tnp;
-#ifdef CCTK_MPI
-  int jj, kk;
-  CCTK_INT np;
-  CCTK_REAL local_result;
-#endif
-
-  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
-
-  pughGH = PUGH_pGH (GH);
-  GA = ((pGA ***) pughGH->variables) [index][timelevel];
-
-  *outVal = 0.0;
-
-#ifndef CCTK_MPI
-  for (ii = 0; ii < GA->extras->npoints; ii++)
   {
-    if (GA->vtype == CCTK_VARIABLE_CHAR)
-    {
-      *outVal += SQR (((CCTK_CHAR *) GA->data) [ii]);
-    }
-    else if (GA->vtype == CCTK_VARIABLE_INT)
-    {
-      *outVal += SQR (((CCTK_INT *) GA->data) [ii]);
-    }
-#ifdef CCTK_INT2
-    else if (GA->vtype == CCTK_VARIABLE_INT2)
-    {
-      *outVal += SQR (((CCTK_INT2 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT4
-    else if (GA->vtype == CCTK_VARIABLE_INT4)
-    {
-      *outVal += SQR (((CCTK_INT4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT8
-    else if (GA->vtype == CCTK_VARIABLE_INT8)
+    for (i = 0; i < total_outvals; i++)
     {
-      *outVal += SQR (((CCTK_INT8 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_REAL)
-    {
-      *outVal += SQR (((CCTK_REAL *) GA->data) [ii]);
-    }
-#ifdef CCTK_REAL4
-    else if (GA->vtype == CCTK_VARIABLE_REAL4)
-    {
-      *outVal += SQR (((CCTK_REAL4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL8
-    else if (GA->vtype == CCTK_VARIABLE_REAL8)
-    {
-      *outVal += SQR (((CCTK_REAL8 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL16
-    else if (GA->vtype == CCTK_VARIABLE_REAL16)
-    {
-      *outVal += SQR (((CCTK_REAL16 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-    {
-      *outVal += 
-        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
-        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im);
-    } 
-    else 
-    {
-      CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
-      return (1);
-    }
-  }
-  tnp = GA->extras->npoints;
-  
-#else
-  /* MPI version depends heavily on the GFExtras->ownership array */
-  np  = 0;
-  local_result = 0.0;
-
-  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
-       kk < GA->extras->ownership[GA->stagger][1][2]; 
-       kk ++)
-  {
-    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
-         jj < GA->extras->ownership[GA->stagger][1][1]; 
-         jj ++)
-    {
-      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
-           ii < GA->extras->ownership[GA->stagger][1][0]; 
-           ii++, np++)
-      {
-        if (GA->vtype == CCTK_VARIABLE_CHAR)
-        {
-          local_result +=
-               SQR (((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-        else if (GA->vtype == CCTK_VARIABLE_INT)
-        {
-          local_result +=
-               SQR (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_INT2
-        else if (GA->vtype == CCTK_VARIABLE_INT2)
-        {
-          local_result +=
-               SQR (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT4
-        else if (GA->vtype == CCTK_VARIABLE_INT4)
-        {
-          local_result += SQR (((CCTK_INT4 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT8
-        else if (GA->vtype == CCTK_VARIABLE_INT8)
-        {
-          local_result += SQR (((CCTK_INT8 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_REAL)
-        {
-          local_result += SQR (((CCTK_REAL *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_REAL4
-        else if (GA->vtype == CCTK_VARIABLE_REAL4)
-        {
-          local_result += SQR (((CCTK_REAL4 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL8
-        else if (GA->vtype == CCTK_VARIABLE_REAL8)
-        {
-          local_result += SQR (((CCTK_REAL8 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL16
-        else if (GA->vtype == CCTK_VARIABLE_REAL16)
-        {
-          local_result +=
-               SQR (((CCTK_REAL16 *) GA->data) 
-		    [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_COMPLEX)
-        {
-          local_result +=
-               SQR((((CCTK_COMPLEX *)GA->data)
-		    [DATINDEX(GA->extras,ii,jj,kk)]).Re)+
-               SQR((((CCTK_COMPLEX *)GA->data)
-		    [DATINDEX(GA->extras,ii,jj,kk)]).Im);        
-	} 
-	else 
-	{
-          CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
-          return (1);
-        }
-      }
+      outvals[i] = sqrt (outvals[i] / num_points);
     }
   }
 
-  if (proc < 0) 
-  {
-    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      pughGH->PUGH_COMM_WORLD));
-  } 
-  else 
-  {
-    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      proc, pughGH->PUGH_COMM_WORLD));
-  }
-#endif
-
-  if (proc < 0 || proc == pughGH->myproc)
-  {
-    *outVal = sqrt (*outVal / tnp);
-  }
-
   return (0);
 }
-
-
diff --git a/src/ReductionSum.c b/src/ReductionSum.c
new file mode 100644
index 0000000..744a6d1
--- /dev/null
+++ b/src/ReductionSum.c
@@ -0,0 +1,325 @@
+ /*@@
+   @file      ReductionSum.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+              Defines the PUGH reduction operator to get the sum
+              of an arbitrary array.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionSum_c)
+
+/* local function prototypes */
+static int PUGH_ReductionSum (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionSumArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the sums
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionSumArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionSum));
+}
+
+
+/*@@
+  @routine PUGH_ReductionSumGVs
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the sum.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           If the number of requested output values equals the size of the
+           variable, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionSumGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionSum));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionSum
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the sum of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+   @enddesc
+@@*/
+static int PUGH_ReductionSum (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, j, total_outvals;
+#ifdef CCTK_MPI
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     0
+#define REDUCTION_OPERATION(sum, scalar)   sum += scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    for (j = 0; j < num_outvals; j++, total_outvals++)
+    {
+      switch (intypes[i])
+      {
+        case CCTK_VARIABLE_CHAR:
+          ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_CHAR, outvals[total_outvals]);
+          break;
+
+        case CCTK_VARIABLE_INT:
+          ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT, outvals[total_outvals]);
+          break;
+
+#ifdef CCTK_INT2
+        case CCTK_VARIABLE_INT2:
+          ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT2, outvals[total_outvals]);
+          break;
+#endif
+
+#ifdef CCTK_INT4
+        case CCTK_VARIABLE_INT4:
+          ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT4, outvals[total_outvals]);
+          break;
+#endif
+
+#ifdef CCTK_INT8
+        case CCTK_VARIABLE_INT8:
+          ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_INT8, outvals[total_outvals]);
+          break;
+#endif
+
+        case CCTK_VARIABLE_REAL:
+          ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL, outvals[total_outvals]);
+          break;
+
+#ifdef CCTK_REAL4
+        case CCTK_VARIABLE_REAL4:
+          ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL4, outvals[total_outvals]);
+          break;
+#endif
+
+#ifdef CCTK_REAL8
+        case CCTK_VARIABLE_REAL8:
+          ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i], from, to, iterator,
+                         points_per_dim, CCTK_REAL8, outvals[total_outvals]);
+          break;
+#endif
+
+#ifdef CCTK_REAL16
+        case CCTK_VARIABLE_REAL16:
+          ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i], from, to,iterator,
+                         points_per_dim, CCTK_REAL16, outvals[total_outvals]);
+          break;
+#endif
+
+        case CCTK_VARIABLE_COMPLEX:
+          CCTK_WARN (1, "PUGH_ReductionSum: Don't know how to compute "
+                        "the sum of complex variables !!!");
+          return (-1);
+
+        default:
+          CCTK_WARN (1, "PUGH_ReductionSum: Unknown variable type");
+          return (-1);
+      }
+    }
+  }
+
+#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
+
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  return (0);
+}
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
index e677991..164b577 100644
--- a/src/pugh_reductions.h
+++ b/src/pugh_reductions.h
@@ -12,24 +12,155 @@
 #define _PUGH_REDUCTIONS_H_
 
 
+#include "cctk.h"
 #include "cctk_Reduction.h"
+#include "pugh.h"
 
-#ifdef _cplusplus
+/***
+   Macro to iterate over every element of an arbitrary sized array
+
+   The array is given by
+     - its CCTK type, cctk_type,
+     - the number of dimensions, vdim,
+     - an (untyped) pointer to its data, vdata,
+     - one-dimensional vectors, from[] and to[], of length vdim,
+       specifying the start- and endpoint to iterate over for each dimension
+
+   An iteration vector, iterator[], is passed in to hold the current indices
+   for each dimension of the array.
+   Another vector of size vdim, points_per_dim[], is used to compute the
+   linear index out of the iterator[] elements.
+
+   Before the iteration loop, a local reduction variable, typed_outval,
+   is defined and initialized by the macro INITIAL_REDUCTION_VALUE(typed_vdata).
+
+   For each iteration the macro REDUCTION_OPERATION(typed_outval, typed_inval)
+   is called, with the local reduction variable and the pointer to the
+   current data element.
+
+   Finally, the local reduction variable, typed_outval, is type-casted
+   to CCTK_REAL and assigned to the output variable, outval.
+
+
+   NOTE: The macro assumes from[i] > to[i] for all 1 <= i < vdim such that
+         at least one iteration will be performed.
+         Make sure this is true when you set up the from[] and to[] vectors.
+
+         For efficiency reasons, the innermost loop (the one for dimension 0)
+         is coded as a for() loop. Thus the overhead to iterate over a
+         one-dimensional array should be kept to a minimum.
+***/
+#define ITERATE_ARRAY(cctk_type, vdim, vdata, from, to, iterator,             \
+                      points_per_dim, outval_type, outval)                    \
+    {                                                                         \
+      int _i, _dim, _vindex;                                                  \
+      cctk_type *typed_vdata = (cctk_type *) (vdata);                         \
+      outval_type typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata);        \
+                                                                              \
+                                                                              \
+      /* set iterator to local startpoint */                                  \
+      memcpy (iterator, from, vdim * sizeof (int));                           \
+                                                                              \
+      /* do the nested loops starting with the second-innermost */            \
+      _dim = 1;                                                               \
+      while (1)                                                               \
+      {                                                                       \
+        /* 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++)                                \
+        {                                                                     \
+          REDUCTION_OPERATION (typed_outval, typed_vdata [_i + _vindex]);     \
+        }                                                                     \
+                                                                              \
+        if (vdim > 1)                                                         \
+        {                                                                     \
+          /* increment current looper and check for end */                    \
+          if (++iterator [_dim] >= to [_dim])                                 \
+          {                                                                   \
+            /* increment outermost loopers */                                 \
+            for (_dim++; _dim < vdim; _dim++)                                 \
+            {                                                                 \
+              if (++iterator [_dim] < to [_dim])                              \
+                break;                                                        \
+            }                                                                 \
+                                                                              \
+            /* done if beyond outermost loop */                               \
+            if (_dim >= vdim)                                                 \
+              break;                                                          \
+                                                                              \
+            /* reset innermost loopers */                                     \
+            for (_dim--; _dim >= 0; _dim--)                                   \
+              iterator [_dim] = from [_dim];                                  \
+            _dim = 1;                                                         \
+          }                                                                   \
+        }                                                                     \
+        else                                                                  \
+        {                                                                     \
+          /* exit loop if array is one-dimensional */                         \
+          break;                                                              \
+        }                                                                     \
+                                                                              \
+      } /* end of nested loops over all dimensions */                         \
+                                                                              \
+      outval = (CCTK_REAL) typed_outval;                                      \
+    }
+
+
+#ifdef __cplusplus
 extern "C" {
 #endif
 
-int PUGH_ReduceMinVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceMaxVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceNorm1(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceNorm2(REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMinValGVs (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValGVs (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumGVs    (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm1GVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm2GVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+
+int PUGH_ReductionMinValArrays (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValArrays (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumArrays    (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm1Arrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm2Arrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+
+typedef int (*reduction_fn_t) (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int from[/* dim */],
+                               int to[/* dim */],
+                               int iterator[/* dim */],
+                               int points_per_dim[/* dim */],
+                               int num_points,
+                               int num_inarrays,
+                               int intypes[/* num_inarrays */],
+                               void *inarrays[/* num_inarrays */],
+                               int num_outvals,
+                               CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
 
-int PUGH_ReduceArrayMinVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceArrayMaxVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceArraySum(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceArrayNorm1(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReduceArrayNorm2(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionGVs (cGH *GH,
+                       int proc,
+                       int num_invars,
+                       int invars[/* num_invars */],
+                       int outtype,
+                       int num_outvals,
+                       void *outvals /* [num_outvals] */,
+                       reduction_fn_t reduction_fn);
+int PUGH_ReductionArrays (cGH *GH,
+                          int proc,
+                          int num_dims,
+                          int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          void *inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn);
 
-#ifdef _cplusplus
+#ifdef __cplusplus
 }
 #endif