This commit was generated by cvs2svn to compensate for changes in r6, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalReduce/trunk@7 7daa882c-dc44-4453-834e-278d26b18e6a

19 files changed, 2915 insertions, 100 deletions

diff --git a/README b/README
index 70df10f..3f9f117 100644
--- a/README
+++ b/README
@@ -1,9 +1,10 @@
-CVS info   : $Header$
-
 Cactus Code Thorn LocalReduce
-Thorn Author(s)     : Ravi Paruchuri <ravi9@bit.csc.lsu.edu>
-Thorn Maintainer(s) : Ravi Paruchuri <ravi9@bit.csc.lsu.edu>
+Authors    : Ravi Paruchuri (ravi9@cct.lsu.edu),
+             Tom Goodale (goodale@cct.lsu.edu), 
+             Yaakoub El Khamra (yye00@cct.lsu.edu)
+CVS info   : $Header$
 --------------------------------------------------------------------------
 
 Purpose of the thorn:
 
+Local array reduction implementation
diff --git a/interface.ccl b/interface.ccl
index cd376e1..749e788 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,5 +1,4 @@
 # Interface definition for thorn LocalReduce
 # $Header$
 
-implements: LocalInterp
-inherits:
+Implements: Reduce
diff --git a/param.ccl b/param.ccl
index 33d07b6..b895612 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,2 +1,3 @@
 # Parameter definitions for thorn LocalReduce
 # $Header$
+
diff --git a/schedule.ccl b/schedule.ccl
index a1c0150..d00ce58 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,7 +1,7 @@
 # Schedule definitions for thorn LocalReduce
 # $Header$
 
-Schedule LocalReduce_Startup at STARTUP after Driver_Startup
+schedule LocalReduce_Startup at STARTUP
 {
-LANG:C
-} "register LocalReduce's reducation operators"
+  LANG:C
+} "Startup routine"
diff --git a/src/Reduction.c b/src/Reduction.c
index d49a468..656ae96 100644
--- a/src/Reduction.c
+++ b/src/Reduction.c
@@ -1,27 +1,41 @@
  /*@@
    @file      Reduction.c
-   @date      Thu  May 27 13:52:53 2004
-   @author    Ravi Paruchuri
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
    @desc
-              Local Reduction operators
+              Various MPI reduction operators.
    @enddesc
    @version   $Id$
  @@*/
 
 #include <stdlib.h>
 
-
+#include "local_reductions.h"
 
 static const char *rcsid = "$Header$";
 
-CCTK_FILEVERSION(CactusBase_LocalReduce_Reduction_c)
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_Reduction_c)
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+static int ReductionScalar (const cGH *GH, int vindex, 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] */);
+
 
 /********************************************************************
  *********************     External Routines   **********************
  ********************************************************************/
 
  /*@@
-   @routine    LocalReductionArrays
+   @routine    PUGH_ReductionArrays
    @author     Thomas Radke
    @date       19 Aug 1999
    @desc
@@ -31,6 +45,17 @@ CCTK_FILEVERSION(CactusBase_LocalReduce_Reduction_c)
    @enddesc
    @calls      copy_real_to_outtype
 
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      const 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_dims
    @vdesc      number of dimensions of input arrays
    @vtype      int
@@ -84,20 +109,23 @@ CCTK_FILEVERSION(CactusBase_LocalReduce_Reduction_c)
                -2  if <num_outvals> is invalid
    @endreturndesc
 @@*/
-int LocalReduce_ReductionArrays ( int num_dims,
+int ReductionArrays (const cGH *GH,
+                          int proc,
+                          int num_dims,
                           const int dims[/* num_dims */],
-                          int intypes[],
+                          int intype,
                           int num_inarrays,
                           const void *const inarrays[/* num_inarrays */],
-                          int outtypes[],
+                          int outtype,
                           int num_outvals,
-                          void *outvals[] /* [num_outvals] */,
+                          void *outvals /* [num_outvals] */,
                           reduction_fn_t reduction_fn)
 {
   int i, num_points, retval;
   int from[1], to[1], iterator[1], points_per_dim[1];
-  int intypes;
-  CCTK_REAL buffer;
+  int *intypes;
+  CCTK_REAL *buffer;
+
 
   points_per_dim[0] = 1;
   from[0] = 0;
@@ -112,7 +140,7 @@ int LocalReduce_ReductionArrays ( int num_dims,
   /* check for zero-sized arrays */
   if (to[0] <= 0)
   {
-    CCTK_WARN (2, "LocalReduce_ReductionArrays: Cannot reduce zero-sized arrays");
+    CCTK_WARN (2, "PUGH_ReductionArrays: Cannot reduce zero-sized arrays");
     return (-1);
   }
 
@@ -121,7 +149,7 @@ int LocalReduce_ReductionArrays ( int num_dims,
     if (num_outvals != to[0])
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "LocalReduce_ReductionArrays: Don't know how to reduce "
+                  "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);
@@ -129,31 +157,36 @@ int LocalReduce_ReductionArrays ( int num_dims,
     }
     to[0] = 1;
   }
-  num_points = to[0];
+  num_points = to[0] * CCTK_nProcs (GH);
 
-  for(array = 0 ; array < num_inarrays; array++)
+  /* set the array types to intype */
+  /* FIXME: could allow to pass in arrays of different types now !!! */
+  intypes = malloc (num_inarrays * sizeof (int));
+  for (i = 0; i < num_inarrays; i++)
   {
-  
-    /* set the array types to intype */
-    /* FIXME: could allow to pass in arrays of different types now !!! */
-      intype = intypes[array];
+    intypes[i] = intype;
+  }
 
-    /* do the reduction on the input arrays */
-    retval = reduction_fn (num_dims, from, to, iterator, points_per_dim,
-                           num_points, 1, 1, intype, inarrays[array],
-                           1, buffer);
+  buffer = malloc (num_outvals * sizeof (CCTK_REAL));
 
-    if (retval == 0 )
-    {
-      /* type-cast the result to the requested datatype */
-      retval = copy_real_to_outtype (1,
-                                     buffer, outtype[array], outvals[array]);
-    }
+  /* do the reduction on the input arrays */
+  retval = reduction_fn (GH, proc, num_dims, from, to, iterator, points_per_dim,
+                         num_points, 1, num_inarrays, intypes, inarrays,
+                         num_outvals, buffer);
+
+  if (retval == 0 && (proc < 0 || proc == CCTK_MyProc (GH)))
+  {
+    /* type-cast the result to the requested datatype */
+    retval = copy_real_to_outtype (num_inarrays * num_outvals,
+                                   buffer, outtype, outvals);
   }
+
+  free (intypes);
+  free (buffer);
+
   return (retval);
 }
 
-
  /*@@
    @routine    copy_real_to_outtype
    @author     Thomas Radke
@@ -164,7 +197,7 @@ int LocalReduce_ReductionArrays ( int num_dims,
    @enddesc
    @calls      CCTK_VarTypeSize
                CCTK_GroupTypeFromVarI
-               ReductionGA
+               ReductionScalar
                copy_real_to_outtype
 
    @var        num_elems
diff --git a/src/ReductionAvg.c b/src/ReductionAvg.c
index 4569307..2f0c602 100644
--- a/src/ReductionAvg.c
+++ b/src/ReductionAvg.c
@@ -1,9 +1,9 @@
  /*@@
    @file      ReductionAvg.c
-   @date      Thu May  27 14:14:53 2004
-   @author    Ravi Paruchuri
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker, Erik Schnetter
    @desc
-              Defines the reduction operator to get the average
+              Defines the PUGH reduction operator to get the average
               of an arbitrary array.
    @enddesc
    @version   $Id$
@@ -12,32 +12,29 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include "local_reductions.h"
 
 static const char *rcsid = "$Id$";
 
-CCTK_FILEVERSION(CactusBase_LocalReduce_ReductionAvg_c)
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionAvg_c)
 
 /* local function prototypes */
-static int ReductionAvg (int num_dims,
-                         const int from[/* dim */],
-                         const int to[/* dim */],
-                         int iterator[/* dim */],
-                         const int points_per_dim[/* dim */],
-                         int num_points,
-                         int have_local_points,
-                         int num_inarrays,
-                         const int intypes[/* num_inarrays */],
-                         const void *const inarrays[/* num_inarrays */],
-                         int num_outvals,
-                         void *outvals[/*num_inarrays*num_outvals*/]);
+static int ReductionAvg (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[]);
 
 
 /*@@
-  @routine LocalReduce_ReductionAvgArrays
-  @author  Ravi Paruchuri
-  @date    May 27 2004
+  @routine PUGH_ReductionAvgArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
   @desc
-           Registered reduction routine for computing the averages
+           Registered PUGH reduction routine for computing the averages
            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.
@@ -99,24 +96,22 @@ static int ReductionAvg (int num_dims,
   @vio     in
   @endvar
 @@*/
-int LocalReduce_ReductionAvgArrays (const cGH *GH,
-                             int proc,
-                             int num_dims,
-                             const int dims[/* num_dims */],
-                             int num_inarrays,
-                             const void *const inarrays[/* num_inarrays */],
-                             int intype,
-                             int num_outvals,
-                             void *outvals[], /* [num_outvals] */,
-                             int outtype)
+int ReductionAvgArrays (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[])
 {
-  return (LocalReduce_ReductionArrays (num_dims, dims,
-                                intype, num_inarrays, inarrays,
-                                outtype, num_outvals, outvals,
-                                ReductionAvg));
+  return (ReductionArrays (N_dims, operator_handle, 
+                         param_table_handle, N_input_arrays,
+                         input_array_dims, input_array_type_codes,
+                         input_arrays, M_output_numbers,
+                         output_number_type_codes, output_numbers));
 }
 
-
 /*****************************************************************************/
 /*                             local functions                               */
 /*****************************************************************************/
@@ -129,25 +124,22 @@ int LocalReduce_ReductionAvgArrays (const cGH *GH,
                (num_outvals == 1) or on the results of the local reductions.
    @enddesc
 @@*/
-static int ReductionAvg (int num_dims,
-                         const int from[/* dim */],
-                         const int to[/* dim */],
-                         int iterator[/* dim */],
-                         const int points_per_dim[/* dim */],
-                         int num_points,
-                         int have_local_points,
-                         int num_inarrays,
-                         const int intypes[/* num_inarrays */],
-                         const void *const inarrays[/* num_inarrays */],
-                         int num_outvals,
-                         CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+static int ReductionAvg (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[])
 {
   int i, total_outvals;
   const char *vtypename;
-  const pGH *pughGH;
 
 
   /* avoid compiler warnings about unused parameters */
+  (void) (GH + 0);
+  (void) (proc + 0);
   (void) (num_points + 0);
 
 /* macros to complete the ITERATE_ARRAY macro */
@@ -156,7 +148,7 @@ static int ReductionAvg (int num_dims,
 
   for (i = total_outvals = 0; i < num_inarrays; i++)
   {
-    switch (intypes[i])
+    switch (input_array_type_codes[i])
     {
       case CCTK_VARIABLE_CHAR:
         ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
@@ -233,24 +225,25 @@ static int ReductionAvg (int num_dims,
 #endif
 
       default:
-        vtypename = CCTK_VarTypeName (intypes[i]);
+        vtypename = CCTK_VarTypeName (input_array_type_codes[i]);
         if (vtypename && strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 21) == 0)
         {
-          CCTK_WARN (1, "LocalReduce_ReductionAvg: Don't know how to compute "
+          CCTK_WARN (1, "PUGH_ReductionAvg: Don't know how to compute "
                         "the average of complex variables !!!");
         }
         else
         {
-          CCTK_WARN (1, "LocalReduce_ReductionAvg: Unknown variable type");
+          CCTK_WARN (1, "PUGH_ReductionAvg: Unknown variable type");
         }
         return (-1);
     }
   }
 
-    for (i = 0; i < total_outvals; i++)
-    {
-      outvals[i] /= num_points;
-    }
+  /* finally assign the return value */
+  for (i = 0; i < total_outvals; i++)
+  {
+    outvals[i] /= num_points;
+  }
 
   return (0);
 }
diff --git a/src/ReductionCount.c b/src/ReductionCount.c
new file mode 100644
index 0000000..71b99b2
--- /dev/null
+++ b/src/ReductionCount.c
@@ -0,0 +1,172 @@
+ /*@@
+   @file      ReductionCount.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke
+   @desc
+              Defines the PUGH reduction operator "count" to get the number of
+              grid points of an arbitrary array.
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionCount_c)
+
+/* local function prototypes */
+static int ReductionCount (const cGH *GH,
+                         int proc,
+                         int num_dims,
+                         const int from[/* dim */],
+                         const int to[/* dim */],
+                         int iterator[/* dim */],
+                         const int points_per_dim[/* dim */],
+                         int num_points,
+                         int have_local_points,
+                         int num_inarrays,
+                         const int intypes[/* num_inarrays */],
+                         const void *const inarrays[/* num_inarrays */],
+                         int num_outvals,
+                         CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionCountArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the counts
+           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
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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   const int *
+  @vio     in
+  @endvar
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionCountArrays (const cGH *GH,
+                             int proc,
+                             int num_dims,
+                             const int dims[/* num_dims */],
+                             int num_inarrays,
+                             const void *const inarrays[/* num_inarrays */],
+                             int intype,
+                             int num_outvals,
+                             void *outvals /* [num_outvals] */,
+                             int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionCount));
+}
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionCount
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc
+               Returns the number of grid points of a distributed array
+   @enddesc
+@@*/
+static int ReductionCount (const cGH *GH,
+                         int proc,
+                         int num_dims,
+                         const int from[/* dim */],
+                         const int to[/* dim */],
+                         int iterator[/* dim */],
+                         const int points_per_dim[/* dim */],
+                         int num_points,
+                         int have_local_points,
+                         int num_inarrays,
+                         const int intypes[/* num_inarrays */],
+                         const void *const inarrays[/* num_inarrays */],
+                         int num_outvals,
+                         CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i;
+
+
+  /* avoid compiler warnings about unused parameters */
+  (void) (GH + 0);
+  (void) (proc + 0);
+  (void) (num_dims + 0);
+  (void) (from + 0);
+  (void) (to + 0);
+  (void) (iterator + 0);
+  (void) (points_per_dim + 0);
+  (void) (have_local_points + 0);
+  (void) (intypes + 0);
+  (void) (inarrays + 0);
+  (void) (num_outvals + 0);
+
+  /* assign the return value */
+  if (proc < 0 || proc == CCTK_MyProc (GH))
+  {
+    for (i = 0; i < num_inarrays; i++)
+    {
+      outvals[i] = num_points;
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionMax.c b/src/ReductionMax.c
new file mode 100644
index 0000000..fb2aae9
--- /dev/null
+++ b/src/ReductionMax.c
@@ -0,0 +1,263 @@
+ /*@@
+   @file      ReductionMax.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc
+              Defines the PUGH reduction operator to get the maximum
+              for a set of arbitrary arrays.
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionMax_c)
+
+/* local function prototypes */
+static int ReductionMaxVal (const cGH *GH,
+                            int proc,
+                            int num_dims,
+                            const int from[/* dim */],
+                            const int to[/* dim */],
+                            int iterator[/* dim */],
+                            const int points_per_dim[/* dim */],
+                            int num_points,
+                            int have_local_points,
+                            int num_inarrays,
+                            const int intypes[/* num_inarrays */],
+                            const void *const inarrays[/* num_inarrays */],
+                            int num_outvals,
+                            CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionMaxValArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the maxima
+           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
+           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   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_inarrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const *
+  @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 ReductionMaxValArrays (const cGH *GH,
+                                int proc,
+                                int num_dims,
+                                const int dims[/* num_dims */],
+                                int num_inarrays,
+                                const void *const inarrays[/* num_inarrays */],
+                                int intype,
+                                int num_outvals,
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionMaxVal));
+}
+
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionMaxVal
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @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
+@@*/
+static int ReductionMaxVal (const cGH *GH,
+                            int proc,
+                            int num_dims,
+                            const int from[/* dim */],
+                            const int to[/* dim */],
+                            int iterator[/* dim */],
+                            const int points_per_dim[/* dim */],
+                            int num_points,
+                            int have_local_points,
+                            int num_inarrays,
+                            const int intypes[/* num_inarrays */],
+                            const void *const inarrays[/* num_inarrays */],
+                            int num_outvals,
+                            CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+  const char *vtypename;
+#ifdef CCTK_MPI
+  const pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+  /* prevent compiler warnings about unused parameters */
+  (void) (GH + 0);
+  (void) (proc + 0);
+  (void) (num_points + 0);
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     ((array)[0])
+#define REDUCTION_OPERATION(max, scalar)   if (max < scalar) max = scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_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, num_outvals, total_outvals);
+      break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        vtypename = CCTK_VarTypeName (intypes[i]);
+        if (vtypename && strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 21) == 0)
+        {
+          CCTK_WARN (1, "PUGH_ReductionMaxVal: Don't know how to compute "
+                        "the maximum of complex variables !!!");
+        }
+        else
+        {
+          CCTK_WARN (1, "PUGH_ReductionMaxVal: Unknown variable type");
+        }
+        return (-1);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionMin.c b/src/ReductionMin.c
new file mode 100644
index 0000000..8d73384
--- /dev/null
+++ b/src/ReductionMin.c
@@ -0,0 +1,264 @@
+ /*@@
+   @file      ReductionMin.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc
+              Defines the PUGH reduction operator to get the minimum
+              for a set of arbitrary arrays.
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionMin_c)
+
+/* local function prototypes */
+static int ReductionMinVal (const cGH *GH,
+                            int proc,
+                            int num_dims,
+                            const int from[/* dim */],
+                            const int to[/* dim */],
+                            int iterator[/* dim */],
+                            const int points_per_dim[/* dim */],
+                            int num_points,
+                            int have_local_points,
+                            int num_inarrays,
+                            const int intypes[/* num_inarrays */],
+                            const void *const inarrays[/* num_inarrays */],
+                            int num_outvals,
+                            CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionMinValArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the minima
+           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
+           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   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_inarrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const *
+  @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 ReductionMinValArrays (const cGH *GH,
+                                int proc,
+                                int num_dims,
+                                const int dims[/* num_dims */],
+                                int num_inarrays,
+                                const void *const inarrays[/* num_inarrays */],
+                                int intype,
+                                int num_outvals,
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionMinVal));
+}
+
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionMinVal
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @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
+@@*/
+static int ReductionMinVal (const cGH *GH,
+                            int proc,
+                            int num_dims,
+                            const int from[/* dim */],
+                            const int to[/* dim */],
+                            int iterator[/* dim */],
+                            const int points_per_dim[/* dim */],
+                            int num_points,
+                            int have_local_points,
+                            int num_inarrays,
+                            const int intypes[/* num_inarrays */],
+                            const void *const inarrays[/* num_inarrays */],
+                            int num_outvals,
+                            CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+  const char *vtypename;
+#ifdef CCTK_MPI
+  const pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+  /* prevent compiler warnings about unused parameters */
+  (void) (GH + 0);
+  (void) (proc + 0);
+  (void) (num_points + 0);
+
+/* 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++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        vtypename = CCTK_VarTypeName (intypes[i]);
+        if (vtypename && strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 21) == 0)
+        {
+          CCTK_WARN (1, "PUGH_ReductionMinVal: Don't know how to compute "
+                        "the minimum of complex variables !!!");
+        }
+        else
+        {
+          CCTK_WARN (1, "PUGH_ReductionMinVal: Unknown variable type");
+        }
+        return (-1);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
new file mode 100644
index 0000000..81df22a
--- /dev/null
+++ b/src/ReductionNorm1.c
@@ -0,0 +1,295 @@
+ /*@@
+   @file      ReductionNorm1.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm1 = $\Sigma |a_i| / np$
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionNorm1_c)
+
+/* local function prototypes */
+static int ReductionNorm1 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm1Arrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm1"
+           of a set of arrays.
+           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 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
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionNorm1Arrays (const cGH *GH,
+                               int proc,
+                               int num_dims,
+                               const int dims[/* num_dims */],
+                               int num_inarrays,
+                               const void *const inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm1));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionNorm1
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm1" 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.
+
+               "norm1" is defined as $\Sigma |a_i| / np$.
+   @enddesc
+@@*/
+static int ReductionNorm1 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+
+/* 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 = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += scalar
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += CCTK_Cmplx8Abs (scalar)
+        ITERATE_ARRAY (CCTK_COMPLEX8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += CCTK_Cmplx16Abs (scalar)
+        ITERATE_ARRAY (CCTK_COMPLEX16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += CCTK_Cmplx32Abs (scalar)
+        ITERATE_ARRAY (CCTK_COMPLEX32, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm1: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  /* finally assign the return value */
+  for (i = 0; i < total_outvals; i++)
+  {
+    outvals[i] /= num_points;
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNorm2.c b/src/ReductionNorm2.c
new file mode 100644
index 0000000..d2c2093
--- /dev/null
+++ b/src/ReductionNorm2.c
@@ -0,0 +1,290 @@
+/*@@
+   @file      ReductionNorm2.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @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   $Id$
+ @@*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionNorm2_c)
+
+/* local function prototypes */
+static int ReductionNorm2 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm2Arrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm2"
+           of a set of arrays.
+           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 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
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionNorm2Arrays (const cGH *GH,
+                               int proc,
+                               int num_dims,
+                               const int dims[/* num_dims */],
+                               int num_inarrays,
+                               const void *const inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm2));
+}
+
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionNorm2
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @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
+@@*/
+static int ReductionNorm2 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+
+/* 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 = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm2, scalar)   norm2 += SQR (scalar)
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+        ITERATE_ARRAY (CCTK_COMPLEX8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+        ITERATE_ARRAY (CCTK_COMPLEX16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+        ITERATE_ARRAY (CCTK_COMPLEX32, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm2: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  /* finally assign the return value */
+  for (i = 0; i < total_outvals; i++)
+  {
+    outvals[i] = sqrt (outvals[i] / num_points);
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNorm3.c b/src/ReductionNorm3.c
new file mode 100644
index 0000000..366bd61
--- /dev/null
+++ b/src/ReductionNorm3.c
@@ -0,0 +1,288 @@
+/*@@
+   @file      ReductionNorm3.c
+   @date      Tue Apr 15 2003
+   @author    Thomas Radke
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm3 = $\sqrt[3]{\Sigma |a_i|^3 / np}$
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionNorm3_c)
+
+/* local function prototypes */
+static int ReductionNorm3 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm3Arrays
+  @author  Thomas Radke
+  @date    Tue Apr 15 2003
+  @desc
+           Registered PUGH reduction routine for computing the "norm3"
+           of a set of arrays.
+           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 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
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionNorm3Arrays (const cGH *GH,
+                               int proc,
+                               int num_dims,
+                               const int dims[/* num_dims */],
+                               int num_inarrays,
+                               const void *const inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm3));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionNorm3
+   @date       Tue Apr 15 2003
+   @author     Thomas Radke
+   @desc       Returns the "norm3" 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.
+
+               "norm3" is defined as $\sqrt{\Sigma |a_i|^3 / np}$.
+   @enddesc
+@@*/
+static int ReductionNorm3 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)       0
+
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define CUBE_ABS(x)      ((x) * (x) * (x))
+#define REDUCTION_OPERATION(norm3, scalar)   norm3 += CUBE_ABS (scalar)
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+#undef  CUBE_ABS
+#define CUBE_ABS(x)      ((x) < 0 ? -((x) * (x) * (x)) : (x) * (x) * (x))
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm3, scalar)   norm3 += CUBE_ABS ((scalar).Re) +\
+                                                      CUBE_ABS ((scalar).Im)
+        ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+        ITERATE_ARRAY (CCTK_COMPLEX8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+        ITERATE_ARRAY (CCTK_COMPLEX16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+        ITERATE_ARRAY (CCTK_COMPLEX32, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm3: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  /* finally assign the return value */
+  for (i = 0; i < total_outvals; i++)
+  {
+    outvals[i] = pow (outvals[i] / num_points, 1.0 / 3.0);
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNorm4.c b/src/ReductionNorm4.c
new file mode 100644
index 0000000..cf09f60
--- /dev/null
+++ b/src/ReductionNorm4.c
@@ -0,0 +1,286 @@
+/*@@
+   @file      ReductionNorm4.c
+   @date      Tue Apr 15 2003
+   @author    Thomas Radke
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm4 = $\sqrt[4]{\Sigma (a_i^4) / np}$
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionNorm4_c)
+
+/* local function prototypes */
+static int ReductionNorm4 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm4Arrays
+  @author  Thomas Radke
+  @date    Tue Apr 15 2003
+  @desc
+           Registered PUGH reduction routine for computing the "norm4"
+           of a set of arrays.
+           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 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
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionNorm4Arrays (const cGH *GH,
+                               int proc,
+                               int num_dims,
+                               const int dims[/* num_dims */],
+                               int num_inarrays,
+                               const void *const inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm4));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionNorm4
+   @date       Tue Apr 15 2003
+   @author     Thomas Radke
+   @desc       Returns the "norm4" 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.
+
+               "norm4" is defined as $\sqrt{\Sigma (a_i^4) / np}$.
+   @enddesc
+@@*/
+static int ReductionNorm4 (const cGH *GH,
+                           int proc,
+                           int num_dims,
+                           const int from[/* dim */],
+                           const int to[/* dim */],
+                           int iterator[/* dim */],
+                           const int points_per_dim[/* dim */],
+                           int num_points,
+                           int have_local_points,
+                           int num_inarrays,
+                           const int intypes[/* num_inarrays */],
+                           const void *const inarrays[/* num_inarrays */],
+                           int num_outvals,
+                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)       0
+#define POWER4(x)      ((x) * (x) * (x) * (x))
+
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm4, scalar)   norm4 += POWER4 (scalar)
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_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, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm4, scalar)   norm4 += POWER4 ((scalar).Re) +  \
+                                                      POWER4 ((scalar).Im)
+        ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+        ITERATE_ARRAY (CCTK_COMPLEX8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+        ITERATE_ARRAY (CCTK_COMPLEX16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+        ITERATE_ARRAY (CCTK_COMPLEX32, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm4: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  /* finally assign the return value */
+  for (i = 0; i < total_outvals; i++)
+  {
+    outvals[i] = pow (outvals[i] / num_points, 1.0/4.0);
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNormInf.c b/src/ReductionNormInf.c
new file mode 100644
index 0000000..f3a473a
--- /dev/null
+++ b/src/ReductionNormInf.c
@@ -0,0 +1,301 @@
+ /*@@
+   @file      ReductionNormInf.c
+   @date      Sat Aug 25, 2001
+   @author    Thomas Radke, Paul Walker, Erik Schnetter
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm_inf = $\max |a_i| $
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionNormInf_c)
+
+/* local function prototypes */
+static int ReductionNormInf (const cGH *GH,
+                             int proc,
+                             int num_dims,
+                             const int from[/* dim */],
+                             const int to[/* dim */],
+                             int iterator[/* dim */],
+                             const int points_per_dim[/* dim */],
+                             int num_points,
+                             int have_local_points,
+                             int num_inarrays,
+                             const int intypes[/* num_inarrays */],
+                             const void *const inarrays[/* num_inarrays */],
+                             int num_outvals,
+                             CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNormInfArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm_inf"
+           of a set of arrays.
+           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 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
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   const 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_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   const int *
+  @vio     in
+  @endvar
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionNormInfArrays (const cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 const int dims[/* num_dims */],
+                                 int num_inarrays,
+                                 const void *const inarrays[/* num_inarrays */],
+                                 int intype,
+                                 int num_outvals,
+                                 void *outvals /* [num_outvals] */,
+                                 int outtype)
+{
+  return (ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNormInf));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReductionNormInf
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm_inf" 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.
+
+               "norm_inf" is defined as $\max |a_i|$.
+   @enddesc
+@@*/
+static int ReductionNormInf (const cGH *GH,
+                             int proc,
+                             int num_dims,
+                             const int from[/* dim */],
+                             const int to[/* dim */],
+                             int iterator[/* dim */],
+                             const int points_per_dim[/* dim */],
+                             int num_points,
+                             int have_local_points,
+                             int num_inarrays,
+                             const int intypes[/* num_inarrays */],
+                             const void *const inarrays[/* num_inarrays */],
+                             int num_outvals,
+                             CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+
+/* 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))
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(x,y)    ((x) >= (y) ? (x) : (y))
+
+
+  /* avoid compiler warnings about unused parameters */
+#ifndef CCTK_MPI
+  (void) (GH + 0);
+  (void) (proc + 0);
+#endif
+  (void) (num_points + 0);
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf,scalar)
+        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_BYTE, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf, ABS(scalar))
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_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, num_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, num_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, num_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, num_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, num_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, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf, CCTK_CmplxAbs(scalar))
+        ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf, CCTK_Cmplx8Abs(scalar))
+        ITERATE_ARRAY (CCTK_COMPLEX8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf, CCTK_Cmplx16Abs(scalar))
+        ITERATE_ARRAY (CCTK_COMPLEX16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm_inf, scalar)   norm_inf = MAX(norm_inf, CCTK_Cmplx32Abs(scalar))
+        ITERATE_ARRAY (CCTK_COMPLEX32, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNormInf: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionSum.c b/src/ReductionSum.c
new file mode 100644
index 0000000..e04dc36
--- /dev/null
+++ b/src/ReductionSum.c
@@ -0,0 +1,261 @@
+ /*@@
+   @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   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "local_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionSum_c)
+
+/* local function prototypes */
+static int ReductionSum (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[]);
+
+
+/*@@
+  @routine 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   const 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   const int *
+  @vio     in
+  @endvar
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     input_arrays
+  @vdesc   field of input arrays
+  @vtype   const void *const
+  @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 ReductionSumArrays (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[])
+{
+  return (ReductionSum (N_dims, operator_handle, 
+                         param_table_handle, N_input_arrays,
+                         input_array_dims, input_array_type_codes,
+                         input_arrays, M_output_numbers,
+                         output_number_type_codes, output_numbers));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    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 ReductionSum (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[])
+{
+  int i, total_outvals;
+  const char *vtypename;
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     0
+#define REDUCTION_OPERATION(sum, scalar)   sum += scalar
+
+  int from3D[3], to3D[3];
+  int from2D[2], to2D[2];
+  int from1D[1], to1D[1];
+
+  from3D[0]=0;
+  from3D[1]=0;
+  from3D[2]=0;
+ 
+  from2D[0]=0;
+  from2D[1]=0;
+
+  from1D[0]=0;
+
+  to3D[0]=10;
+  to3D[1]=10;
+  to3D[2]=10;
+
+  to2D[0]=10;
+  to2D[1]=10;
+
+  to31[0]=10;
+
+  for (i = total_outvals = 0; i < N_input_arrays; i++)
+  {
+    switch (input_array_type_codes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_BYTE, input_array_dims[i], input_arrays[i],
+                       from3D, to3D, iterator, points_per_dim,
+                       CCTK_BYTE, output_numbers, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        ITERATE_ARRAY (CCTK_INT, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        ITERATE_ARRAY (CCTK_INT1, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT1, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, input_arrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      default:
+        vtypename = CCTK_VarTypeName (input_array_type_codes[i]);
+        if (vtypename && strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 21) == 0)
+        {
+          CCTK_WARN (1, "ReductionSum: Don't know how to compute "
+                        "the sum of complex variables !!!");
+        }
+        else
+        {
+          CCTK_WARN (1, "ReductionSum: Unknown variable type");
+        }
+        return (-1);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/Startup.c b/src/Startup.c
index ff8c0f7..02188c9 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -9,6 +9,7 @@
 @@*/
 
 #include "cctk.h"
+#include "local_reductions.h"
 
 static const char *rcsid="$Id$";
 
@@ -23,7 +24,21 @@ int LocalReduce_Startup (void)
 {
   /* Register the reduction operators */
 
-  CCTK_RegisterReductionArrayOperator (Local_ReductionAvgArrays,     "average");
-
+/*  CCTK_RegisterReductionLocalArrayOperator (ReductionAvgArrays,     "average");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionAvgArrays,     "mean");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionCountArrays,   "count");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionMaxValArrays,  "maximum");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionMinValArrays,  "minimum");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm1Arrays,   "norm1");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm1Arrays,   "L1Norm");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm2Arrays,   "norm2");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm2Arrays,   "L2Norm");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm3Arrays,   "norm3");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm3Arrays,   "L3Norm");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm4Arrays,   "norm4");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNorm4Arrays,   "L4Norm");
+  CCTK_RegisterReductionLocalArrayOperator (ReductionNormInfArrays, "norm_inf");*/
+  CCTK_RegisterReductionLocalArrayOperator (ReductionSumArrays,     "sum");
+  printf("\n success\n");
   return (0);
 }
diff --git a/src/local_reductions.h b/src/local_reductions.h
new file mode 100644
index 0000000..015dbc1
--- /dev/null
+++ b/src/local_reductions.h
@@ -0,0 +1,168 @@
+ /*@@
+   @header    CCTK_reductions.h
+   @date      April 29 1999
+   @author    Gabrielle Allen
+   @desc
+              Prototypes for pugh reduction operators
+   @enddesc
+   @version   $Header$
+ @@*/
+
+#ifndef _LOCAL_REDUCTIONS_H_
+#define _LOCAL_REDUCTIONS_H_
+
+
+#include "cctk.h"
+#include "cctk_Reduction.h"
+/* #include "CactusPUGH/PUGH/src/include/pugh.h"*/
+
+/***
+   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
+
+   The result of the reduction is specified by
+     - the array where to store the results, outvals[]
+     - the type of the intermediate result scalar, outvals_type
+     - the number of scalars to reduce from this input array, num_outvals
+     - the total number of scalars for all input arrays, total_outvals
+       (incremented after each iteration)
+
+   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 current output array value outvals[].
+   The total output scalar counter, total_outvals, is incremented.
+
+
+   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,                                         \
+                      outvals_type, outvals, num_outvals, total_outvals)      \
+    {                                                                         \
+      int _i, _j, _dim, _vindex;                                              \
+      const cctk_type *typed_vdata = (vdata);                                 \
+      outvals_type typed_outval;                                              \
+                                                                              \
+      for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                   \
+      {                                                                     \
+        /* get the linear index of the element to start with */             \
+        _vindex = from[0];                                                  \
+        for (_i = 1; _i < vdim; _i++)                                       \
+          _vindex += from [_i] * points_per_dim [_i];                       \
+        typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);      \
+                                                                            \
+        /* set iterator to local startpoint */                              \
+        memcpy (iterator, from, vdim * sizeof (int));                       \
+                                                                            \
+        /* do the nested loops starting with the second-innermost */        \
+        _dim = 1;                                                           \
+        while (1)                                                           \
+        {                                                                   \
+          /* 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 */                     \
+                                                                            \
+        outvals [total_outvals++] = (CCTK_REAL) typed_outval;               \
+                                                                            \
+      } /* end of loop over num_outvals */                                  \
+    }                                                                       \
+  }
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+/*
+int ReductionAvgArrays     (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionCountArrays   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionMaxValArrays  (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionMinValArrays  (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionNorm1Arrays   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionNorm2Arrays   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionNorm3Arrays   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionNorm4Arrays   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int ReductionNormInfArrays (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);*/
+int ReductionSumArrays     (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+
+typedef int (*reduction_fn_t) (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[]);
+
+int ReductionArrays (int N_dims, int operator_handle, 
+                          int param_table_handle,   int N_input_arrays,
+                          const CCTK_INT input_array_dims[], 
+                          const CCTK_INT input_array_type_codes[],
+                          const void *const input_arrays[],
+                          int M_output_numbers,
+                          const CCTK_INT output_number_type_codes[],
+                          void *const output_numbers[]);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/src/make.code.defn b/src/make.code.defn
index a0ac3c2..5b65f84 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Startup.c Reduction.c ReductionAvg.c 
+SRCS = ReductionSum.c Startup.c Reduction.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
new file mode 100644
index 0000000..7350d72
--- /dev/null
+++ b/src/pugh_reductions.h
@@ -0,0 +1,185 @@
+ /*@@
+   @header    pugh_reductions.h
+   @date      April 29 1999
+   @author    Gabrielle Allen
+   @desc
+              Prototypes for pugh reduction operators
+   @enddesc
+   @version   $Header$
+ @@*/
+
+#ifndef _PUGH_REDUCTIONS_H_
+#define _PUGH_REDUCTIONS_H_
+
+
+#include "cctk.h"
+#include "cctk_Reduction.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+
+/***
+   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
+
+   The result of the reduction is specified by
+     - the array where to store the results, outvals[]
+     - the type of the intermediate result scalar, outvals_type
+     - the number of scalars to reduce from this input array, num_outvals
+     - the total number of scalars for all input arrays, total_outvals
+       (incremented after each iteration)
+
+   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 current output array value outvals[].
+   The total output scalar counter, total_outvals, is incremented.
+
+
+   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,                                         \
+                      outvals_type, outvals, num_outvals, total_outvals)      \
+    {                                                                         \
+      int _i, _j, _dim, _vindex;                                              \
+      const cctk_type *typed_vdata = (vdata);                                 \
+      outvals_type typed_outval;                                              \
+                                                                              \
+                                                                              \
+      if (have_local_points)                                                  \
+      {                                                                       \
+        for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                   \
+        {                                                                     \
+          /* get the linear index of the element to start with */             \
+          _vindex = from[0];                                                  \
+          for (_i = 1; _i < vdim; _i++)                                       \
+            _vindex += from [_i] * points_per_dim [_i];                       \
+          typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);      \
+                                                                              \
+          /* set iterator to local startpoint */                              \
+          memcpy (iterator, from, vdim * sizeof (int));                       \
+                                                                              \
+          /* do the nested loops starting with the second-innermost */        \
+          _dim = 1;                                                           \
+          while (1)                                                           \
+          {                                                                   \
+            /* 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 */                     \
+                                                                              \
+          outvals [total_outvals++] = (CCTK_REAL) typed_outval;               \
+                                                                              \
+        } /* end of loop over num_outvals */                                  \
+      }                                                                       \
+      else                                                                    \
+      {                                                                       \
+        total_outvals += num_outvals;                                         \
+      }                                                                       \
+    }
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int PUGH_ReductionAvgArrays     (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionCountArrays   (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValArrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMinValArrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm1Arrays   (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm2Arrays   (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm3Arrays   (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm4Arrays   (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNormInfArrays (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumArrays     (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+
+typedef int (*reduction_fn_t) (const cGH *GH,
+                               int proc,
+                               int num_dims,
+                               const int from[/* dim */],
+                               const int to[/* dim */],
+                               int iterator[/* dim */],
+                               const int points_per_dim[/* dim */],
+                               int num_points,
+                               int have_local_points,
+                               int num_inarrays,
+                               const int intypes[/* num_inarrays */],
+                               const void *const inarrays[/* num_inarrays */],
+                               int num_outvals,
+                               CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+int PUGH_ReductionArrays (const cGH *GH,
+                          int proc,
+                          int num_dims,
+                          const int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          const void *const inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif