Registered reduction routines return an error code, and split reduction.c

into several files


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

10 files changed, 2940 insertions, 2050 deletions

diff --git a/src/Reduction.c b/src/Reduction.c
index 3cf224c..c0315b8 100644
--- a/src/Reduction.c
+++ b/src/Reduction.c
@@ -11,6 +11,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+
 #include "cctk.h"
 #include "pugh.h"
 
@@ -18,283 +19,14 @@ static char *rcsid = "$Id$";
 
 CCTK_FILEVERSION(CactusPUGH_PUGH_Reduction_c)
 
-/* some useful macros */
-#ifdef MIN
-#undef MIN
-#endif
-#ifdef MAX
-#undef MAX
-#endif
-#ifdef ABS
-#undef ABS
-#endif
-#ifdef SQR
-#undef SQR
-#endif
-#define MIN(x, y)   ((x) < (y) ? (x) : (y))
-#define MAX(x, y)   ((x) > (y) ? (x) : (y))
-#define ABS(x)      ((x) < 0   ?-(x) : (x))
-#define SQR(x)      ((x) * (x))
-
-
 /* local function prototypes */
-int PUGH_Norm1_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
-int PUGH_Norm2_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
-int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals []);
-int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals []);
 int PUGH_Sum_Array (cGH *GH, int proc, int nOutVals, int nPoints,
                     void *inArray, int inType, CCTK_REAL outVals []);
-int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals []);
-int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals []);
 void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
 int getSizeofElement (int varType);
 
 
 /*@@
-  @routine PUGH_ReduceArrayMinVal
-  @author  Thomas Radke
-  @date    19 Aug 1999
-  @desc
-           Registered PUGH reduction routine for computing the minima
-           of a set of arrays.
-           The arrays are described by their dimensions and variable type.
-           For the number of output values only 1 is accepted.
-           Type casting of the result is provided by specifying the
-           requested output datatype. The intermediate reduction value
-           is always computed as a CCTK_REAL value internally.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   int *
-  @vio     in
-  @endvar 
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     inArrays
-  @vdesc   field of input arrays
-  @vtype   void **
-  @vio     in
-  @endvar 
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @vdesc   number of values per output array
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     outVals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
-  @vio     in
-  @endvar 
-  @var     outType
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceArrayMinVal (cGH *GH, int proc, int nDims, int dims [],
-                             int nArrays, void *inArrays [], int inType,
-                             int nOutVals, void *outVals, int outType)
-{
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a %d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-      CCTK_WARN (1, "PUGH_ReduceArrayMinVal: Unsupported return type");
-      return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArrays [i],
-                           inType, buffer) == 0) 
-    {
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-      {
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-      }
-    }
-  }
-  
-  free (buffer);
-}
-
-
-/*@@
-  @routine PUGH_ReduceArrayMaxVal
-  @author  Thomas Radke
-  @date    19 Aug 1999
-  @desc
-           Registered PUGH reduction routine for computing the maxima
-           of a set of arrays.
-           The arrays are described by their dimensions and variable type.
-           For the number of output values only 1 is accepted.
-           Type casting of the result is provided by specifying the
-           requested output datatype. The intermediate reduction value
-           is always computed as a CCTK_REAL value internally.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   int *
-  @vio     in
-  @endvar 
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     inArrays
-  @vdesc   field of input arrays
-  @vtype   void **
-  @vio     in
-  @endvar 
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @vdesc   number of values per output array
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     outVals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
-  @vio     in
-  @endvar 
-  @var     outType
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceArrayMaxVal (cGH *GH, int proc, int nDims, int dims [],
-                             int nArrays, void *inArrays [], int inType,
-                             int nOutVals, void *outVals, int outType)
-{
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a %d-dimensional array "
-                  "with %d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayMaxVal: Unsupported return type");
-    return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                           buffer) == 0) {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-    }
-  }
-
-  free (buffer);
-}
-
-
-/*@@
   @routine PUGH_ReduceArraySum
   @author  Thomas Radke
   @date    19 Aug 1999
@@ -361,10 +93,11 @@ void PUGH_ReduceArrayMaxVal (cGH *GH, int proc, int nDims, int dims [],
   @vio     in
   @endvar 
 @@*/
-void PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
+int PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
                           int nArrays, void *inArrays [], int inType,
                           int nOutVals, void *outVals, int outType)
 {
+  int retval=0;
   int i, elemSize, nPoints;
   CCTK_REAL *buffer;
 
@@ -380,10 +113,11 @@ void PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
     if (nOutVals != nPoints) 
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a %d-dimensional array with "
+                  "PUGH_ReduceArraySum: Don't know how to reduce a "
+		  "%d-dimensional array with "
                   "%d elements to an output array of %d elements",
                   nDims, nPoints, nOutVals);
-      return;
+      return -1;
     }
     nPoints = 1;
   }
@@ -391,7 +125,7 @@ void PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
   if ((elemSize = getSizeofElement (outType)) < 0) 
   {
     CCTK_WARN (1, "PUGH_ReduceArraySum: Unsupported return type");
-    return;
+    return -1;
   }
 
   buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
@@ -409,1105 +143,22 @@ void PUGH_ReduceArraySum (cGH *GH, int proc, int nDims, int dims [],
         copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
                               outType, nOutVals);
       }
+      else
+      {
+	retval = -1;
+      }
     }
   }
 
   free (buffer);
+  return retval;
 }
 
 
-/*@@
-  @routine PUGH_ReduceArrayNorm1
-  @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 nOutVals 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   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   int *
-  @vio     in
-  @endvar 
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     inArrays
-  @vdesc   field of input arrays
-  @vtype   void **
-  @vio     in
-  @endvar 
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @vdesc   number of values per output array
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     outVals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
-  @vio     in
-  @endvar 
-  @var     outType
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
-                            int nArrays, void *inArrays [], int inType,
-                            int nOutVals, void *outVals, int outType)
-{
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "PUGH_ReduceArrayNorm1: Don't know how to reduce "
-                  "a %d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type");
-    return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-
-    /* do the reduction on input array [i] */
-    if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                           buffer) == 0) 
-    {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-    }
-  }
-
-  free (buffer);
-}
-
-
-/*@@
-  @routine PUGH_ReduceArrayNorm2
-  @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 nOutVals 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   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   int *
-  @vio     in
-  @endvar 
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     inArrays
-  @vdesc   field of input arrays
-  @vtype   void **
-  @vio     in
-  @endvar 
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @vdesc   number of values per output array
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     outVals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
-  @vio     in
-  @endvar 
-  @var     outType
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceArrayNorm2 (cGH *GH, int proc, int nDims, int dims [],
-                            int nArrays, void *inArrays [], int inType,
-                            int nOutVals, void *outVals, int outType)
-{
-  int i, elemSize, nPoints;
-  CCTK_REAL *buffer;
-
-
-  nPoints = dims [0];
-  for (i = 1; i < nDims; i++)
-  {
-    nPoints *= dims [i];
-  }
-
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a %d-dimensional array with "
-                  "%d elements to an output array of %d elements",
-                  nDims, nPoints, nOutVals);
-      return;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceArrayNorm2: Unsupported return type");
-    return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < nArrays; i++) 
-  {
-    /* do the reduction on input array [i] */
-    if (PUGH_Norm2_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
-                          buffer) == 0) {
-
-      /* type-cast the result to the requested datatype */
-      if (proc < 0 || proc == CCTK_MyProc (GH))
-        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                              outType, nOutVals);
-    }
-  }
-
-  free (buffer);
-}
-
-
-/*@@
-  @routine PUGH_ReduceMinVal
-  @author  Thomas Radke
-  @date    23 Apr 1999
-  @desc
-           PUGH reduction routine to be registered for calculating the minimum.
-           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 nOutVals must be 1, and global reduction is done on the
-           results of the local reductions.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @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     numinvars
-  @vdesc   number of variables in the list
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     varlist
-  @vdesc   list of variables (given as their global indices)
-  @vtype   int *
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceMinVal (cGH *GH, 
-                        int proc, 
-                        int nOutVals,
-                        int outType, 
-                        void *outVals, 
-                        int numinvars, 
-                        int *varlist)
-{
-
-  int        i;
-  int        elemSize;
-  int        nPoints;
-  pGH       *pughGH;
-  CCTK_REAL *buffer;
-  pGExtras  *extras;
-
-  pughGH = PUGH_pGH (GH);
-
-  /* Get the extras pointer from first variable index */
-  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
-
-  nPoints = extras->npoints;
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints) 
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a GF to an output array "
-                  "of %d elements", nOutVals);
-      return;
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-    CCTK_WARN (1, "PUGH_ReduceMinVal: Unsupported output datatype");
-    return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < numinvars; i++)
-  {
-    int result;
-    int timelevel;
-    void *inArray;
-    int inType;
-
-    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
-    if (timelevel > 0)
-      timelevel--;
-    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
-    inType = CCTK_VarTypeI (varlist [i]);
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
-    {
-      case GROUP_GF:
-        result = PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArray,
-                                    inType, buffer);
-        break;
-
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceMinVal: Reducing arrays not yet supported");
-        break;
-
-      case GROUP_SCALAR:
-        CCTK_WARN(1, "PUGH_ReduceMinVal: Reducing a scalar doesn't make sense");
-        break;
-
-      default:
-        CCTK_WARN (1, "PUGH_ReduceMinVal: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                            outType, nOutVals);
-    }
-  }
-
-  free (buffer);
-}
-
-
-/*@@
-  @routine PUGH_ReduceMaxVal
-  @author  Thomas Radke
-  @date    23 Apr 1999
-  @desc
-           PUGH reduction routine to be registered for calculating the maximum.
-           It just goes through the list of variables and calls the appropriate
-           grouptype reduction routine.
-           If the number of requested output values equals the size of the
-           variable, global reduction is done element-wise.
-           Otherwise nOutVals must be 1, and global reduction is done on the
-           results of the local reductions.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @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     numinvars
-  @vdesc   number of variables in the list
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     varlist
-  @vdesc   list of variables (given as their global indices)
-  @vtype   int *
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceMaxVal (cGH *GH, int proc,  int nOutVals,
-                       int outType, void *outVals, int numinvars, int *varlist)
-{
-  int i, elemSize, nPoints;
-  pGH *pughGH;
-  CCTK_REAL *buffer;
-  pGExtras  *extras;
-
-  pughGH = PUGH_pGH (GH);
-
-  /* Get the extras pointer from first variable index */
-  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
-
-  nPoints = extras->npoints;
-  if (nOutVals != 1) 
-  {
-    if (nOutVals != nPoints)
-    {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Don't know how to reduce a GF to an output array "
-                  "of %d elements", nOutVals);
-    }
-    nPoints = 1;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceMaxVal: Unsupported return type");
-    return;
-  }
-
-  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-  for (i = 0; i < numinvars; i++) 
-  {
-    int result;
-    int timelevel;
-    void *inArray;
-    int inType;
-
-    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
-    if (timelevel > 0)
-    {
-      timelevel--;
-    }
-    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
-    inType = CCTK_VarTypeI (varlist [i]);
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
-    {
-      case GROUP_GF:
-        result = PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArray,
-                                    inType, buffer);
-        break;
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing arrays not yet supported");
-        break;
-      case GROUP_SCALAR:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing a scalar doesn't make sense");
-        break;
-      default:
-        CCTK_WARN (1, "PUGH_ReduceMaxVal: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
-                            outType, nOutVals);
-    }
-  }
-
-  free (buffer);
-}
-
-
-/*@@
-  @routine PUGH_ReduceNorm1
-  @author  Thomas Radke
-  @date    23 Apr 1999
-  @desc
-           PUGH reduction routine to be registered for calculating the norm1.
-           It just goes through the list of variables and calls the appropriate
-           grouptype reduction routine.
-           Global reduction is always done on the results of the local
-           reductions, so nOutVals must be 1.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     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     numinvars
-  @vdesc   number of variables in the list
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     varlist
-  @vdesc   list of variables (given as their global indices)
-  @vtype   int *
-  @vio     in
-  @endvar 
-@@*/
-void PUGH_ReduceNorm1 (cGH *GH, int proc,  
-		       int nOutVals,
-                       int outType, 
-		       void *outVals, 
-		       int numinvars, 
-		       int *varlist)
-{
-  int i, elemSize;
-  pGH *pughGH;
-  CCTK_REAL intermediate_result;
-
-  if (nOutVals != 1) 
-  {
-    CCTK_WARN (1, "PUGH_ReduceNorm1 accepts only one output value per "
-                  "input array");
-    return;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-      CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type");
-      return;
-  }
-
-  pughGH = PUGH_pGH (GH);
-
-  for (i = 0; i < numinvars; i++) 
-  {
-    int result;
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i]))
-    {
-      case GROUP_GF:
-        result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result);
-        break;
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported");
-        break;
-      case GROUP_SCALAR:
-        CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense");
-        break;
-      default:
-        CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (&intermediate_result, (char *) outVals +
-                            i*nOutVals*elemSize, outType, nOutVals);
-    }
-  }
-  
-}
-
-
-/*@@
-  @routine PUGH_ReduceNorm2
-  @author  Thomas Radke
-  @date    23 Apr 1999
-  @desc
-           PUGH reduction routine to be registered for calculating the norm2.
-           It just goes through the list of variables and calls the appropriate
-           grouptype reduction routine.
-           Global reduction is always done on the results of the local
-           reductions, so nOutVals must be 1.
-  @enddesc
-  @history
-  @endhistory 
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   cGH
-  @vio     in
-  @endvar 
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     nOutVals
-  @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     numinvars
-  @vdesc   number of variables in the list
-  @vtype   int
-  @vio     in
-  @endvar 
-  @var     varlist
-  @vdesc   list of variables (given as their global indices)
-  @vtype   int *
-  @vio     in
-  @endvar 
-@@*/
-
-void PUGH_ReduceNorm2 (cGH *GH, int proc, int nOutVals,
-                       int outType, void *outVals, int numinvars, int *varlist)
-{
-  int i, elemSize;
-  pGH *pughGH;
-  CCTK_REAL intermediate_result;
-
-
-  if (nOutVals != 1)
-  {
-    CCTK_WARN (1, "PUGH_ReduceNorm2 accepts only one output value per "
-                  "input array");
-    return;
-  }
-
-  if ((elemSize = getSizeofElement (outType)) < 0)
-  {
-      CCTK_WARN (1, "PUGH_ReduceNorm2: Unsupported return type");
-      return;
-  }
-
-  pughGH = PUGH_pGH (GH);
-
-  for (i = 0; i < numinvars; i++)
-  {
-    int result;
-
-    result = 1;
-
-    switch (CCTK_GroupTypeFromVarI (varlist [i])) {
-      case GROUP_GF:
-        result = PUGH_Norm2_GF (GH, varlist [i], proc, &intermediate_result);
-        break;
-      case GROUP_ARRAY:
-        CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing arrays not yet supported");
-        break;
-      case GROUP_SCALAR:
-        CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing a scalar doesn't make sense");
-        break;
-      default:
-        CCTK_WARN (1, "PUGH_ReduceNorm2: Invalid group type");
-        break;
-    }
-
-    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
-    {
-      copy_real_to_outtype (&intermediate_result, (char *) outVals +
-                            i*nOutVals*elemSize, outType, nOutVals);
-    }
-  }
-}
-
-
-
 /*********************************************************************/
 /*                        local functions                            */
 /*********************************************************************/
 
-/*@@
-   @routine    PUGH_MinVal_Array
-   @date       Aug 19 1999
-   @author     Thomas Radke
-   @desc 
-               Returns the minimum value of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-   @enddesc 
-@@*/
-int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals [])
-{
-  int i, j;
-#ifdef CCTK_MPI
-  pGH *pughGH = PUGH_pGH (GH);
-#endif
-  
-  switch (inType) 
-  {
-    case CCTK_VARIABLE_CHAR:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr)
-          outVals [j] = MIN (*inPtr, (CCTK_CHAR) (outVals [j]));
-      }
-      break;
-
-    case CCTK_VARIABLE_INT:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT) (outVals [j]));
-      }
-      break;
-
-#ifdef CCTK_INT2
-    case CCTK_VARIABLE_INT2:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT2) (outVals [j]));
-      }
-      break;
-#endif
-
-#ifdef CCTK_INT4
-    case CCTK_VARIABLE_INT4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT4) (outVals [j]));
-      }
-      break;
-#endif 
-
-#ifdef CCTK_INT8
-    case CCTK_VARIABLE_INT8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_INT8) (outVals [j]));
-      }
-      break;
-#endif
-
-    case CCTK_VARIABLE_REAL:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, outVals [j]);
-      }
-      break;
-
-#ifdef CCTK_REAL4
-    case CCTK_VARIABLE_REAL4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL4)( outVals [j]));
-      }
-      break;
-#endif
-
-#ifdef CCTK_REAL8
-    case CCTK_VARIABLE_REAL8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL8)( outVals [j]));
-      }
-      break;
-#endif
-
-#ifdef CCTK_REAL16
-    case CCTK_VARIABLE_REAL16:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-          outVals [j] = MIN (*inPtr, (CCTK_REAL16)( outVals [j]));
-      }
-      break;
-#endif
-
-    case CCTK_VARIABLE_COMPLEX:
-      CCTK_WARN(1,"Don't know how to compute minimum of complex variables !!!");
-      return (1);
-
-    default:
-      CCTK_WARN (1, "Unsupported variable type in PUGH_MinVal_Array");
-      return (1);
-  }
-
-
-#ifdef CCTK_MPI
-
-  if (pughGH->nprocs > 1) 
-  {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-    /* outVals[] contains now the local min */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
-    if (proc < 0)
-    {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_MIN, pughGH->PUGH_COMM_WORLD));
-    }
-    else
-    {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_MIN, proc, pughGH->PUGH_COMM_WORLD));
-    }
-    free (localOutVals);
-  }
-
-#endif
-
-  return (0);
-}
-
-
-/*@@
-   @routine    PUGH_MaxVal_Array
-   @date       Aug 19 1999
-   @author     Thomas Radke
-   @desc 
-               Returns the maximum value of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-   @enddesc 
-@@*/
-int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                       void *inArray, int inType, CCTK_REAL outVals [])
-{
-  int i, j;
-#ifdef CCTK_MPI
-  pGH *pughGH = PUGH_pGH (GH);
-#endif
-  
-
-  switch (inType) 
-  {
-    case CCTK_VARIABLE_CHAR:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_CHAR) (outVals [j]));
-        }
-      }
-      break;
-
-    case CCTK_VARIABLE_INT:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT) (outVals [j]));
-        }
-      }
-      break;
-
-#ifdef CCTK_INT2
-    case CCTK_VARIABLE_INT2:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT2) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-#ifdef CCTK_INT4
-    case CCTK_VARIABLE_INT4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT4) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-#ifdef CCTK_INT8
-    case CCTK_VARIABLE_INT8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_INT8) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-    case CCTK_VARIABLE_REAL:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, outVals [j]);
-        }
-      }
-      break;
-
-#ifdef CCTK_REAL4
-    case CCTK_VARIABLE_REAL4:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL4) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-#ifdef CCTK_REAL8
-    case CCTK_VARIABLE_REAL8:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL8) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-#ifdef CCTK_REAL16
-    case CCTK_VARIABLE_REAL16:
-      for (j = 0; j < nOutVals; j++) 
-      {
-        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
-
-        outVals [j] = (CCTK_REAL) *inPtr++;
-        for (i = 1; i < nPoints; i++, inPtr++)
-        {
-          outVals [j] = MAX (*inPtr, (CCTK_REAL16) (outVals [j]));
-        }
-      }
-      break;
-#endif
-
-    case CCTK_VARIABLE_COMPLEX:
-      CCTK_WARN(1,"Don't know how to compute maximum of complex variables !!!");
-      return (1);
-
-    default:
-      CCTK_WARN (1, "Unsupported variable type in PUGH_MaxVal_Array");
-      return (1);
-  }
-
-#ifdef CCTK_MPI
-
-  if (pughGH->nprocs > 1) 
-  {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-    /* outVals[] contains now the local max */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
-    if (proc < 0)
-    {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_MAX, pughGH->PUGH_COMM_WORLD));
-    }
-    else
-    {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_MAX, proc, pughGH->PUGH_COMM_WORLD));
-    }
-
-    free (localOutVals);
-  }
-
-#endif
-
-  return (0);
-}
-
 
 /*@@
    @routine    PUGH_Sum_Array
@@ -1621,683 +272,6 @@ int PUGH_Sum_Array (cGH *GH, int proc, int nOutVals, int nPoints,
 }
 
 
-/*@@
-   @routine    PUGH_Norm1_Array
-   @date       Aug 19 1999
-   @author     Thomas Radke
-   @desc       Returns the "norm1" of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-
-               "norm1" is defined as $\Sigma |a_i| / np$.
-   @enddesc 
-@@*/
-int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals [])
-{
-  int i, j;
-  pGH *pughGH = PUGH_pGH (GH);
-
-
-  for (i = 0; i < nOutVals; i++) 
-  {
-    outVals [i] = 0.0;
-    for (j = 0; j < nPoints; j++) 
-    {
-      if (inType == CCTK_VARIABLE_CHAR)           /* CCTK_CHAR is unsigned */
-      {
-        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
-      }
-      else if (inType == CCTK_VARIABLE_INT)
-      {
-        outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]);
-      }
-#ifdef CCTK_INT2
-      else if (inType == CCTK_VARIABLE_INT2)
-      {
-        outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_INT4
-      else if (inType == CCTK_VARIABLE_INT4)
-      {
-        outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_INT8
-      else if (inType == CCTK_VARIABLE_INT8)
-      {
-        outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]);
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_REAL)
-      {
-        outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]);
-      }
-#ifdef CCTK_REAL4
-      else if (inType == CCTK_VARIABLE_REAL4)
-      {
-        outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_REAL8
-      else if (inType == CCTK_VARIABLE_REAL8)
-      {
-        outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_REAL16
-      else if (inType == CCTK_VARIABLE_REAL16)
-      {
-        outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_COMPLEX) 
-      {
-        outVals [i] += 
-          sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
-                 SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) );
-      } 
-      else 
-      {
-        CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type");
-        return (1);
-      }
-    }
-  }
-
-#ifdef CCTK_MPI
-
-  if (pughGH->nprocs > 1) 
-  {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-    /* outVals[] contains now the local sum */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
-    if (proc < 0)
-    {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    }
-    else
-    {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    }
-    free (localOutVals);
-
-    nPoints *= pughGH->nprocs;
-  }
-
-#endif
-
-  if (proc < 0 || proc == pughGH->myproc)
-  {
-    for (i = 0; i < nOutVals; i++)
-    {
-      outVals [i] /= nPoints;
-    }
-  }
-
-  return (0);
-}
-
-
-/*@@
-   @routine    PUGH_Norm2_Array
-   @date       Aug 19 1999
-   @author     Thomas Radke
-   @desc 
-               Returns the "norm2" of a distributed array with
-               'nPoints' elements. Global reduction is done element-wise
-               (nOutVals == 1) or on the results of the local reductions.
-
-               "norm2" is defined as $\sqrt{\Sigma (a_i * a_i) / np}$.
-   @enddesc 
-@@*/
-int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
-                      void *inArray, int inType, CCTK_REAL outVals [])
-{
-  int i, j;
-  pGH *pughGH;
-
-
-  pughGH = PUGH_pGH (GH);
-
-  for (i = 0; i < nOutVals; i++) 
-  {
-    outVals [i] = 0.0;
-    for (j = 0; j < nPoints; j++) 
-    {
-      if (inType == CCTK_VARIABLE_CHAR)
-      {
-        outVals [i] += SQR (((CCTK_CHAR *) inArray) [i*nPoints + j]);
-      }
-      else if (inType == CCTK_VARIABLE_INT)
-      {
-        outVals [i] += SQR (((CCTK_INT *) inArray) [i*nPoints + j]);
-      }
-#ifdef CCTK_INT2
-      else if (inType == CCTK_VARIABLE_INT2)
-      {
-        outVals [i] += SQR (((CCTK_INT2 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_INT4
-      else if (inType == CCTK_VARIABLE_INT4)
-      {
-        outVals [i] += SQR (((CCTK_INT4 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_INT8
-      else if (inType == CCTK_VARIABLE_INT8)
-      {
-        outVals [i] += SQR (((CCTK_INT8 *) inArray) [i*nPoints + j]);
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_REAL)
-      {
-        outVals [i] += SQR (((CCTK_REAL *) inArray) [i*nPoints + j]);
-      }
-#ifdef CCTK_REAL4
-      else if (inType == CCTK_VARIABLE_REAL4)
-      {
-        outVals [i] += SQR (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_REAL8
-      else if (inType == CCTK_VARIABLE_REAL8)
-      {
-        outVals [i] += SQR (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
-      }
-#endif
-#ifdef CCTK_REAL16
-      else if (inType == CCTK_VARIABLE_REAL16)
-      {
-        outVals [i] += SQR (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
-      }
-#endif
-      else if (inType == CCTK_VARIABLE_COMPLEX) 
-      {
-        outVals [i] += 
-          ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
-            SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im)                
-            );
-      } 
-      else 
-      {
-        CCTK_WARN (1, "PUGH_Norm2_Array: Unknown variable type");
-        return (1);
-      }
-    }
-  }
-
-#ifdef CCTK_MPI
-
-  if (pughGH->nprocs > 1) 
-  {
-    CCTK_REAL *localOutVals;
-
-    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
-
-    /* outVals[] contains now the local sum */
-    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
-    if (proc < 0)
-    {
-      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    }
-    else
-    {
-      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
-                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    }
-
-    free (localOutVals);
-
-    nPoints *= pughGH->nprocs;
-  }
-
-#endif
-
-  if (proc < 0 || proc == pughGH->myproc)
-    for (i = 0; i < nOutVals; i++)
-      outVals [i] = sqrt (outVals [i] / nPoints);
-
-  return (0);
-}
-
-
-/*@@
-   @routine    PUGH_Norm1_GF
-   @date       Thu Apr  3 11:56:05 1997
-   @author     Paul Walker
-   @desc 
-   Returns the average of a distributed grid function.
-   @enddesc 
-@@*/
-int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
-{
-  int timelevel;
-  int ii;
-  CCTK_INT tnp;
-  pGH *pughGH;
-  pGA *GA;
-#ifdef CCTK_MPI
-  int jj, kk;
-  CCTK_INT np;
-  CCTK_REAL local_result;
-#endif
-
-  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
-
-  pughGH = PUGH_pGH (GH);
-  GA = ((pGA ***) pughGH->variables) [index][timelevel];
-
-  *outVal = 0.0;
-
-#ifndef CCTK_MPI
-  for (ii = 0; ii < GA->extras->npoints; ii++)
-  {
-    if (GA->vtype == CCTK_VARIABLE_CHAR)
-    {
-      *outVal += ((CCTK_CHAR *) GA->data) [ii];     /* CCTK_CHAR is unsigned */
-    }
-    else if (GA->vtype == CCTK_VARIABLE_INT)
-    {
-      *outVal += ABS (((CCTK_INT *) GA->data) [ii]);
-    }
-#ifdef CCTK_INT2
-    else if (GA->vtype == CCTK_VARIABLE_INT2)
-    {
-      *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT4
-    else if (GA->vtype == CCTK_VARIABLE_INT4)
-    {
-      *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT8
-    else if (GA->vtype == CCTK_VARIABLE_INT8)
-    {
-      *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_REAL)
-    {
-      *outVal += ABS (((CCTK_REAL *) GA->data) [ii]);
-    }
-#ifdef CCTK_REAL4
-    else if (GA->vtype == CCTK_VARIABLE_REAL4)
-    {
-      *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL8
-    else if (GA->vtype == CCTK_VARIABLE_REAL8)
-    {
-      *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL16
-    else if (GA->vtype == CCTK_VARIABLE_REAL16)
-    {
-      *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-    {
-      *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
-                       SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) );
-    } 
-    else 
-    {
-      CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type");
-      return (1);
-    }
-  }
-
-  tnp = GA->extras->npoints;
-
-#else
-
-  /* MPI version depends heavily on the ownership array */
-  np  = 0;
-  local_result = 0.0;
-
-  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
-       kk < GA->extras->ownership[GA->stagger][1][2]; 
-       kk ++)
-  {
-    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
-         jj < GA->extras->ownership[GA->stagger][1][1]; 
-         jj ++)
-    {
-      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
-           ii < GA->extras->ownership[GA->stagger][1][0]; 
-           ii++, np++)
-      {
-        if (GA->vtype == CCTK_VARIABLE_CHAR)
-        {
-          local_result +=              /* CCTK_CHAR is unsigned */
-	    ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)];
-        }
-        else if (GA->vtype == CCTK_VARIABLE_INT)
-        {
-          local_result +=
-	    ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_INT2
-        else if (GA->vtype == CCTK_VARIABLE_INT2)
-        {
-          local_result +=
-	    ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT4
-        else if (GA->vtype == CCTK_VARIABLE_INT4)
-        {
-          local_result +=
-	    ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT8
-        else if (GA->vtype == CCTK_VARIABLE_INT8)
-        {
-          local_result +=
-	    ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_REAL)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_REAL4
-        else if (GA->vtype == CCTK_VARIABLE_REAL4)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL8
-        else if (GA->vtype == CCTK_VARIABLE_REAL8)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL16
-        else if (GA->vtype == CCTK_VARIABLE_REAL16)
-        {
-          local_result +=
-	    ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-        {
-          local_result +=
-            sqrt( SQR((((CCTK_COMPLEX *)GA->data)
-		       [DATINDEX (GA->extras,ii,jj,kk)]).Re)+
-		  SQR((((CCTK_COMPLEX *)GA->data)
-		       [DATINDEX (GA->extras,ii,jj,kk)]).Im) );  
-        } 
-        else 
-        {
-          CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type");
-          return (1);
-        }
-      }
-    }
-  }
-
-  if (proc < 0) 
-  {
-    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      pughGH->PUGH_COMM_WORLD));
-  } 
-  else 
-  {
-    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      proc, pughGH->PUGH_COMM_WORLD));
-  }
-#endif
- 
-  if (proc < 0 || proc == pughGH->myproc)
-  {
-    *outVal /= tnp;
-  }
-
-  return (0);
-}
-
-
-/*@@
-   @routine    PUGH_Norm2_GF
-   @date       Thu Apr  3 11:56:36 1997
-   @author     Paul Walker
-   @desc 
-   Returns the "norm2" of a distributed grid function. That
-   is, returns $\sqrt{\Sigma (a_i * a_i) / np}$.
-   @enddesc 
-@@*/
-int PUGH_Norm2_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
-{
-  int timelevel;
-  pGH *pughGH;
-  pGA *GA;
-  int ii;
-  CCTK_INT tnp;
-#ifdef CCTK_MPI
-  int jj, kk;
-  CCTK_INT np;
-  CCTK_REAL local_result;
-#endif
-
-  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
-  if (timelevel > 0)
-  {
-    timelevel--;
-  }
-
-  pughGH = PUGH_pGH (GH);
-  GA = ((pGA ***) pughGH->variables) [index][timelevel];
-
-  *outVal = 0.0;
-
-#ifndef CCTK_MPI
-  for (ii = 0; ii < GA->extras->npoints; ii++)
-  {
-    if (GA->vtype == CCTK_VARIABLE_CHAR)
-    {
-      *outVal += SQR (((CCTK_CHAR *) GA->data) [ii]);
-    }
-    else if (GA->vtype == CCTK_VARIABLE_INT)
-    {
-      *outVal += SQR (((CCTK_INT *) GA->data) [ii]);
-    }
-#ifdef CCTK_INT2
-    else if (GA->vtype == CCTK_VARIABLE_INT2)
-    {
-      *outVal += SQR (((CCTK_INT2 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT4
-    else if (GA->vtype == CCTK_VARIABLE_INT4)
-    {
-      *outVal += SQR (((CCTK_INT4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_INT8
-    else if (GA->vtype == CCTK_VARIABLE_INT8)
-    {
-      *outVal += SQR (((CCTK_INT8 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_REAL)
-    {
-      *outVal += SQR (((CCTK_REAL *) GA->data) [ii]);
-    }
-#ifdef CCTK_REAL4
-    else if (GA->vtype == CCTK_VARIABLE_REAL4)
-    {
-      *outVal += SQR (((CCTK_REAL4 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL8
-    else if (GA->vtype == CCTK_VARIABLE_REAL8)
-    {
-      *outVal += SQR (((CCTK_REAL8 *) GA->data) [ii]);
-    }
-#endif
-#ifdef CCTK_REAL16
-    else if (GA->vtype == CCTK_VARIABLE_REAL16)
-    {
-      *outVal += SQR (((CCTK_REAL16 *) GA->data) [ii]);
-    }
-#endif
-    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
-    {
-      *outVal += 
-        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
-        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im);
-    } 
-    else 
-    {
-      CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
-      return (1);
-    }
-  }
-  tnp = GA->extras->npoints;
-  
-#else
-  /* MPI version depends heavily on the GFExtras->ownership array */
-  np  = 0;
-  local_result = 0.0;
-
-  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
-       kk < GA->extras->ownership[GA->stagger][1][2]; 
-       kk ++)
-  {
-    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
-         jj < GA->extras->ownership[GA->stagger][1][1]; 
-         jj ++)
-    {
-      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
-           ii < GA->extras->ownership[GA->stagger][1][0]; 
-           ii++, np++)
-      {
-        if (GA->vtype == CCTK_VARIABLE_CHAR)
-        {
-          local_result +=
-               SQR (((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-        else if (GA->vtype == CCTK_VARIABLE_INT)
-        {
-          local_result +=
-               SQR (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_INT2
-        else if (GA->vtype == CCTK_VARIABLE_INT2)
-        {
-          local_result +=
-               SQR (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT4
-        else if (GA->vtype == CCTK_VARIABLE_INT4)
-        {
-          local_result += SQR (((CCTK_INT4 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_INT8
-        else if (GA->vtype == CCTK_VARIABLE_INT8)
-        {
-          local_result += SQR (((CCTK_INT8 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_REAL)
-        {
-          local_result += SQR (((CCTK_REAL *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#ifdef CCTK_REAL4
-        else if (GA->vtype == CCTK_VARIABLE_REAL4)
-        {
-          local_result += SQR (((CCTK_REAL4 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL8
-        else if (GA->vtype == CCTK_VARIABLE_REAL8)
-        {
-          local_result += SQR (((CCTK_REAL8 *) GA->data) 
-			       [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-#ifdef CCTK_REAL16
-        else if (GA->vtype == CCTK_VARIABLE_REAL16)
-        {
-          local_result +=
-               SQR (((CCTK_REAL16 *) GA->data) 
-		    [DATINDEX (GA->extras, ii, jj, kk)]);
-        }
-#endif
-        else if (GA->vtype == CCTK_VARIABLE_COMPLEX)
-        {
-          local_result +=
-               SQR((((CCTK_COMPLEX *)GA->data)
-		    [DATINDEX(GA->extras,ii,jj,kk)]).Re)+
-               SQR((((CCTK_COMPLEX *)GA->data)
-		    [DATINDEX(GA->extras,ii,jj,kk)]).Im);        
-	} 
-	else 
-	{
-          CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
-          return (1);
-        }
-      }
-    }
-  }
-
-  if (proc < 0) 
-  {
-    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      pughGH->PUGH_COMM_WORLD));
-  } 
-  else 
-  {
-    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
-                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
-    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
-                      proc, pughGH->PUGH_COMM_WORLD));
-  }
-#endif
-
-  if (proc < 0 || proc == pughGH->myproc)
-  {
-    *outVal = sqrt (*outVal / tnp);
-  }
-
-  return (0);
-}
-
 
 void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n)
 {
diff --git a/src/ReductionMax.c b/src/ReductionMax.c
new file mode 100644
index 0000000..7ce59fa
--- /dev/null
+++ b/src/ReductionMax.c
@@ -0,0 +1,497 @@
+ /*@@
+   @file      Reduction.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionMax_c)
+
+/* some useful macros */
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(x, y)   ((x) > (y) ? (x) : (y))
+
+/* local function prototypes */
+int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                       void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+/*@@
+  @routine PUGH_ReduceArrayMaxVal
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the maxima
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceArrayMaxVal (cGH *GH, int proc, int nDims, int dims [],
+                             int nArrays, void *inArrays [], int inType,
+                             int nOutVals, void *outVals, int outType)
+{
+  int retval=1;
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayMaxVal: Don't know how to reduce "
+		  "a %d-dimensional array "
+                  "with %d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return -1;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceArrayMaxVal: Unsupported return type");
+    return -1;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+    /* do the reduction on input array [i] */
+    if (PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
+				    buffer) == 0) 
+    {
+
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+      {
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+	retval = 0;
+      }
+    }
+    else 
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+  
+  return retval;
+}
+
+
+
+/*@@
+  @routine PUGH_ReduceMaxVal
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the maximum.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           If the number of requested output values equals the size of the
+           variable, global reduction is done element-wise.
+           Otherwise nOutVals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @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     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceMaxVal (cGH *GH, int proc,  int nOutVals,
+                       int outType, void *outVals, int numinvars, int *varlist)
+{
+  int i, elemSize, nPoints;
+  int retval=1;
+  pGH *pughGH;
+  CCTK_REAL *buffer;
+  pGExtras  *extras;
+
+  pughGH = PUGH_pGH (GH);
+
+  /* Get the extras pointer from first variable index */
+  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
+
+  nPoints = extras->npoints;
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceMaxVal: Don't know how to reduce a GF to "
+		  "an output array "
+                  "of %d elements", nOutVals);
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceMaxVal: Unsupported return type");
+    return -2;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < numinvars; i++) 
+  {
+    int result;
+    int timelevel;
+    void *inArray;
+    int inType;
+
+    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
+    inType = CCTK_VarTypeI (varlist [i]);
+
+    result = 1;
+
+    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
+    {
+      case GROUP_GF:
+        result = PUGH_MaxVal_Array (GH, proc, nOutVals, nPoints, inArray,
+                                    inType, buffer);
+        break;
+      case GROUP_ARRAY:
+        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing arrays not yet supported");
+        break;
+      case GROUP_SCALAR:
+        CCTK_WARN (1, "PUGH_ReduceMaxVal: Reducing a scalar doesn't make sense");
+        break;
+      default:
+        CCTK_WARN (1, "PUGH_ReduceMaxVal: Invalid group type");
+        break;
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                            outType, nOutVals);
+      retval = 0;
+    }
+    else if (result == 1)
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+
+  return retval;
+}
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+/*@@
+   @routine    PUGH_MaxVal_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc 
+               Returns the maximum value of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+   @enddesc 
+@@*/
+int PUGH_MaxVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                       void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+#ifdef CCTK_MPI
+  pGH *pughGH = PUGH_pGH (GH);
+#endif
+  
+
+  switch (inType) 
+  {
+    case CCTK_VARIABLE_CHAR:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_CHAR) (outVals [j]));
+        }
+      }
+      break;
+
+    case CCTK_VARIABLE_INT:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_INT) (outVals [j]));
+        }
+      }
+      break;
+
+#ifdef CCTK_INT2
+    case CCTK_VARIABLE_INT2:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_INT2) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+#ifdef CCTK_INT4
+    case CCTK_VARIABLE_INT4:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_INT4) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+#ifdef CCTK_INT8
+    case CCTK_VARIABLE_INT8:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_INT8) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+    case CCTK_VARIABLE_REAL:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, outVals [j]);
+        }
+      }
+      break;
+
+#ifdef CCTK_REAL4
+    case CCTK_VARIABLE_REAL4:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_REAL4) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+#ifdef CCTK_REAL8
+    case CCTK_VARIABLE_REAL8:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_REAL8) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+#ifdef CCTK_REAL16
+    case CCTK_VARIABLE_REAL16:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+        {
+          outVals [j] = MAX (*inPtr, (CCTK_REAL16) (outVals [j]));
+        }
+      }
+      break;
+#endif
+
+    case CCTK_VARIABLE_COMPLEX:
+      CCTK_WARN(1,"Don't know how to compute maximum of complex variables !!!");
+      return (1);
+
+    default:
+      CCTK_WARN (1, "Unsupported variable type in PUGH_MaxVal_Array");
+      return (1);
+  }
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local max */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_MAX, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_MAX, proc, pughGH->PUGH_COMM_WORLD));
+    }
+
+    free (localOutVals);
+  }
+
+#endif
+
+  return (0);
+}
+
+
diff --git a/src/ReductionMin.c b/src/ReductionMin.c
new file mode 100644
index 0000000..3e05821
--- /dev/null
+++ b/src/ReductionMin.c
@@ -0,0 +1,492 @@
+ /*@@
+   @file      ReductionMin.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   MPI reduction operator to calculate minimum value
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionMin_c)
+
+/* some useful macros */
+#ifdef MIN
+#undef MIN
+#endif
+#define MIN(x, y)   ((x) < (y) ? (x) : (y))
+
+/* local function prototypes */
+int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                       void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+
+/*@@
+  @routine PUGH_ReduceArrayMinVal
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the minima
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceArrayMinVal (cGH *GH, int proc, int nDims, int dims [],
+                             int nArrays, void *inArrays [], int inType,
+                             int nOutVals, void *outVals, int outType)
+{
+  int retval =1;
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayMinVal: Don't know how to reduce a "
+		  "%d-dimensional array with "
+                  "%d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return -1;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+      CCTK_WARN (1, "PUGH_ReduceArrayMinVal: Unsupported return type");
+      return -1;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+    /* do the reduction on input array [i] */
+    if (PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArrays [i],
+                           inType, buffer) == 0) 
+    {
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+      {
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+	retval = 0;
+      }
+    }
+    else
+    {
+      retval = -1;
+    }
+  }
+  
+  free (buffer);
+
+  return 0;
+
+}
+
+
+
+
+/*@@
+  @routine PUGH_ReduceMinVal
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the minimum.
+           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 nOutVals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @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     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceMinVal (cGH *GH, 
+		       int proc, 
+		       int nOutVals,
+		       int outType, 
+		       void *outVals, 
+		       int numinvars, 
+		       int *varlist)
+{
+
+  int        i;
+  int        retval=1;
+  int        elemSize;
+  int        nPoints;
+  pGH       *pughGH;
+  CCTK_REAL *buffer;
+  pGExtras  *extras;
+
+  pughGH = PUGH_pGH (GH);
+
+  /* Get the extras pointer from first variable index */
+  extras = ((pGA ***) pughGH->variables)[varlist[0]][0]->extras;  
+
+  nPoints = extras->npoints;
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceMinVal: Don't know how to reduce a GF "
+		  "to an output array "
+                  "of %d elements", nOutVals);
+      return -2;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0)
+  {
+    CCTK_WARN (1, "PUGH_ReduceMinVal: Unsupported output datatype");
+    return -3;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < numinvars; i++)
+  {
+    int result;
+    int timelevel;
+    void *inArray;
+    int inType;
+
+    timelevel = CCTK_NumTimeLevelsFromVarI (varlist [i]) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    inArray = CCTK_VarDataPtrI (GH, timelevel, varlist [i]);
+    inType = CCTK_VarTypeI (varlist [i]);
+
+    result = 1;
+
+    switch (CCTK_GroupTypeFromVarI (varlist [i])) 
+    {
+      case GROUP_GF:
+        result = PUGH_MinVal_Array (GH, proc, nOutVals, nPoints, inArray,
+                                    inType, buffer);
+        break;
+
+      case GROUP_ARRAY:
+        CCTK_WARN (1, "PUGH_ReduceMinVal: Reducing arrays not yet supported");
+        break;
+
+      case GROUP_SCALAR:
+        CCTK_WARN(1, "PUGH_ReduceMinVal: Reducing a scalar doesn't make sense");
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReduceMinVal: Invalid group type");
+        break;
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                            outType, nOutVals);
+      retval = 0;
+    }
+    else if (result == 1)
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+
+  return retval;
+}
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+/*@@
+   @routine    PUGH_MinVal_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc 
+               Returns the minimum value of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+   @enddesc 
+@@*/
+int PUGH_MinVal_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                       void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+#ifdef CCTK_MPI
+  pGH *pughGH = PUGH_pGH (GH);
+#endif
+  
+  switch (inType) 
+  {
+    case CCTK_VARIABLE_CHAR:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_CHAR *inPtr = (CCTK_CHAR *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr)
+          outVals [j] = MIN (*inPtr, (CCTK_CHAR) (outVals [j]));
+      }
+      break;
+
+    case CCTK_VARIABLE_INT:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT *inPtr = (CCTK_INT *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_INT) (outVals [j]));
+      }
+      break;
+
+#ifdef CCTK_INT2
+    case CCTK_VARIABLE_INT2:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT2 *inPtr = (CCTK_INT2 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_INT2) (outVals [j]));
+      }
+      break;
+#endif
+
+#ifdef CCTK_INT4
+    case CCTK_VARIABLE_INT4:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT4 *inPtr = (CCTK_INT4 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_INT4) (outVals [j]));
+      }
+      break;
+#endif 
+
+#ifdef CCTK_INT8
+    case CCTK_VARIABLE_INT8:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_INT8 *inPtr = (CCTK_INT8 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_INT8) (outVals [j]));
+      }
+      break;
+#endif
+
+    case CCTK_VARIABLE_REAL:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL *inPtr = (CCTK_REAL *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, outVals [j]);
+      }
+      break;
+
+#ifdef CCTK_REAL4
+    case CCTK_VARIABLE_REAL4:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL4 *inPtr = (CCTK_REAL4 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_REAL4)( outVals [j]));
+      }
+      break;
+#endif
+
+#ifdef CCTK_REAL8
+    case CCTK_VARIABLE_REAL8:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL8 *inPtr = (CCTK_REAL8 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_REAL8)( outVals [j]));
+      }
+      break;
+#endif
+
+#ifdef CCTK_REAL16
+    case CCTK_VARIABLE_REAL16:
+      for (j = 0; j < nOutVals; j++) 
+      {
+        CCTK_REAL16 *inPtr = (CCTK_REAL16 *) inArray + j;
+
+        outVals [j] = (CCTK_REAL) *inPtr++;
+        for (i = 1; i < nPoints; i++, inPtr++)
+          outVals [j] = MIN (*inPtr, (CCTK_REAL16)( outVals [j]));
+      }
+      break;
+#endif
+
+    case CCTK_VARIABLE_COMPLEX:
+      CCTK_WARN(1,"Don't know how to compute minimum of complex variables !!!");
+      return (1);
+
+    default:
+      CCTK_WARN (1, "Unsupported variable type in PUGH_MinVal_Array");
+      return (1);
+  }
+
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local min */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_MIN, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_MIN, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (localOutVals);
+  }
+
+#endif
+
+  return (0);
+}
+
+
diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
new file mode 100644
index 0000000..94a9d3f
--- /dev/null
+++ b/src/ReductionNorm1.c
@@ -0,0 +1,650 @@
+ /*@@
+   @file      Reduction_nm1.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm1_c)
+
+/* some useful macros */
+#ifdef MIN
+#undef MIN
+#endif
+#ifdef ABS
+#undef ABS
+#endif
+#ifdef SQR
+#undef SQR
+#endif
+#define MIN(x, y)   ((x) < (y) ? (x) : (y))
+#define ABS(x)      ((x) < 0   ?-(x) : (x))
+#define SQR(x)      ((x) * (x))
+
+
+/* local function prototypes */
+int PUGH_Norm1_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+
+/*@@
+  @routine PUGH_ReduceArrayNorm1
+  @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 nOutVals 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   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
+                            int nArrays, void *inArrays [], int inType,
+                            int nOutVals, void *outVals, int outType)
+{
+  int retval = 1;
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayNorm1: Don't know how to reduce "
+                  "a %d-dimensional array with "
+                  "%d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return -1;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type");
+    return -1;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+
+    /* do the reduction on input array [i] */
+    if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
+                           buffer) == 0) 
+    {
+
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+      {
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+	retval = 0;
+      }
+    }
+    else
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+  return retval;
+}
+
+
+
+
+/*@@
+  @routine PUGH_ReduceNorm1
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the norm1.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so nOutVals must be 1.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     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     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceNorm1 (cGH *GH, int proc,  
+		       int nOutVals,
+                       int outType, 
+		       void *outVals, 
+		       int numinvars, 
+		       int *varlist)
+{
+  int retval=1;
+  int i, elemSize;
+  pGH *pughGH;
+  CCTK_REAL intermediate_result;
+
+  if (nOutVals != 1) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceNorm1: Only one output value per "
+                  "input array allowed");
+    return -2;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0)
+  {
+      CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type");
+      return -3;
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+  for (i = 0; i < numinvars; i++) 
+  {
+    int result;
+
+    result = 1;
+
+    if (CCTK_GroupDimFromVarI(varlist[i])<3)
+    {
+      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
+		 "PUGH_ReduceNorm1: Reducing dimension %d not yet supported",
+		 CCTK_GroupDimFromVarI(varlist[i]));
+    }
+    else
+    {
+		 
+      switch (CCTK_GroupTypeFromVarI (varlist [i]))
+      {
+        case GROUP_GF:
+	  result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result);
+	  break;
+        case GROUP_ARRAY:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported");
+	  break;
+        case GROUP_SCALAR:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense");
+	  break;
+        default:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type");
+	  break;
+      }
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (&intermediate_result, (char *) outVals +
+                            i*nOutVals*elemSize, outType, nOutVals);
+      retval = 0;
+    }
+    else if (result == 0)
+    {
+      retval = -1;
+    }
+  }
+
+  return retval;
+}
+
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+
+/*@@
+   @routine    PUGH_Norm1_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm1" of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+
+               "norm1" is defined as $\Sigma |a_i| / np$.
+   @enddesc 
+@@*/
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+  pGH *pughGH = PUGH_pGH (GH);
+
+
+  for (i = 0; i < nOutVals; i++) 
+  {
+    outVals [i] = 0.0;
+    for (j = 0; j < nPoints; j++) 
+    {
+      if (inType == CCTK_VARIABLE_CHAR)           /* CCTK_CHAR is unsigned */
+      {
+        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
+      }
+      else if (inType == CCTK_VARIABLE_INT)
+      {
+        outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_INT2
+      else if (inType == CCTK_VARIABLE_INT2)
+      {
+        outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT4
+      else if (inType == CCTK_VARIABLE_INT4)
+      {
+        outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT8
+      else if (inType == CCTK_VARIABLE_INT8)
+      {
+        outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_REAL)
+      {
+        outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_REAL4
+      else if (inType == CCTK_VARIABLE_REAL4)
+      {
+        outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (inType == CCTK_VARIABLE_REAL8)
+      {
+        outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (inType == CCTK_VARIABLE_REAL16)
+      {
+        outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_COMPLEX) 
+      {
+        outVals [i] += 
+          sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
+                 SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) );
+      } 
+      else 
+      {
+        CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type");
+        return (1);
+      }
+    }
+  }
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local sum */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (localOutVals);
+
+    nPoints *= pughGH->nprocs;
+  }
+
+#endif
+
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < nOutVals; i++)
+    {
+      outVals [i] /= nPoints;
+    }
+  }
+
+  return (0);
+}
+
+
+
+/*@@
+   @routine    PUGH_Norm1_GF
+   @date       Thu Apr  3 11:56:05 1997
+   @author     Paul Walker
+   @desc 
+   Returns the average of a distributed grid function.
+   @enddesc 
+@@*/
+int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
+{
+  int timelevel;
+  int ii;
+  CCTK_INT tnp;
+  pGH *pughGH;
+  pGA *GA;
+#ifdef CCTK_MPI
+  int jj, kk;
+  CCTK_INT np;
+  CCTK_REAL local_result;
+#endif
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+
+  pughGH = PUGH_pGH (GH);
+  GA = ((pGA ***) pughGH->variables) [index][timelevel];
+
+  *outVal = 0.0;
+
+#ifndef CCTK_MPI
+  for (ii = 0; ii < GA->extras->npoints; ii++)
+  {
+    if (GA->vtype == CCTK_VARIABLE_CHAR)
+    {
+      *outVal += ((CCTK_CHAR *) GA->data) [ii];     /* CCTK_CHAR is unsigned */
+    }
+    else if (GA->vtype == CCTK_VARIABLE_INT)
+    {
+      *outVal += ABS (((CCTK_INT *) GA->data) [ii]);
+    }
+#ifdef CCTK_INT2
+    else if (GA->vtype == CCTK_VARIABLE_INT2)
+    {
+      *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT4
+    else if (GA->vtype == CCTK_VARIABLE_INT4)
+    {
+      *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT8
+    else if (GA->vtype == CCTK_VARIABLE_INT8)
+    {
+      *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_REAL)
+    {
+      *outVal += ABS (((CCTK_REAL *) GA->data) [ii]);
+    }
+#ifdef CCTK_REAL4
+    else if (GA->vtype == CCTK_VARIABLE_REAL4)
+    {
+      *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL8
+    else if (GA->vtype == CCTK_VARIABLE_REAL8)
+    {
+      *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL16
+    else if (GA->vtype == CCTK_VARIABLE_REAL16)
+    {
+      *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+    {
+      *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
+                       SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) );
+    } 
+    else 
+    {
+      CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type");
+      return (1);
+    }
+  }
+
+  tnp = GA->extras->npoints;
+
+#else
+
+  /* MPI version depends heavily on the ownership array */
+  np  = 0;
+  local_result = 0.0;
+
+  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
+       kk < GA->extras->ownership[GA->stagger][1][2]; 
+       kk ++)
+  {
+    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
+         jj < GA->extras->ownership[GA->stagger][1][1]; 
+         jj ++)
+    {
+      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
+           ii < GA->extras->ownership[GA->stagger][1][0]; 
+           ii++, np++)
+      {
+        if (GA->vtype == CCTK_VARIABLE_CHAR)
+        {
+          local_result +=              /* CCTK_CHAR is unsigned */
+	    ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)];
+        }
+        else if (GA->vtype == CCTK_VARIABLE_INT)
+        {
+          local_result +=
+	    ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_INT2
+        else if (GA->vtype == CCTK_VARIABLE_INT2)
+        {
+          local_result +=
+	    ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT4
+        else if (GA->vtype == CCTK_VARIABLE_INT4)
+        {
+          local_result +=
+	    ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT8
+        else if (GA->vtype == CCTK_VARIABLE_INT8)
+        {
+          local_result +=
+	    ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_REAL)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_REAL4
+        else if (GA->vtype == CCTK_VARIABLE_REAL4)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (GA->vtype == CCTK_VARIABLE_REAL8)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (GA->vtype == CCTK_VARIABLE_REAL16)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+        {
+          local_result +=
+            sqrt( SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Re)+
+		  SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Im) );  
+        } 
+        else 
+        {
+          CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type");
+          return (1);
+        }
+      }
+    }
+  }
+
+  if (proc < 0) 
+  {
+    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      pughGH->PUGH_COMM_WORLD));
+  } 
+  else 
+  {
+    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      proc, pughGH->PUGH_COMM_WORLD));
+  }
+#endif
+ 
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    *outVal /= tnp;
+  }
+
+  return (0);
+}
+
+
diff --git a/src/ReductionNorm2.c b/src/ReductionNorm2.c
new file mode 100644
index 0000000..797f8f5
--- /dev/null
+++ b/src/ReductionNorm2.c
@@ -0,0 +1,636 @@
+ /*@@
+   @file      ReductionNorm2.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm2_c)
+
+/* some useful macros */
+#ifdef SQR
+#undef SQR
+#endif
+#define SQR(x)      ((x) * (x))
+
+
+/* local function prototypes */
+int PUGH_Norm2_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
+int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+
+/*@@
+  @routine PUGH_ReduceArrayNorm2
+  @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 nOutVals 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   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceArrayNorm2 (cGH *GH, int proc, int nDims, int dims [],
+                            int nArrays, void *inArrays [], int inType,
+                            int nOutVals, void *outVals, int outType)
+{
+  int retval = 1;
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayNorm2: Don't know how to reduce a "
+		  "%d-dimensional array with "
+                  "%d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return -1;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceArrayNorm2: Unsupported return type");
+    return -1;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+    /* do the reduction on input array [i] */
+    if (PUGH_Norm2_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
+                          buffer) == 0) 
+    {
+
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+      {
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+	retval = 0;
+      }
+    }
+    else
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+
+  return retval;
+}
+
+
+/*@@
+  @routine PUGH_ReduceNorm2
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the norm2.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so nOutVals must be 1.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @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     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+
+int PUGH_ReduceNorm2 (cGH *GH, int proc, int nOutVals,
+                       int outType, void *outVals, int numinvars, int *varlist)
+{
+  int retval = 1;
+  int i, elemSize;
+  pGH *pughGH;
+  CCTK_REAL intermediate_result;
+
+
+  if (nOutVals != 1)
+  {
+    CCTK_WARN (1, "PUGH_ReduceNorm2 accepts only one output value per "
+                  "input array");
+    return -2;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0)
+  {
+      CCTK_WARN (1, "PUGH_ReduceNorm2: Unsupported return type");
+      return -3;
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+  for (i = 0; i < numinvars; i++)
+  {
+    int result;
+
+    result = 1;
+
+    if (CCTK_GroupDimFromVarI(varlist[i])<3)
+    {
+      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
+		 "PUGH_ReduceNorm2: Reducing dimension %d not yet supported",
+		 CCTK_GroupDimFromVarI(varlist[i]));
+    }
+    else
+    {
+      switch (CCTK_GroupTypeFromVarI (varlist [i])) {
+        case GROUP_GF:
+	  result = PUGH_Norm2_GF (GH, varlist [i], proc, &intermediate_result);
+	  break;
+        case GROUP_ARRAY:
+	  CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing arrays not yet supported");
+	  break;
+        case GROUP_SCALAR:
+	  CCTK_WARN (1, "PUGH_ReduceNorm2: Reducing a scalar doesn't make sense");
+	  break;
+        default:
+	  CCTK_WARN (1, "PUGH_ReduceNorm2: Invalid group type");
+	  break;
+      }
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (&intermediate_result, (char *) outVals +
+                            i*nOutVals*elemSize, outType, nOutVals);
+      retval = 0;
+    }
+    else if (result == 1)
+    {
+      retval = -1;
+    }
+  }
+
+  return retval;
+}
+
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+/*@@
+   @routine    PUGH_Norm2_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc 
+               Returns the "norm2" of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+
+               "norm2" is defined as $\sqrt{\Sigma (a_i * a_i) / np}$.
+   @enddesc 
+@@*/
+int PUGH_Norm2_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+  pGH *pughGH;
+
+
+  pughGH = PUGH_pGH (GH);
+
+  for (i = 0; i < nOutVals; i++) 
+  {
+    outVals [i] = 0.0;
+    for (j = 0; j < nPoints; j++) 
+    {
+      if (inType == CCTK_VARIABLE_CHAR)
+      {
+        outVals [i] += SQR (((CCTK_CHAR *) inArray) [i*nPoints + j]);
+      }
+      else if (inType == CCTK_VARIABLE_INT)
+      {
+        outVals [i] += SQR (((CCTK_INT *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_INT2
+      else if (inType == CCTK_VARIABLE_INT2)
+      {
+        outVals [i] += SQR (((CCTK_INT2 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT4
+      else if (inType == CCTK_VARIABLE_INT4)
+      {
+        outVals [i] += SQR (((CCTK_INT4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT8
+      else if (inType == CCTK_VARIABLE_INT8)
+      {
+        outVals [i] += SQR (((CCTK_INT8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_REAL)
+      {
+        outVals [i] += SQR (((CCTK_REAL *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_REAL4
+      else if (inType == CCTK_VARIABLE_REAL4)
+      {
+        outVals [i] += SQR (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (inType == CCTK_VARIABLE_REAL8)
+      {
+        outVals [i] += SQR (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (inType == CCTK_VARIABLE_REAL16)
+      {
+        outVals [i] += SQR (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_COMPLEX) 
+      {
+        outVals [i] += 
+          ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
+            SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im)                
+            );
+      } 
+      else 
+      {
+        CCTK_WARN (1, "PUGH_Norm2_Array: Unknown variable type");
+        return (1);
+      }
+    }
+  }
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local sum */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+
+    free (localOutVals);
+
+    nPoints *= pughGH->nprocs;
+  }
+
+#endif
+
+  if (proc < 0 || proc == pughGH->myproc)
+    for (i = 0; i < nOutVals; i++)
+      outVals [i] = sqrt (outVals [i] / nPoints);
+
+  return (0);
+}
+
+
+/*@@
+   @routine    PUGH_Norm2_GF
+   @date       Thu Apr  3 11:56:36 1997
+   @author     Paul Walker
+   @desc 
+   Returns the "norm2" of a distributed grid function. That
+   is, returns $\sqrt{\Sigma (a_i * a_i) / np}$.
+   @enddesc 
+@@*/
+int PUGH_Norm2_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
+{
+  int timelevel;
+  pGH *pughGH;
+  pGA *GA;
+  int ii;
+  CCTK_INT tnp;
+#ifdef CCTK_MPI
+  int jj, kk;
+  CCTK_INT np;
+  CCTK_REAL local_result;
+#endif
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+
+  pughGH = PUGH_pGH (GH);
+  GA = ((pGA ***) pughGH->variables) [index][timelevel];
+
+  *outVal = 0.0;
+
+#ifndef CCTK_MPI
+  for (ii = 0; ii < GA->extras->npoints; ii++)
+  {
+    if (GA->vtype == CCTK_VARIABLE_CHAR)
+    {
+      *outVal += SQR (((CCTK_CHAR *) GA->data) [ii]);
+    }
+    else if (GA->vtype == CCTK_VARIABLE_INT)
+    {
+      *outVal += SQR (((CCTK_INT *) GA->data) [ii]);
+    }
+#ifdef CCTK_INT2
+    else if (GA->vtype == CCTK_VARIABLE_INT2)
+    {
+      *outVal += SQR (((CCTK_INT2 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT4
+    else if (GA->vtype == CCTK_VARIABLE_INT4)
+    {
+      *outVal += SQR (((CCTK_INT4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT8
+    else if (GA->vtype == CCTK_VARIABLE_INT8)
+    {
+      *outVal += SQR (((CCTK_INT8 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_REAL)
+    {
+      *outVal += SQR (((CCTK_REAL *) GA->data) [ii]);
+    }
+#ifdef CCTK_REAL4
+    else if (GA->vtype == CCTK_VARIABLE_REAL4)
+    {
+      *outVal += SQR (((CCTK_REAL4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL8
+    else if (GA->vtype == CCTK_VARIABLE_REAL8)
+    {
+      *outVal += SQR (((CCTK_REAL8 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL16
+    else if (GA->vtype == CCTK_VARIABLE_REAL16)
+    {
+      *outVal += SQR (((CCTK_REAL16 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+    {
+      *outVal += 
+        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
+        SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im);
+    } 
+    else 
+    {
+      CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
+      return (1);
+    }
+  }
+  tnp = GA->extras->npoints;
+  
+#else
+  /* MPI version depends heavily on the GFExtras->ownership array */
+  np  = 0;
+  local_result = 0.0;
+
+  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
+       kk < GA->extras->ownership[GA->stagger][1][2]; 
+       kk ++)
+  {
+    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
+         jj < GA->extras->ownership[GA->stagger][1][1]; 
+         jj ++)
+    {
+      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
+           ii < GA->extras->ownership[GA->stagger][1][0]; 
+           ii++, np++)
+      {
+        if (GA->vtype == CCTK_VARIABLE_CHAR)
+        {
+          local_result +=
+               SQR (((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+        else if (GA->vtype == CCTK_VARIABLE_INT)
+        {
+          local_result +=
+               SQR (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_INT2
+        else if (GA->vtype == CCTK_VARIABLE_INT2)
+        {
+          local_result +=
+               SQR (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT4
+        else if (GA->vtype == CCTK_VARIABLE_INT4)
+        {
+          local_result += SQR (((CCTK_INT4 *) GA->data) 
+			       [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT8
+        else if (GA->vtype == CCTK_VARIABLE_INT8)
+        {
+          local_result += SQR (((CCTK_INT8 *) GA->data) 
+			       [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_REAL)
+        {
+          local_result += SQR (((CCTK_REAL *) GA->data) 
+			       [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_REAL4
+        else if (GA->vtype == CCTK_VARIABLE_REAL4)
+        {
+          local_result += SQR (((CCTK_REAL4 *) GA->data) 
+			       [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (GA->vtype == CCTK_VARIABLE_REAL8)
+        {
+          local_result += SQR (((CCTK_REAL8 *) GA->data) 
+			       [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (GA->vtype == CCTK_VARIABLE_REAL16)
+        {
+          local_result +=
+               SQR (((CCTK_REAL16 *) GA->data) 
+		    [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_COMPLEX)
+        {
+          local_result +=
+               SQR((((CCTK_COMPLEX *)GA->data)
+		    [DATINDEX(GA->extras,ii,jj,kk)]).Re)+
+               SQR((((CCTK_COMPLEX *)GA->data)
+		    [DATINDEX(GA->extras,ii,jj,kk)]).Im);        
+	} 
+	else 
+	{
+          CCTK_WARN (1, "PUGH_Norm2_GF: Unknown variable type");
+          return (1);
+        }
+      }
+    }
+  }
+
+  if (proc < 0) 
+  {
+    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      pughGH->PUGH_COMM_WORLD));
+  } 
+  else 
+  {
+    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      proc, pughGH->PUGH_COMM_WORLD));
+  }
+#endif
+
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    *outVal = sqrt (*outVal / tnp);
+  }
+
+  return (0);
+}
+
+
diff --git a/src/Reduction_nm1.c b/src/Reduction_nm1.c
new file mode 100644
index 0000000..db7c079
--- /dev/null
+++ b/src/Reduction_nm1.c
@@ -0,0 +1,638 @@
+ /*@@
+   @file      Reduction_nm1.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_Reduction_c)
+
+/* some useful macros */
+#ifdef MIN
+#undef MIN
+#endif
+#ifdef MAX
+#undef MAX
+#endif
+#ifdef ABS
+#undef ABS
+#endif
+#ifdef SQR
+#undef SQR
+#endif
+#define MIN(x, y)   ((x) < (y) ? (x) : (y))
+#define MAX(x, y)   ((x) > (y) ? (x) : (y))
+#define ABS(x)      ((x) < 0   ?-(x) : (x))
+#define SQR(x)      ((x) * (x))
+
+
+/* local function prototypes */
+int PUGH_Norm1_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+/*@@
+  @routine PUGH_ReduceArrayNorm1
+  @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 nOutVals 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   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+void PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
+                            int nArrays, void *inArrays [], int inType,
+                            int nOutVals, void *outVals, int outType)
+{
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayNorm1: Don't know how to reduce "
+                  "a %d-dimensional array with "
+                  "%d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type");
+    return;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+
+    /* do the reduction on input array [i] */
+    if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
+                           buffer) == 0) 
+    {
+
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+    }
+  }
+
+  free (buffer);
+}
+
+
+
+
+/*@@
+  @routine PUGH_ReduceNorm1
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the norm1.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so nOutVals must be 1.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     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     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+void PUGH_ReduceNorm1 (cGH *GH, int proc,  
+		       int nOutVals,
+                       int outType, 
+		       void *outVals, 
+		       int numinvars, 
+		       int *varlist)
+{
+  int i, elemSize;
+  pGH *pughGH;
+  CCTK_REAL intermediate_result;
+  outVals = NULL;
+
+  if (nOutVals != 1) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceNorm1 accepts only one output value per "
+                  "input array");
+    return;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0)
+  {
+      CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type");
+      return;
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+  for (i = 0; i < numinvars; i++) 
+  {
+    int result;
+
+    result = 1;
+
+    if (CCTK_GroupDimFromVarI(varlist[i])==1)
+    {
+      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
+		 "PUGH_ReduceNorm1: Reducing dimension %d not yet supported",
+		 CCTK_GroupDimFromVarI(varlist[i]));
+    }
+    else
+    {
+		 
+      switch (CCTK_GroupTypeFromVarI (varlist [i]))
+      {
+        case GROUP_GF:
+	  result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result);
+	  break;
+        case GROUP_ARRAY:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported");
+	  break;
+        case GROUP_SCALAR:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense");
+	  break;
+        default:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type");
+	  break;
+      }
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (&intermediate_result, (char *) outVals +
+                            i*nOutVals*elemSize, outType, nOutVals);
+    }
+  }
+  
+}
+
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+
+/*@@
+   @routine    PUGH_Norm1_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm1" of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+
+               "norm1" is defined as $\Sigma |a_i| / np$.
+   @enddesc 
+@@*/
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+  pGH *pughGH = PUGH_pGH (GH);
+
+
+  for (i = 0; i < nOutVals; i++) 
+  {
+    outVals [i] = 0.0;
+    for (j = 0; j < nPoints; j++) 
+    {
+      if (inType == CCTK_VARIABLE_CHAR)           /* CCTK_CHAR is unsigned */
+      {
+        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
+      }
+      else if (inType == CCTK_VARIABLE_INT)
+      {
+        outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_INT2
+      else if (inType == CCTK_VARIABLE_INT2)
+      {
+        outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT4
+      else if (inType == CCTK_VARIABLE_INT4)
+      {
+        outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT8
+      else if (inType == CCTK_VARIABLE_INT8)
+      {
+        outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_REAL)
+      {
+        outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_REAL4
+      else if (inType == CCTK_VARIABLE_REAL4)
+      {
+        outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (inType == CCTK_VARIABLE_REAL8)
+      {
+        outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (inType == CCTK_VARIABLE_REAL16)
+      {
+        outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_COMPLEX) 
+      {
+        outVals [i] += 
+          sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
+                 SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) );
+      } 
+      else 
+      {
+        CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type");
+        return (1);
+      }
+    }
+  }
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local sum */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (localOutVals);
+
+    nPoints *= pughGH->nprocs;
+  }
+
+#endif
+
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < nOutVals; i++)
+    {
+      outVals [i] /= nPoints;
+    }
+  }
+
+  return (0);
+}
+
+
+
+/*@@
+   @routine    PUGH_Norm1_GF
+   @date       Thu Apr  3 11:56:05 1997
+   @author     Paul Walker
+   @desc 
+   Returns the average of a distributed grid function.
+   @enddesc 
+@@*/
+int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
+{
+  int timelevel;
+  int ii;
+  CCTK_INT tnp;
+  pGH *pughGH;
+  pGA *GA;
+#ifdef CCTK_MPI
+  int jj, kk;
+  CCTK_INT np;
+  CCTK_REAL local_result;
+#endif
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+
+  pughGH = PUGH_pGH (GH);
+  GA = ((pGA ***) pughGH->variables) [index][timelevel];
+
+  *outVal = 0.0;
+
+#ifndef CCTK_MPI
+  for (ii = 0; ii < GA->extras->npoints; ii++)
+  {
+    if (GA->vtype == CCTK_VARIABLE_CHAR)
+    {
+      *outVal += ((CCTK_CHAR *) GA->data) [ii];     /* CCTK_CHAR is unsigned */
+    }
+    else if (GA->vtype == CCTK_VARIABLE_INT)
+    {
+      *outVal += ABS (((CCTK_INT *) GA->data) [ii]);
+    }
+#ifdef CCTK_INT2
+    else if (GA->vtype == CCTK_VARIABLE_INT2)
+    {
+      *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT4
+    else if (GA->vtype == CCTK_VARIABLE_INT4)
+    {
+      *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT8
+    else if (GA->vtype == CCTK_VARIABLE_INT8)
+    {
+      *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_REAL)
+    {
+      *outVal += ABS (((CCTK_REAL *) GA->data) [ii]);
+    }
+#ifdef CCTK_REAL4
+    else if (GA->vtype == CCTK_VARIABLE_REAL4)
+    {
+      *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL8
+    else if (GA->vtype == CCTK_VARIABLE_REAL8)
+    {
+      *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL16
+    else if (GA->vtype == CCTK_VARIABLE_REAL16)
+    {
+      *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+    {
+      *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
+                       SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) );
+    } 
+    else 
+    {
+      CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type");
+      return (1);
+    }
+  }
+
+  tnp = GA->extras->npoints;
+
+#else
+
+  /* MPI version depends heavily on the ownership array */
+  np  = 0;
+  local_result = 0.0;
+
+  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
+       kk < GA->extras->ownership[GA->stagger][1][2]; 
+       kk ++)
+  {
+    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
+         jj < GA->extras->ownership[GA->stagger][1][1]; 
+         jj ++)
+    {
+      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
+           ii < GA->extras->ownership[GA->stagger][1][0]; 
+           ii++, np++)
+      {
+        if (GA->vtype == CCTK_VARIABLE_CHAR)
+        {
+          local_result +=              /* CCTK_CHAR is unsigned */
+	    ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)];
+        }
+        else if (GA->vtype == CCTK_VARIABLE_INT)
+        {
+          local_result +=
+	    ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_INT2
+        else if (GA->vtype == CCTK_VARIABLE_INT2)
+        {
+          local_result +=
+	    ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT4
+        else if (GA->vtype == CCTK_VARIABLE_INT4)
+        {
+          local_result +=
+	    ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT8
+        else if (GA->vtype == CCTK_VARIABLE_INT8)
+        {
+          local_result +=
+	    ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_REAL)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_REAL4
+        else if (GA->vtype == CCTK_VARIABLE_REAL4)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (GA->vtype == CCTK_VARIABLE_REAL8)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (GA->vtype == CCTK_VARIABLE_REAL16)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+        {
+          local_result +=
+            sqrt( SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Re)+
+		  SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Im) );  
+        } 
+        else 
+        {
+          CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type");
+          return (1);
+        }
+      }
+    }
+  }
+
+  if (proc < 0) 
+  {
+    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      pughGH->PUGH_COMM_WORLD));
+  } 
+  else 
+  {
+    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      proc, pughGH->PUGH_COMM_WORLD));
+  }
+#endif
+ 
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    *outVal /= tnp;
+  }
+
+  return (0);
+}
+
+
diff --git a/src/SetupPGH.c b/src/SetupPGH.c
index 1e6e76e..74cf461 100644
--- a/src/SetupPGH.c
+++ b/src/SetupPGH.c
@@ -232,7 +232,8 @@ void PUGH_DestroyPGH(pGH **GHin)
 
     for (this_var = 0; this_var < pgroup.numvars; this_var++, variable++) 
     {
-      for(i = 0; i < pgroup.numtimelevels; i++)
+      /*      for(i = 0; i < pgroup.numtimelevels; i++)*/
+      for(i = 0; i < 1; i++)
       {
         switch(pgroup.grouptype)
         {
diff --git a/src/include/pughi.h b/src/include/pughi.h
index 6847bdc..5998960 100644
--- a/src/include/pughi.h
+++ b/src/include/pughi.h
@@ -147,6 +147,8 @@ int PUGH_EnableGArrayDataStorage(pGA *GA,
 int PUGH_EnableGArrayComm(pGA *GA,
                           int commflag);
 
+int PUGH_DisableGArrayDataStorage(pGA *GA);
+
 int PUGH_DisableGArrayComm(pGA *GA);
 
 int PUGH_DisableGArrayGroupComm(pGH *pughGH,
diff --git a/src/make.code.defn b/src/make.code.defn
index b512cf1..518d170 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,4 +1,4 @@
 # Main make.code.defn file for thorn pugh
 # : /usr/users/cactus/CCTK/lib/make/new_thorn.pl,v 1.1 1999/02/03 17:00:50 goodale Exp n
 # Source files in this directory
-SRCS =  Overloadables.c Evolve.c Startup.c GHExtension.c Comm.c Storage.c SetupPGH.c SetupGroup.c  Reduction.c PughUtils.c SetupPGV.c PostSendGA.c PostReceiveGA.c FinishReceiveGA.c LoadAware.c
+SRCS =  ReductionNorm1.c ReductionNorm2.c ReductionMin.c ReductionMax.c Overloadables.c Evolve.c Startup.c GHExtension.c Comm.c Storage.c SetupPGH.c SetupGroup.c  Reduction.c PughUtils.c SetupPGV.c PostSendGA.c PostReceiveGA.c FinishReceiveGA.c LoadAware.c
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
index 4ae823c..e677991 100644
--- a/src/pugh_reductions.h
+++ b/src/pugh_reductions.h
@@ -18,16 +18,16 @@
 extern "C" {
 #endif
 
-void PUGH_ReduceMinVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceMaxVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceNorm1(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceNorm2(REDUCTION_OPERATOR_REGISTER_ARGLIST);
-
-void PUGH_ReduceArrayMinVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceArrayMaxVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceArraySum(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceArrayNorm1(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
-void PUGH_ReduceArrayNorm2(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceMinVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceMaxVal(REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceNorm1(REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceNorm2(REDUCTION_OPERATOR_REGISTER_ARGLIST);
+
+int PUGH_ReduceArrayMinVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceArrayMaxVal(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceArraySum(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceArrayNorm1(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReduceArrayNorm2(REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
 
 #ifdef _cplusplus
 }