Added new reduction operator "count" which returns the number of grid points

of a grid variable which a reduction operator is working on. This closes
PR CactusPUGH/1477.
Added new reduction operators "L3Norm" and "L4Norm" to compute the L3 and L4
norm resp.
Also added synonyms for some existing reduction operators: "average" = "mean",
"norm1" = "L1Norm", "norm2" = "L2Norm", "norm_inf" = "LinfNorm".


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

6 files changed, 1019 insertions, 13 deletions

diff --git a/src/ReductionCount.c b/src/ReductionCount.c
new file mode 100644
index 0000000..ee4c627
--- /dev/null
+++ b/src/ReductionCount.c
@@ -0,0 +1,234 @@
+ /*@@
+   @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 "pugh_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_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 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 PUGH_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 (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionCount));
+}
+
+
+/*@@
+  @routine PUGH_ReductionCountGVs
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the count.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           If the number of requested output values equals the size of the
+           variable, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @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_outvals
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   const int *
+  @vio     in
+  @endvar
+@@*/
+int PUGH_ReductionCountGVs (const cGH *GH,
+                          int proc,
+                          int num_outvals,
+                          int outtype,
+                          void *outvals /* [num_outvals] */,
+                          int num_invars,
+                          const int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             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 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) (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/ReductionNorm3.c b/src/ReductionNorm3.c
new file mode 100644
index 0000000..2ccb67b
--- /dev/null
+++ b/src/ReductionNorm3.c
@@ -0,0 +1,375 @@
+/*@@
+   @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 "pugh_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_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 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 PUGH_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 (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm3));
+}
+
+
+/*@@
+  @routine PUGH_ReductionNorm3GVs
+  @author  Thomas Radke
+  @date    Tue Apr 15 2003
+  @desc
+           PUGH reduction routine to be registered for calculating the norm3.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so num_outvals must be 1.
+  @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     nOutVars
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   const int *
+  @vio     in
+  @endvar
+@@*/
+int PUGH_ReductionNorm3GVs (const cGH *GH,
+                            int proc,
+                            int num_outvals,
+                            int outtype,
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            const int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             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 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 pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* 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_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);
+    }
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  /* finally assign the return value */
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    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..1476eb3
--- /dev/null
+++ b/src/ReductionNorm4.c
@@ -0,0 +1,373 @@
+/*@@
+   @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 "pugh_reductions.h"
+
+static const char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_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 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 PUGH_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 (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                ReductionNorm4));
+}
+
+
+/*@@
+  @routine PUGH_ReductionNorm4GVs
+  @author  Thomas Radke
+  @date    Tue Apr 15 2003
+  @desc
+           PUGH reduction routine to be registered for calculating the norm4.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so num_outvals must be 1.
+  @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     nOutVars
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   const int *
+  @vio     in
+  @endvar
+@@*/
+int PUGH_ReductionNorm4GVs (const cGH *GH,
+                            int proc,
+                            int num_outvals,
+                            int outtype,
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            const int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             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 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 pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* 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_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);
+    }
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  /* finally assign the return value */
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < total_outvals; i++)
+    {
+      outvals[i] = pow (outvals[i] / num_points, 1.0/4.0);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/Startup.c b/src/Startup.c
index 2422b58..e40f2b9 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -38,21 +38,37 @@ int PUGHReduce_Startup(void);
 int PUGHReduce_Startup (void)
 {
   /* Register the reduction operators provided by PUGH */
-  CCTK_RegisterReductionOperator (PUGH_ReductionMinValGVs,  "minimum");
-  CCTK_RegisterReductionOperator (PUGH_ReductionMaxValGVs,  "maximum");
-  CCTK_RegisterReductionOperator (PUGH_ReductionSumGVs,     "sum");
   CCTK_RegisterReductionOperator (PUGH_ReductionAvgGVs,     "average");
+  CCTK_RegisterReductionOperator (PUGH_ReductionAvgGVs,     "mean");
+  CCTK_RegisterReductionOperator (PUGH_ReductionCountGVs,   "count");
+  CCTK_RegisterReductionOperator (PUGH_ReductionMaxValGVs,  "maximum");
+  CCTK_RegisterReductionOperator (PUGH_ReductionMinValGVs,  "minimum");
   CCTK_RegisterReductionOperator (PUGH_ReductionNorm1GVs,   "norm1");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm1GVs,   "L1Norm");
   CCTK_RegisterReductionOperator (PUGH_ReductionNorm2GVs,   "norm2");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm2GVs,   "L2Norm");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm3GVs,   "norm3");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm3GVs,   "L3Norm");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm4GVs,   "norm4");
+  CCTK_RegisterReductionOperator (PUGH_ReductionNorm4GVs,   "L4Norm");
   CCTK_RegisterReductionOperator (PUGH_ReductionNormInfGVs, "norm_inf");
+  CCTK_RegisterReductionOperator (PUGH_ReductionSumGVs,     "sum");
 
-  CCTK_RegisterReductionArrayOperator (PUGH_ReductionMinValArrays,  "minimum");
-  CCTK_RegisterReductionArrayOperator (PUGH_ReductionMaxValArrays,  "maximum");
-  CCTK_RegisterReductionArrayOperator (PUGH_ReductionSumArrays,     "sum");
   CCTK_RegisterReductionArrayOperator (PUGH_ReductionAvgArrays,     "average");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionAvgArrays,     "mean");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionCountArrays,   "count");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionMaxValArrays,  "maximum");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionMinValArrays,  "minimum");
   CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm1Arrays,   "norm1");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm1Arrays,   "L1Norm");
   CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm2Arrays,   "norm2");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm2Arrays,   "L2Norm");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm3Arrays,   "norm3");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm3Arrays,   "L3Norm");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm4Arrays,   "norm4");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionNorm4Arrays,   "L4Norm");
   CCTK_RegisterReductionArrayOperator (PUGH_ReductionNormInfArrays, "norm_inf");
+  CCTK_RegisterReductionArrayOperator (PUGH_ReductionSumArrays,     "sum");
 
   return (0);
 }
diff --git a/src/make.code.defn b/src/make.code.defn
index b53a2d3..7e73b26 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,4 +2,6 @@
 # $Header$
 
 # Source files in this directory
-SRCS = ReductionAvg.c ReductionMax.c ReductionMin.c ReductionNorm1.c ReductionNorm2.c ReductionNormInf.c ReductionSum.c Startup.c Reduction.c
+SRCS = ReductionAvg.c ReductionCount.c ReductionMax.c ReductionMin.c          \
+       ReductionNorm1.c ReductionNorm2.c ReductionNorm3.c ReductionNorm4.c    \
+       ReductionNormInf.c ReductionSum.c Startup.c Reduction.c
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
index 6bb56c1..b65cdfb 100644
--- a/src/pugh_reductions.h
+++ b/src/pugh_reductions.h
@@ -133,21 +133,27 @@
 extern "C" {
 #endif
 
-int PUGH_ReductionMinValGVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReductionMaxValGVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReductionSumGVs     (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 int PUGH_ReductionAvgGVs     (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionCountGVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValGVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMinValGVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 int PUGH_ReductionNorm1GVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 int PUGH_ReductionNorm2GVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm3GVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm4GVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 int PUGH_ReductionNormInfGVs (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumGVs     (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 
-int PUGH_ReductionMinValArrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReductionMaxValArrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-int PUGH_ReductionSumArrays     (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int PUGH_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,