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

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

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

14 files changed, 2550 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..057cc49
--- /dev/null
+++ b/README
@@ -0,0 +1,8 @@
+Cactus Code Thorn PUGHReduce
+Authors    : Cactus Team
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Standard reduction operators which use PUGH for communications
diff --git a/doc/documention.tex b/doc/documention.tex
new file mode 100644
index 0000000..68f584b
--- /dev/null
+++ b/doc/documention.tex
@@ -0,0 +1,118 @@
+\documentclass{article}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% MANPAGE like description setting for options, use as 
+% \begin{Lentry} \item[text] text  \end{Lentry}
+
+\usepackage{ifthen,calc}
+
+\newcommand{\entrylabel}[1]{\mbox{\textsf{#1}}\hfil}
+\newenvironment{entry}
+  {\begin{list}{}
+    {\renewcommand{\makelabel}{\entrylabel}
+      \setlength{\labelwidth}{90pt}
+      \setlength{\leftmargin}{\labelwidth+\labelsep}
+    }
+  }
+  {\end{list}}
+
+\newlength{\Mylen}
+\newcommand{\Lentrylabel}[1]{%
+  \settowidth{\Mylen}{\textsf{#1}}%
+  \ifthenelse{\lengthtest{\Mylen > \labelwidth}}%
+    {\parbox[b]{\labelwidth} %  term > labelwidth
+      {\makebox[0pt][l]{\textsf{#1}}\\}} %
+    {\textsf{#1}} %
+
+  \hfil\relax}
+\newenvironment{Lentry}
+  {\renewcommand{\entrylabel}{\Lentrylabel}
+   \begin{entry}}
+  {\end{entry}}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+\begin{document}
+
+\title{PUGHReduce}
+\author{Gabrielle Allen, Thomas Radke}
+\date{2000}
+\maketitle
+
+\abstract{Reductions operations which are performed using the PUGH driver}
+
+\section{Purpose}
+
+This thorn registers a number of reduction operators with the flesh. The 
+reductions are performed using internals of the PUGH driver, so that this
+thorn can only be used when {\tt CactusPUGH/PUGH} is active.
+
+The reduction operations this thorn registers are
+
+\begin{Lentry}
+\item[maximum]
+	Calculates the maximum of a grid variable
+\item[minimum]
+	Calculates the minimum of a grid variable
+\item[sum]
+	Calculates the sum of the elements of a grid variable
+\item[norm1]
+	Calculates the norm of a grid variable
+	$$
+	\left(\Sigma | GV | \right)/n
+	$$
+	where $n$ is the number of elements
+	
+\item[norm2]
+	Calculates the norm of a grid variable
+	$$
+	\sqrt{\Sigma GV^2}/n
+	$$
+	where $n$ is the number of elements
+	
+\end{Lentry}
+
+\section{Examples}
+
+\begin{itemize}
+
+\item{} Maximum Value of a Grid Function
+
+{\tt
+\begin{verbatim}
+    int ierr;
+    int handle;         /* handle for reduction operator */
+    int index; 	        /* grid variable index */
+    CCTK_REAL *value;   /* value returned from reduction */
+    const char *reduction_name = "maximum";
+                        /* reduction operation to use */
+
+    handle = CCTK_ReductionHandle (reduction_name);
+    if (handle >= 0) 
+    {
+      if ( CCTK_Reduce (cctkGH, 0, handle, 
+                        1, CCTK_VARIABLE_REAL,
+                        &value, 1, index) == 0 
+           && CCTK_MyProc(cctkGH) == 0)
+      {
+        printf("Reduction (%s) value is %f\n",reduction_name,value);
+      }
+    }
+    else
+    {
+      CCTK_WARN(1,"Invalid reduction operator"); 
+    }
+\end{verbatim}
+}
+
+\end{itemize}
+
+
+% Automatically created from the ccl files by using gmake thorndoc
+\include{interface}
+\include{param}
+\include{schedule}
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..4e5d4a6
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn PUGHReduce
+# $Header$
+
+implements: reduce
+inherits: driver
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..0c18bee
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,3 @@
+# Parameter definitions for thorn PUGHReduce
+# $Header$
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..eaed53c
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,8 @@
+# Schedule definitions for thorn PUGHReduce
+# $Header$
+
+schedule PUGHReduce_Startup at STARTUP
+{
+  LANG:C
+} "Startup routine."
+
diff --git a/src/Reduction.c b/src/Reduction.c
new file mode 100644
index 0000000..80266d8
--- /dev/null
+++ b/src/Reduction.c
@@ -0,0 +1,461 @@
+ /*@@
+   @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 <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_Reduction_c)
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+static int PUGH_ReductionGA (cGH *GH, 
+			     int index, 
+			     int proc, 
+			     CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn);
+
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray [/* num_elems */],
+                                 int outtype,
+                                 void *outarray     /* [num_elems] */);
+
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    PUGH_ReductionArrays
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+int PUGH_ReductionArrays (cGH *GH,
+                          int proc,
+                          int num_dims,
+                          int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          void *inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn)
+{
+  int retval;
+  int i;
+  int num_points;
+  int from[1];
+  int to[1];
+  int iterator[1];
+  int points_per_dim[1];
+  int *intypes;
+  CCTK_REAL *buffer;
+
+  points_per_dim[0] = 1;
+  from[0] = 0;
+  to[0] = dims[0];
+
+  for (i = 1; i < num_dims; i++)
+  {
+    to[0] *= dims[i];
+  }
+
+  if (num_outvals != 1)
+  {
+    if (num_outvals != to[0])
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReductionArrays: Don't know how to reduce "
+                  "a %d-dimensional array with %d elements "
+                  "to an output array of %d elements",
+                  num_dims, to[0], num_outvals);
+      return (-1);
+    }
+    to[0] = 1;
+  }
+  num_points = to[0] * CCTK_nProcs (GH);
+
+  /* set the array types to intype */
+  /* FIXME: could allow to pass in arrays of different types now !!! */
+  intypes = (int *) malloc (num_inarrays * sizeof (int));
+  for (i = 0; i < num_inarrays; i++)
+  {
+    intypes[i] = intype;
+  }
+
+  buffer = (CCTK_REAL *) malloc (num_outvals * sizeof (CCTK_REAL));
+
+  /* do the reduction on the input arrays */
+  retval = reduction_fn (GH, 
+			 proc, 
+			 num_dims, 
+			 from, 
+			 to, 
+			 iterator, 
+			 points_per_dim,
+                         num_points, 
+			 num_inarrays, 
+			 intypes, 
+			 inarrays,
+                         num_outvals, 
+			 buffer);
+
+  if (retval == 0 && (proc < 0 || proc == CCTK_MyProc (GH)))
+  {
+    /* type-cast the result to the requested datatype */
+    retval = copy_real_to_outtype (num_inarrays * num_outvals, 
+				   buffer,
+                                   outtype, 
+				   outvals);
+  }
+
+  free (intypes);
+  free (buffer);
+
+  return (retval);
+}
+
+
+ /*@@
+   @routine    PUGH_ReductionArrays
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+int PUGH_ReductionGVs (cGH *GH,
+                       int proc,
+                       int num_invars,
+                       int invars[/* num_invars */],
+                       int outtype,
+                       int num_outvals,
+                       void *outvals /* [num_outvals] */,
+                       reduction_fn_t reduction_fn)
+{
+  int i, myproc, outtypesize, this_retval, retval;
+  CCTK_REAL result;
+
+
+  if (num_outvals != 1)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Only one output value per "
+                  "input array allowed");
+    return (-1);
+  }
+
+  outtypesize = CCTK_VarTypeSize (outtype);
+  if (outtypesize <= 0)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Invalid output variable type");
+    return (-1);
+  }
+
+  retval = 0;
+  myproc = CCTK_MyProc (GH);
+
+  for (i = 0; i < num_invars; i++)
+  {
+    switch (CCTK_GroupTypeFromVarI (invars[i]))
+    {
+      case CCTK_GF:
+      case CCTK_ARRAY:
+        this_retval = PUGH_ReductionGA (GH, invars[i], proc, &result,
+                                        reduction_fn);
+        break;
+
+      case CCTK_SCALAR:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Reducing scalars doesn't make sense");
+        this_retval = -1;
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Invalid variable index");
+        this_retval = -1;
+        break;
+    }
+
+    if (this_retval == 0 && (proc < 0 || proc == myproc))
+    {
+      /* type-cast the result to the requested datatype */
+      this_retval = copy_real_to_outtype (1, &result, outtype,
+                                          (char *) outvals + i*outtypesize);
+    }
+
+    retval |= this_retval;
+  }
+
+  return (retval);
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+
+ /*@@
+   @routine    PUGH_ReductionGA
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+static int PUGH_ReductionGA (cGH *GH, int index, int proc, CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn)
+{
+  int i, retval, stagger_index, num_points;
+  pGA *GA;
+  int *from, *to, *iterator, *points_per_dim;
+
+
+  /* get the GA structure for the current timelevel */
+  i = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (i > 0)
+  {
+    i--;
+  }
+  GA = ((pGA ***) PUGH_pGH (GH)->variables)[index][i];
+
+  /* allocate temporary vectors */
+  from           = (int *) malloc (4 * GA->connectivity->dim * sizeof (int));
+  to             = from + 1*GA->connectivity->dim;
+  iterator       = from + 2*GA->connectivity->dim;
+  points_per_dim = from + 3*GA->connectivity->dim;
+
+  /* get the start and endpoint to iterate over, the total number of points,
+     and the points_per_dim[] vector */
+  num_points = 1;
+  points_per_dim[0] = 1;
+  for (i = 0; i < GA->connectivity->dim; i++)
+  {
+    stagger_index = CCTK_StaggerDirIndex (i, GA->stagger);
+
+    from[i] = GA->extras->ownership[stagger_index][0][i];
+    to[i]   = GA->extras->ownership[stagger_index][1][i];
+
+    if (stagger_index == 0)
+    {
+      num_points *= GA->extras->nsize[i];
+    }
+    else
+    {
+      num_points *= GA->extras->nsize[i] - 1;
+    }
+
+    if (i > 0)
+    {
+      points_per_dim[i] = points_per_dim[i-1] * GA->extras->lnsize[i-1];
+    }
+  }
+
+  /* now do the reduction */
+  retval = reduction_fn (GH, proc, GA->connectivity->dim, from, to, iterator,
+                         points_per_dim, num_points, 1, &GA->vtype, &GA->data,
+                         1, outval);
+
+  /* free temporary vectors */
+  free (from);
+
+  return (retval);
+}
+
+
+ /*@@
+   @routine    copy_real_to_outtype
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray[/* num_elems */],
+                                 int outtype,
+                                 void *outarray /* [num_elems] */)
+{
+  int i, retval;
+
+
+  retval = 0;
+
+  if (outtype == CCTK_VARIABLE_CHAR)
+  {
+    CCTK_CHAR *_outarray = (CCTK_CHAR *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_CHAR) inarray[i];
+    }
+  }
+  else if (outtype == CCTK_VARIABLE_INT)
+  {
+    CCTK_INT *_outarray = (CCTK_INT *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT) inarray[i];
+    }
+  }
+#ifdef CCTK_INT2
+  else if (outtype == CCTK_VARIABLE_INT2)
+  {
+    CCTK_INT2 *_outarray = (CCTK_INT2 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT2) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_INT4
+  else if (outtype == CCTK_VARIABLE_INT4)
+  {
+    CCTK_INT4 *_outarray = (CCTK_INT4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT4) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_INT8
+  else if (outtype == CCTK_VARIABLE_INT8)
+  {
+    CCTK_INT8 *_outarray = (CCTK_INT8 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT8) inarray[i];
+    }
+  }
+#endif
+  else if (outtype == CCTK_VARIABLE_REAL)
+  {
+    CCTK_REAL *_outarray = (CCTK_REAL *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL) inarray[i];
+    }
+  }
+#ifdef CCTK_REAL4
+  else if (outtype == CCTK_VARIABLE_REAL4)
+  {
+    CCTK_REAL4 *_outarray = (CCTK_REAL4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL4) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_REAL8
+  else if (outtype == CCTK_VARIABLE_REAL8)
+  {
+	CCTK_REAL8 *_outarray = (CCTK_REAL8 *) outarray;
+
+
+	for (i = 0; i < num_elems; i++)
+	{
+	  _outarray[i] = (CCTK_REAL8) inarray[i];
+	}
+  }
+#endif
+#ifdef CCTK_REAL16
+  else if (outtype == CCTK_VARIABLE_REAL16)
+  {
+    CCTK_REAL16 *_outarray = (CCTK_REAL16 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL16) inarray[i];
+    }
+  }
+#endif
+  else
+  {
+    CCTK_WARN (1, "copy_real_to_outtype: Unsupported output type");
+    retval = -1;
+  }
+
+  return (retval);
+}
diff --git a/src/ReductionMax.c b/src/ReductionMax.c
new file mode 100644
index 0000000..7d2c53d
--- /dev/null
+++ b/src/ReductionMax.c
@@ -0,0 +1,331 @@
+ /*@@
+   @file      ReductionMax.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+              Defines the PUGH reduction operator to get the maximum
+              for a set of arbitrary arrays.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReductionMax_c)
+
+/* local function prototypes */
+static int PUGH_ReductionMaxVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionMaxValArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the maxima
+           for a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     num_inarrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     intype
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionMaxValArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionMaxVal));
+}
+
+
+/*@@
+  @routine PUGH_ReductionMaxValGVs
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the maximum.
+           It just calls the general routine PUGH_ReductionGVs() passing it
+           the callback to compute the maxima for a set of arrays.
+           If the number of requested output values equals the size of the
+           variable, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionMaxValGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionMaxVal));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionMaxVal
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the maximum of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+   @enddesc
+@@*/
+static int PUGH_ReductionMaxVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+#ifdef CCTK_MPI
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     ((array)[0])
+#define REDUCTION_OPERATION(max, scalar)   if (max < scalar) max = scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_CHAR, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+      break;
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+        CCTK_WARN (1, "PUGH_ReductionMaxVal: Don't know how to compute "
+                      "the maximum of complex variables !!!");
+        return (-1);
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionMaxVal: Unknown variable type");
+        return (-1);
+    }
+  }
+
+#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
+
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local maximum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_MAX, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_MAX, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  return (0);
+}
diff --git a/src/ReductionMin.c b/src/ReductionMin.c
new file mode 100644
index 0000000..5e08255
--- /dev/null
+++ b/src/ReductionMin.c
@@ -0,0 +1,331 @@
+ /*@@
+   @file      ReductionMin.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+              Defines the PUGH reduction operator to get the minimum
+              for a set of arbitrary arrays.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReductionMin_c)
+
+/* local function prototypes */
+static int PUGH_ReductionMinVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionMinValArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the minima
+           for a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     num_inarrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     intype
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionMinValArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionMinVal));
+}
+
+
+/*@@
+  @routine PUGH_ReductionMinValGVs
+  @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 num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionMinValGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionMinVal));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionMinVal
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the minimum for a set of distributed array with
+               'num_points' elements each.
+               Global reduction is done element-wise (num_outvals == 1)
+               or on the results of the local reductions.
+   @enddesc
+@@*/
+static int PUGH_ReductionMinVal (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+#ifdef CCTK_MPI
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     ((array)[0])
+#define REDUCTION_OPERATION(min, scalar)   if (min > scalar) min = scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_CHAR, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+        CCTK_WARN (1, "PUGH_ReductionMinVal: Don't know how to compute "
+                      "the minimum of complex variables !!!");
+        return (-1);
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionMinVal: Unknown variable type");
+        return (-1);
+    }
+  }
+
+#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
+
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local minimum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_MIN, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_MIN, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+#endif
+
+  return (0);
+}
diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
new file mode 100644
index 0000000..efa2f0b
--- /dev/null
+++ b/src/ReductionNorm1.c
@@ -0,0 +1,354 @@
+ /*@@
+   @file      ReductionNorm1.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm1 = $\Sigma |a_i| / np$
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReductionNorm1_c)
+
+/* local function prototypes */
+static int PUGH_ReductionNorm1 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm1Arrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm1"
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           If the number of requested output values equals the number of
+           elements in an input array, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     intype
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionNorm1Arrays (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int dims[/* num_dims */],
+                               int num_inarrays,
+                               void *inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionNorm1));
+}
+
+
+/*@@
+  @routine PUGH_ReductionNorm1GVs
+  @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 num_outvals 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     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionNorm1GVs (cGH *GH,
+                            int proc,  
+                            int num_outvals,
+                            int outtype, 
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionNorm1));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionNorm1
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm1" of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+
+               "norm1" is defined as $\Sigma |a_i| / np$.
+   @enddesc 
+@@*/
+static int PUGH_ReductionNorm1 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+  pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(data)       0
+#ifdef ABS
+#undef ABS
+#endif
+#define ABS(x)      ((x) < 0 ? -(x) : (x))
+
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += scalar
+        ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_CHAR, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += ABS (scalar)
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm1, scalar)   norm1 += CCTK_CmplxAbs (scalar)
+        ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm1: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1) 
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  /* finally assign the return value */
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < total_outvals; i++)
+    {
+      outvals[i] /= num_points;
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionNorm2.c b/src/ReductionNorm2.c
new file mode 100644
index 0000000..beb2cad
--- /dev/null
+++ b/src/ReductionNorm2.c
@@ -0,0 +1,354 @@
+/*@@
+   @file      ReductionNorm2.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc
+              Defines the PUGH reduction operator to get the norm
+              of an arbitrary array which is defined as:
+                  norm2 = $\sqrt{\Sigma (a_i * a_i) / np}$
+   @enddesc
+   @version $Header$
+ @@*/
+
+#include <math.h>
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReductionNorm2_c)
+
+/* local function prototypes */
+static int PUGH_ReductionNorm2 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionNorm2Arrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm2"
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           If the number of requested output values equals the number of
+           elements in an input array, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_dims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     num_arrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inarrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     intype
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionNorm2Arrays (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int dims[/* num_dims */],
+                               int num_inarrays,
+                               void *inarrays[/* num_inarrays */],
+                               int intype,
+                               int num_outvals,
+                               void *outvals /* [num_outvals] */,
+                               int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionNorm2));
+}
+
+
+/*@@
+  @routine PUGH_ReductionNorm2GVs
+  @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 num_outvals 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     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionNorm2GVs (cGH *GH,
+                            int proc,  
+                            int num_outvals,
+                            int outtype, 
+                            void *outvals /* [num_outvals] */,
+                            int num_invars,
+                            int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars,
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionNorm2));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionNorm2
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm2" of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+
+               "norm2" is defined as $\sqrt{\Sigma (a_i * a_i) / np}$.
+   @enddesc 
+@@*/
+static int PUGH_ReductionNorm2 (cGH *GH,
+                                int proc,
+                                int num_dims,
+                                int from[/* dim */],
+                                int to[/* dim */],
+                                int iterator[/* dim */],
+                                int points_per_dim[/* dim */],
+                                int num_points,
+                                int num_inarrays,
+                                int intypes[/* num_inarrays */],
+                                void *inarrays[/* num_inarrays */],
+                                int num_outvals,
+                                CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+  pGH *pughGH;
+#ifdef CCTK_MPI
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)       0
+#ifdef SQR
+#undef SQR
+#endif
+#define SQR(x)      ((x) * (x))
+
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+#define REDUCTION_OPERATION(norm2, scalar)   norm2 += SQR (scalar)
+        ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_CHAR, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+#undef  REDUCTION_OPERATION
+#define REDUCTION_OPERATION(norm2, scalar)   norm2 += SQR ((scalar).Re) +     \
+                                                      SQR ((scalar).Im)
+        ITERATE_ARRAY (CCTK_COMPLEX, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionNorm2: Unknown variable type");
+        return (-1);
+    }
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+#ifdef CCTK_MPI
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1) 
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  /* finally assign the return value */
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < total_outvals; i++)
+    {
+      outvals[i] = sqrt (outvals[i] / num_points);
+    }
+  }
+
+  return (0);
+}
diff --git a/src/ReductionSum.c b/src/ReductionSum.c
new file mode 100644
index 0000000..e4bde2b
--- /dev/null
+++ b/src/ReductionSum.c
@@ -0,0 +1,331 @@
+ /*@@
+   @file      ReductionSum.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+              Defines the PUGH reduction operator to get the sum
+              of an arbitrary array.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReductionSum_c)
+
+/* local function prototypes */
+static int PUGH_ReductionSum (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+
+/*@@
+  @routine PUGH_ReductionSumArrays
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the sums
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           For the number of output values only 1 is accepted.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionSumArrays (cGH *GH, 
+                                int proc, 
+                                int num_dims, 
+                                int dims[/* num_dims */],
+                                int num_inarrays, 
+                                void *inarrays[/* num_inarrays */], 
+                                int intype,
+                                int num_outvals, 
+                                void *outvals /* [num_outvals] */,
+                                int outtype)
+{
+  return (PUGH_ReductionArrays (GH, proc, num_dims, dims,
+                                intype, num_inarrays, inarrays,
+                                outtype, num_outvals, outvals,
+                                PUGH_ReductionSum));
+}
+
+
+/*@@
+  @routine PUGH_ReductionSumGVs
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the sum.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           If the number of requested output values equals the size of the
+           variable, global reduction is done element-wise.
+           Otherwise num_outvals must be 1, and global reduction is done on the
+           results of the local reductions.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     num_outvals
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outtype
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outvals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     num_invars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     invars
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReductionSumGVs (cGH *GH, 
+                             int proc, 
+                             int num_outvals,
+                             int outtype, 
+                             void *outvals /* [num_outvals] */,
+                             int num_invars, 
+                             int invars[/* num_invars */])
+{
+  return (PUGH_ReductionGVs (GH, proc, num_invars, invars, 
+                             outtype, num_outvals, outvals,
+                             PUGH_ReductionSum));
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    PUGH_ReductionSum
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the sum of a distributed array with
+               'num_points' elements. Global reduction is done element-wise
+               (num_outvals == 1) or on the results of the local reductions.
+   @enddesc
+@@*/
+static int PUGH_ReductionSum (cGH *GH,
+                                 int proc,
+                                 int num_dims,
+                                 int from[/* dim */],
+                                 int to[/* dim */],
+                                 int iterator[/* dim */],
+                                 int points_per_dim[/* dim */],
+                                 int num_points,
+                                 int num_inarrays,
+                                 int intypes[/* num_inarrays */],
+                                 void *inarrays[/* num_inarrays */],
+                                 int num_outvals,
+                                 CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+{
+  int i, total_outvals;
+#ifdef CCTK_MPI
+  pGH *pughGH;
+  CCTK_REAL *local_outvals;
+#endif
+
+
+/* macros to complete the ITERATE_ARRAY macro */
+#define INITIAL_REDUCTION_VALUE(array)     0
+#define REDUCTION_OPERATION(sum, scalar)   sum += scalar
+
+  for (i = total_outvals = 0; i < num_inarrays; i++)
+  {
+    switch (intypes[i])
+    {
+      case CCTK_VARIABLE_CHAR:
+        ITERATE_ARRAY (CCTK_CHAR, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_CHAR, outvals, num_outvals, total_outvals);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        ITERATE_ARRAY (CCTK_INT, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        ITERATE_ARRAY (CCTK_INT2, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT2, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        ITERATE_ARRAY (CCTK_INT4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        ITERATE_ARRAY (CCTK_INT8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_INT8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_REAL:
+        ITERATE_ARRAY (CCTK_REAL, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL, outvals, num_outvals, total_outvals);
+        break;
+
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        ITERATE_ARRAY (CCTK_REAL4, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL4, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        ITERATE_ARRAY (CCTK_REAL8, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL8, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        ITERATE_ARRAY (CCTK_REAL16, num_dims, inarrays[i],
+                       from, to, iterator, points_per_dim,
+                       CCTK_REAL16, outvals, num_outvals, total_outvals);
+        break;
+#endif
+
+      case CCTK_VARIABLE_COMPLEX:
+        CCTK_WARN (1, "PUGH_ReductionSum: Don't know how to compute "
+                      "the sum of complex variables !!!");
+        return (-1);
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionSum: Unknown variable type");
+        return (-1);
+    }
+  }
+
+#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
+
+  /* do the global reduction from the processors' local reduction results */
+  if (pughGH->nprocs > 1)
+  {
+    local_outvals = (CCTK_REAL *) malloc (total_outvals * sizeof (CCTK_REAL));
+
+    /* outvals[] contains now the local sum */
+    memcpy (local_outvals, outvals, total_outvals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (local_outvals, outvals, total_outvals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (local_outvals);
+  }
+
+#endif
+
+  return (0);
+}
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..56040da
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,51 @@
+ /*@@
+   @file      Startup.c
+   @date      Wed Feb  3 23:10:19 1999
+   @author    Tom Goodale
+   @desc 
+   Startup routines for PUGHReduce.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include "cctk.h"
+
+#include "pugh_reductions.h"
+
+static char *rcsid="$Header$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_Startup_c)
+
+
+ /*@@
+   @routine    PUGHReduce_Startup
+   @date       Wed Feb  3 23:14:38 1999
+   @author     Tom Goodale
+   @desc 
+   The startup registration routine for PUGHReduce.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+ 
+@@*/
+int PUGHReduce_Startup(void)
+{
+
+  /* Register the reduction operators provided by PUGH */
+  CCTK_RegisterReductionOperator(PUGH_ReductionMinValGVs, "minimum");
+  CCTK_RegisterReductionOperator(PUGH_ReductionMaxValGVs, "maximum");
+  CCTK_RegisterReductionOperator(PUGH_ReductionSumGVs,    "sum");
+  CCTK_RegisterReductionOperator(PUGH_ReductionNorm1GVs,  "norm1");
+  CCTK_RegisterReductionOperator(PUGH_ReductionNorm2GVs,  "norm2");
+
+  CCTK_RegisterReductionArrayOperator(PUGH_ReductionMinValArrays, "minimum");
+  CCTK_RegisterReductionArrayOperator(PUGH_ReductionMaxValArrays, "maximum");
+  CCTK_RegisterReductionArrayOperator(PUGH_ReductionSumArrays,    "sum");
+  CCTK_RegisterReductionArrayOperator(PUGH_ReductionNorm1Arrays,  "norm1");
+  CCTK_RegisterReductionArrayOperator(PUGH_ReductionNorm2Arrays,  "norm2");
+
+  return 0;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..0389a78
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn PUGHReduce
+# $Header$
+
+# Source files in this directory
+SRCS = ReductionMax.c ReductionMin.c ReductionNorm1.c ReductionNorm2.c ReductionSum.c Startup.c Reduction.c
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
new file mode 100644
index 0000000..a1a3c52
--- /dev/null
+++ b/src/pugh_reductions.h
@@ -0,0 +1,186 @@
+ /*@@
+   @header    pugh_reductions.h
+   @date      April 29 1999
+   @author    Gabrielle
+   @desc 
+   Prototypes for pugh reduction operators
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#ifndef _PUGH_REDUCTIONS_H_
+#define _PUGH_REDUCTIONS_H_
+
+
+#include "cctk.h"
+#include "cctk_Reduction.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+
+/***
+   Macro to iterate over every element of an arbitrary sized array
+
+   The array is given by
+     - its CCTK type, cctk_type,
+     - the number of dimensions, vdim,
+     - an (untyped) pointer to its data, vdata,
+     - one-dimensional vectors, from[] and to[], of length vdim,
+       specifying the start- and endpoint to iterate over for each dimension
+
+   The result of the reduction is specified by
+     - the array where to store the results, outvals[]
+     - the type of the intermediate result scalar, outvals_type
+     - the number of scalars to reduce from this input array, num_outvals
+     - the total number of scalars for all input arrays, total_outvals
+       (incremented after each iteration)
+
+   An iteration vector, iterator[], is passed in to hold the current indices
+   for each dimension of the array.
+   Another vector of size vdim, points_per_dim[], is used to compute the
+   linear index out of the iterator[] elements.
+
+   Before the iteration loop, a local reduction variable, typed_outval,
+   is defined and initialized by the macro INITIAL_REDUCTION_VALUE(typed_vdata).
+
+   For each iteration the macro REDUCTION_OPERATION(typed_outval, typed_inval)
+   is called, with the local reduction variable and the pointer to the
+   current data element.
+
+   Finally, the local reduction variable, typed_outval, is type-casted
+   to CCTK_REAL and assigned to the current output array value outvals[].
+   The total output scalar counter, total_outvals, is incremented.
+
+
+   NOTE: The macro assumes from[i] > to[i] for all 1 <= i < vdim such that
+         at least one iteration will be performed.
+         Make sure this is true when you set up the from[] and to[] vectors.
+
+         For efficiency reasons, the innermost loop (the one for dimension 0)
+         is coded as a for() loop. Thus the overhead to iterate over a
+         one-dimensional array should be kept to a minimum.
+***/
+#define ITERATE_ARRAY(cctk_type, vdim, vdata, from, to, iterator,             \
+                      points_per_dim,                                         \
+                      outvals_type, outvals, num_outvals, total_outvals)      \
+    {                                                                         \
+      int _i, _j, _dim, _vindex;                                              \
+      cctk_type *typed_vdata = (cctk_type *) (vdata);                         \
+      outvals_type typed_outval;                                              \
+                                                                              \
+                                                                              \
+      for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                     \
+      {                                                                       \
+        /* get the linear index of the element to start with */               \
+        _vindex = from[0];                                                    \
+        for (_i = 1; _i < vdim; _i++)                                         \
+          _vindex += from [_i] * points_per_dim [_i];                         \
+        typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);        \
+                                                                              \
+        /* set iterator to local startpoint */                                \
+        memcpy (iterator, from, vdim * sizeof (int));                         \
+                                                                              \
+        /* do the nested loops starting with the second-innermost */          \
+        _dim = 1;                                                             \
+        while (1)                                                             \
+        {                                                                     \
+          /* get the linear index */                                          \
+          _vindex = 0;                                                        \
+          for (_i = 1; _i < vdim; _i++)                                       \
+            _vindex += iterator [_i] * points_per_dim [_i];                   \
+                                                                              \
+          /* do the reduction for the innermost loop (lowest dimension) */    \
+          for (_i = from [0]; _i < to [0]; _i++)                              \
+          {                                                                   \
+            REDUCTION_OPERATION (typed_outval, typed_vdata [_i + _vindex]);   \
+          }                                                                   \
+                                                                              \
+          if (vdim > 1)                                                       \
+          {                                                                   \
+            /* increment current looper and check for end */                  \
+            if (++iterator [_dim] >= to [_dim])                               \
+            {                                                                 \
+              /* increment outermost loopers */                               \
+              for (_dim++; _dim < vdim; _dim++)                               \
+              {                                                               \
+                if (++iterator [_dim] < to [_dim])                            \
+                  break;                                                      \
+              }                                                               \
+                                                                              \
+              /* done if beyond outermost loop */                             \
+              if (_dim >= vdim)                                               \
+                break;                                                        \
+                                                                              \
+              /* reset innermost loopers */                                   \
+              for (_dim--; _dim >= 0; _dim--)                                 \
+                iterator [_dim] = from [_dim];                                \
+              _dim = 1;                                                       \
+            }                                                                 \
+          }                                                                   \
+          else                                                                \
+          {                                                                   \
+            /* exit loop if array is one-dimensional */                       \
+            break;                                                            \
+          }                                                                   \
+                                                                              \
+        } /* end of nested loops over all dimensions */                       \
+                                                                              \
+        outvals [total_outvals++] = (CCTK_REAL) typed_outval;                 \
+                                                                              \
+      } /* end of loop over num_outvals */                                    \
+    }
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int PUGH_ReductionMinValGVs (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValGVs (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumGVs    (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm1GVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm2GVs  (REDUCTION_OPERATOR_REGISTER_ARGLIST);
+
+int PUGH_ReductionMinValArrays (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionMaxValArrays (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionSumArrays    (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm1Arrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int PUGH_ReductionNorm2Arrays  (REDUCTION_ARRAY_OPERATOR_REGISTER_ARGLIST);
+
+typedef int (*reduction_fn_t) (cGH *GH,
+                               int proc,
+                               int num_dims,
+                               int from[/* dim */],
+                               int to[/* dim */],
+                               int iterator[/* dim */],
+                               int points_per_dim[/* dim */],
+                               int num_points,
+                               int num_inarrays,
+                               int intypes[/* num_inarrays */],
+                               void *inarrays[/* num_inarrays */],
+                               int num_outvals,
+                               CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+
+int PUGH_ReductionGVs (cGH *GH,
+                       int proc,
+                       int num_invars,
+                       int invars[/* num_invars */],
+                       int outtype,
+                       int num_outvals,
+                       void *outvals /* [num_outvals] */,
+                       reduction_fn_t reduction_fn);
+int PUGH_ReductionArrays (cGH *GH,
+                          int proc,
+                          int num_dims,
+                          int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          void *inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif