New global reduction stuff.

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

4 files changed, 563 insertions, 4 deletions

diff --git a/src/ReduceGA.c b/src/ReduceGA.c
new file mode 100644
index 0000000..5e4e9ac
--- /dev/null
+++ b/src/ReduceGA.c
@@ -0,0 +1,555 @@
+ /*@@
+   @file      ReduceGA.c
+   @date      Tue Aug 31 13:29:04 2004
+   @author    Tom Goodale, Yaakoub Y El Khamra
+   @desc
+              New grid arrays implementation
+   @enddesc
+   @version   $Header$
+ @@*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "util_Table.h"
+#include "pugh_reductions.h"
+
+static const char *rcsid = "$Id:";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReduceGA_c);
+
+/* local function prototypes */
+static int ReduceGridArrays (const cGH *GH,
+                          int local_reduce_handle,
+                          int param_table_handle,
+                          int N_input_arrays,
+                          const CCTK_INT input_array_variable_indices[],
+                          int M_output_values,
+                          const CCTK_INT output_value_type_codes[],
+                          void* const output_values[]);
+
+
+ /*@@
+   @routine    PUGH_ReduceGlobalArrays
+   @date       
+   @author     Gabrielle Allen, Yaakoub El Khamra
+   @desc
+               Generic routine for doing a reduction operation on a set of
+               Cactus variables.
+   @enddesc
+   @var     GH
+   @vdesc   pointer to the grid hierarchy
+   @vtype   cGH *
+   @vio     in
+   @endvar
+   @var     local_reduce_handle
+   @vdesc   the handle specifying the reduction operator
+   @vtype   int
+   @vio     in
+   @endvar
+   @var     param_reduce_handle
+   @vdesc   the parameter table handle
+   @vtype   int
+   @vio     in
+   @endvar
+   @var     N_input_arrays
+   @vdesc   number of elements in the reduction input
+   @vtype   int
+   @vio     in
+   @endvar
+   @var     input_array_variable_indices
+   @vdesc   input arrays
+   @vtype   const CCTK_INT
+   @vio     in
+   @endvar
+   @var     M_output_values
+   @vdesc   number of output values
+   @vtype   int
+   @vio     in
+   @endvar
+   @var     output_value_type_codes
+   @vdesc   array containing output value type codes
+   @vtype   const CCTK_IN
+   @vio     in
+   @endvar
+   @var     output_values
+   @vdesc   array of output values
+   @vtype   void * const
+   @vio     in
+   @endvar
+@@*/
+int PUGH_ReduceGridArrays (const cGH *GH,
+                           int local_reduce_handle,
+                           int param_table_handle,
+                           int N_input_arrays,
+                           const CCTK_INT input_array_variable_indices[],
+                           int M_output_values,
+                           const CCTK_INT output_value_type_codes[],
+                           void* const output_values[])
+{
+  int retval;
+  retval = ReduceGridArrays(GH, local_reduce_handle,
+                          param_table_handle, N_input_arrays,
+                          input_array_variable_indices,
+                          M_output_values, output_value_type_codes, 
+                          output_values);
+  return retval;
+}
+
+
+/*****************************************************************************/
+/*                             local functions                               */
+/*****************************************************************************/
+/*@@
+   @routine    ReduceGridArrays
+   @date       
+   @author     Tom Goodale, Yaakoub Y El Khamra
+   @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 ReduceGridArrays (const cGH *GH,
+                          int local_reduce_handle,
+                          int param_table_handle,
+                          int N_input_arrays,
+                          const CCTK_INT input_array_variable_indices[],
+                          int M_output_values,
+                          const CCTK_INT output_value_type_codes[],
+                          void* const output_values[])
+{
+  /* Utility variables */
+  int i, ierr = 0;
+  int dim = 0;
+  int proc;
+  int num_points=0;
+  int total_num_points = 0;
+  int perform_division = 1;
+
+  CCTK_INT * lower_array_bounds;
+  CCTK_INT * upper_array_bounds;
+  CCTK_INT * input_array_dims;
+  CCTK_POINTER * input_arrays;
+  CCTK_INT * input_array_type_codes;
+
+
+#ifdef CCTK_MPI
+  int nprocs = 0, myproc =0, mpi_operation = 0;
+  const pGH *pughGH = NULL;
+  void *local_outvals = NULL;
+  nprocs = CCTK_nProcs(GH);
+  myproc = CCTK_MyProc(GH);
+#endif
+
+
+  /* prevent compiler warnings about unused parameters */
+  (void) (GH + 0);
+  
+  /* allocate memory for type codes and input array pointer */
+  input_array_type_codes = (CCTK_INT *)malloc (N_input_arrays *sizeof(CCTK_INT));
+  input_arrays           = (CCTK_POINTER) malloc (N_input_arrays *sizeof(CCTK_POINTER));
+ 
+  /* allocate memory for the array of the dimensions of the input arrays */
+  lower_array_bounds = (CCTK_INT *)malloc (dim *sizeof(CCTK_INT));
+  upper_array_bounds = (CCTK_INT *)malloc (dim *sizeof(CCTK_INT));
+  input_array_dims = (CCTK_INT *)malloc (dim *sizeof(CCTK_INT));
+
+  /* set the number of dimensions */
+  dim = CCTK_GroupDimFromVarI(input_array_variable_indices[0]);
+
+  /* find out the types of the input arrays and put that */
+  /* in an array */
+  for ( i = 0; i < N_input_arrays; i++)
+  {
+    input_array_type_codes[i] = CCTK_VarTypeI( input_array_variable_indices[i] );
+    input_arrays[i]           = CCTK_VarDataPtrI(GH, 0, input_array_variable_indices[i] );
+  }
+
+  /* fill out the dimensions of the input arrays */
+  ierr = CCTK_GrouplbndVI(GH, dim, lower_array_bounds, CCTK_GroupIndexFromVarI(input_array_variable_indices[0]) );
+  ierr = CCTK_GroupubndVI(GH, dim, upper_array_bounds, CCTK_GroupIndexFromVarI(input_array_variable_indices[0]) );
+
+  /* incement max subscript to use < instead of <= */
+  for ( i=0; i < dim; i++)
+  {
+    upper_array_bounds[i]++;
+  }
+  /* Create the parameter table if it is not there and add the bounds to it */
+  if ( Util_TableQueryNKeys(param_table_handle) == 0)
+  {
+    param_table_handle = Util_TableCreate (UTIL_TABLE_FLAGS_DEFAULT);
+    ierr = Util_TableSetGenericArray (param_table_handle, CCTK_VARIABLE_INT, dim,lower_array_bounds, "input_array_min_subscripts"); 
+    ierr = Util_TableSetGenericArray (param_table_handle, CCTK_VARIABLE_INT, dim,upper_array_bounds, "input_array_max_subscripts"); 
+  }
+  else
+  {
+    if (!Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "input_array_min_subscripts"))
+    {
+      ierr = Util_TableSetGenericArray (param_table_handle, CCTK_VARIABLE_INT, dim,lower_array_bounds, "input_array_min_subscripts"); 
+    }
+    if (!Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "input_array_max_subscripts"))
+    {
+      ierr = Util_TableSetGenericArray (param_table_handle, CCTK_VARIABLE_INT, dim,lower_array_bounds, "input_array_max_subscripts"); 
+    }
+  }
+
+  for (i = 0; i<dim; i++)
+  {
+    input_array_dims[i] = upper_array_bounds[i] - lower_array_bounds[i] + 1;
+  }
+
+  /* Set flag to tell local reduction not to divide by num_points */
+  ierr = Util_TableSetInt(param_table_handle, 1, "global_calling");
+
+  ierr = CCTK_ReduceLocalArrays(dim, local_reduce_handle, param_table_handle,
+                                N_input_arrays, input_array_dims, input_array_type_codes,
+                                (void *)input_arrays, M_output_values, output_value_type_codes, output_values);
+
+  /* Get flag to perform division by number of points or not */
+  ierr = Util_TableGetInt(param_table_handle, &perform_division, "perform_division");
+
+#ifdef CCTK_MPI
+  pughGH = PUGH_pGH (GH);
+  ierr = Util_TableGetInt(param_table_handle, &num_points, "num_points");
+  ierr = Util_TableGetInt(param_table_handle, &proc, "proc");
+  ierr = Util_TableGetInt(param_table_handle, &mpi_operation, "mpi_operation");
+
+  CACTUS_MPI_ERROR (MPI_Allreduce (&num_points, &total_num_points, 1,
+                    PUGH_MPI_INT, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+
+  for (i = 0; i< M_output_values; i++)
+  {
+    switch (output_value_type_codes[i])
+    {
+     /* out values type switches*/
+      case CCTK_VARIABLE_CHAR:
+        local_outvals = malloc (1 * sizeof (CCTK_CHAR));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_CHAR *) output_values[i], 1 * sizeof (CCTK_CHAR));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_CHAR *) output_values[i], 1,
+                        PUGH_MPI_CHAR, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_CHAR *) output_values[i], 1,
+                        PUGH_MPI_CHAR, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_CHAR *) output_values[i], 1,
+                        PUGH_MPI_CHAR, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      case CCTK_VARIABLE_INT:
+        local_outvals = malloc (1 * sizeof (CCTK_INT));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_INT *) output_values[i], 1 * sizeof (CCTK_INT));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT *) output_values[i], 1,
+                        PUGH_MPI_INT, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT *) output_values[i], 1,
+                        PUGH_MPI_INT, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT *) output_values[i], 1,
+                        PUGH_MPI_INT, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #ifdef CCTK_INT1
+      case CCTK_VARIABLE_INT1:
+        local_outvals = malloc (1 * sizeof (CCTK_INT1));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_INT1 *) output_values[i], 1 * sizeof (CCTK_INT1));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT1 *) output_values[i], 1,
+                        PUGH_MPI_INT1, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT1 *) output_values[i], 1,
+                        PUGH_MPI_INT1, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT1 *) output_values[i], 1,
+                        PUGH_MPI_INT1, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif                                                              
+      #ifdef CCTK_INT2
+      case CCTK_VARIABLE_INT2:
+        local_outvals = malloc (1 * sizeof (CCTK_INT2));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_INT2 *) output_values[i], 1 * sizeof (CCTK_INT2));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT2 *) output_values[i], 1,
+                        PUGH_MPI_INT2, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT2 *) output_values[i], 1,
+                        PUGH_MPI_INT2, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT2 *) output_values[i], 1,
+                        PUGH_MPI_INT2, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_INT4
+      case CCTK_VARIABLE_INT4:
+        local_outvals = malloc (1 * sizeof (CCTK_INT4));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_INT4 *) output_values[i], 1 * sizeof (CCTK_INT4));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT4 *) output_values[i], 1,
+                        PUGH_MPI_INT4, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT4 *) output_values[i], 1,
+                        PUGH_MPI_INT4, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT4 *) output_values[i], 1,
+                        PUGH_MPI_INT4, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_INT8
+      case CCTK_VARIABLE_INT8:
+        local_outvals = malloc (1 * sizeof (CCTK_INT8));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_INT8 *) output_values[i], 1 * sizeof (CCTK_INT8));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT8 *) output_values[i], 1,
+                        PUGH_MPI_INT8, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT8 *) output_values[i], 1,
+                        PUGH_MPI_INT8, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_INT8 *) output_values[i], 1,
+                        PUGH_MPI_INT8, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;                                                                
+      #endif
+      case CCTK_VARIABLE_REAL:
+        local_outvals = malloc (1 * sizeof (CCTK_REAL));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_REAL *) output_values[i], 1 * sizeof (CCTK_REAL));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL *) output_values[i], 1,
+                        PUGH_MPI_REAL, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL *) output_values[i], 1,
+                        PUGH_MPI_REAL, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL *) output_values[i], 1,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        local_outvals = malloc (1 * sizeof (CCTK_REAL4));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_REAL4 *) output_values[i], 1 * sizeof (CCTK_REAL4));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL4 *) output_values[i], 1,
+                        PUGH_MPI_REAL4, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL4 *) output_values[i], 1,
+                        PUGH_MPI_REAL4, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL4 *) output_values[i], 1,
+                        PUGH_MPI_REAL4, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        local_outvals = malloc (1 * sizeof (CCTK_REAL8));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_REAL8 *) output_values[i], 1 * sizeof (CCTK_REAL8));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL8 *) output_values[i], 1,
+                        PUGH_MPI_REAL8, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL8 *) output_values[i], 1,
+                        PUGH_MPI_REAL8, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL8 *) output_values[i], 1,
+                        PUGH_MPI_REAL8, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL6:
+        local_outvals = malloc (1 * sizeof (CCTK_REAL16));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_REAL16 *) output_values[i], 1 * sizeof (CCTK_REAL16));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL16 *) output_values[i], 1,
+                        PUGH_MPI_REAL16, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL16 *) output_values[i], 1,
+                        PUGH_MPI_REAL16, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_REAL16 *) output_values[i], 1,
+                        PUGH_MPI_REAL16, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      case CCTK_VARIABLE_COMPLEX:
+        local_outvals = malloc (1 * sizeof (CCTK_COMPLEX));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_COMPLEX *) output_values[i], 1 * sizeof (CCTK_COMPLEX));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #ifdef CCTK_COMPLEX8
+      case CCTK_VARIABLE_COMPLEX8:
+        local_outvals = malloc (1 * sizeof (CCTK_COMPLEX8));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_COMPLEX8 *) output_values[i], 1 * sizeof (CCTK_COMPLEX8));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX8 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex8, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX8 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex8, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX8 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex8, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_COMPLEX16
+      case CCTK_VARIABLE_COMPLEX16:
+        local_outvals = malloc (1 * sizeof (CCTK_COMPLEX16));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_COMPLEX16 *) output_values[i], 1 * sizeof (CCTK_COMPLEX16));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX16 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex16, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX16 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex16, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX16 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex16, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+      #ifdef CCTK_COMPLEX32
+      case CCTK_VARIABLE_COMPLEX32:
+        local_outvals = malloc (1 * sizeof (CCTK_COMPLEX32));
+
+        /* outvals[] contains now the local reduction values */
+        memcpy ( local_outvals, (CCTK_COMPLEX32 *) output_values[i], 1 * sizeof (CCTK_COMPLEX32));
+
+        if (mpi_operation == 1)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX32 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex32, MPI_MAX, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 2)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX32 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex32, MPI_MIN, pughGH->PUGH_COMM_WORLD));         
+        else if (mpi_operation == 3)
+          CACTUS_MPI_ERROR (MPI_Allreduce (local_outvals, (CCTK_COMPLEX32 *) output_values[i], 1,
+                        pughGH->PUGH_mpi_complex32, MPI_SUM, pughGH->PUGH_COMM_WORLD));         
+      break;
+      #endif
+    }
+  }
+    num_points = total_num_points;
+    free (local_outvals);
+#endif
+
+
+  if (perform_division == 0)
+  {
+    for (i = 0; i< M_output_values; i++)
+    {
+      switch (output_value_type_codes[i])
+      {
+       /* out values type switches*/
+        case CCTK_VARIABLE_CHAR:
+          *( (CCTK_CHAR *) output_values[i]) = *( (CCTK_CHAR *) output_values[i]) / num_points; 
+        break;
+        case CCTK_VARIABLE_INT:
+          *( (CCTK_INT *) output_values[i]) = *( (CCTK_INT *) output_values[i]) / num_points; 
+        break;
+        #ifdef CCTK_INT1
+        case CCTK_VARIABLE_INT1:
+          *( (CCTK_INT1 *) output_values[i]) = *( (CCTK_INT1 *) output_values[i]) / num_points; 
+        break;
+        #endif                                                              
+        #ifdef CCTK_INT2
+        case CCTK_VARIABLE_INT2:
+          *( (CCTK_INT2 *) output_values[i]) = *( (CCTK_INT2 *) output_values[i]) / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_INT4
+        case CCTK_VARIABLE_INT4:
+          *( (CCTK_INT4 *) output_values[i]) = *( (CCTK_INT4 *) output_values[i]) / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_INT8
+        case CCTK_VARIABLE_INT8:
+          *( (CCTK_INT8 *) output_values[i]) = *( (CCTK_INT8 *) output_values[i]) / num_points; 
+        break;                                                                
+        #endif
+        case CCTK_VARIABLE_REAL:
+          *( (CCTK_REAL *) output_values[i]) = *( (CCTK_REAL *) output_values[i]) / num_points; 
+        break;
+        #ifdef CCTK_REAL4
+        case CCTK_VARIABLE_REAL4:
+          *( (CCTK_REAL4 *) output_values[i]) = *( (CCTK_REAL4 *) output_values[i]) / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_REAL8
+        case CCTK_VARIABLE_REAL8:
+          *( (CCTK_REAL8 *) output_values[i]) = *( (CCTK_REAL8 *) output_values[i]) / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_REAL16
+        case CCTK_VARIABLE_REAL6:
+          *( (CCTK_REAL16 *) output_values[i]) = *( (CCTK_REAL16 *) output_values[i]) / num_points; 
+        break;
+        #endif
+        case CCTK_VARIABLE_COMPLEX:
+          (*( (CCTK_COMPLEX *) output_values[i])).Re = (*( (CCTK_COMPLEX *) output_values[i])).Re / num_points; 
+          (*( (CCTK_COMPLEX *) output_values[i])).Im = (*( (CCTK_COMPLEX *) output_values[i])).Im / num_points; 
+        break;
+        #ifdef CCTK_COMPLEX8
+        case CCTK_VARIABLE_COMPLEX8:
+          (*( (CCTK_COMPLEX8 *) output_values[i])).Re = (*( (CCTK_COMPLEX8 *) output_values[i])).Re / num_points; 
+          (*( (CCTK_COMPLEX8 *) output_values[i])).Im = (*( (CCTK_COMPLEX8 *) output_values[i])).Im / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_COMPLEX16
+        case CCTK_VARIABLE_COMPLEX16:
+          (*( (CCTK_COMPLEX16 *) output_values[i])).Re = (*( (CCTK_COMPLEX16 *) output_values[i])).Re / num_points; 
+          (*( (CCTK_COMPLEX16 *) output_values[i])).Im = (*( (CCTK_COMPLEX16 *) output_values[i])).Im / num_points; 
+        break;
+        #endif
+        #ifdef CCTK_COMPLEX32
+        case CCTK_VARIABLE_COMPLEX32:
+          (*( (CCTK_COMPLEX32 *) output_values[i])).Re = (*( (CCTK_COMPLEX32 *) output_values[i])).Re / num_points; 
+          (*( (CCTK_COMPLEX32 *) output_values[i])).Im = (*( (CCTK_COMPLEX32 *) output_values[i])).Im / num_points; 
+        break;
+        #endif
+      }
+    }
+  }
+
+  return ierr;
+}
diff --git a/src/Startup.c b/src/Startup.c
index e40f2b9..872e27f 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -5,7 +5,7 @@
    @desc
               Startup routines for PUGHReduce.
    @enddesc
-   @version   $Id$
+   @version   $Header$
 @@*/
 
 #include "cctk.h"
@@ -13,7 +13,7 @@
 
 static const char *rcsid="$Header$";
 
-CCTK_FILEVERSION(CactusPUGH_PUGHReduce_Startup_c)
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_Startup_c);
 
 
 /* prototypes of routines defined in this source file */
@@ -37,6 +37,9 @@ int PUGHReduce_Startup(void);
 @@*/
 int PUGHReduce_Startup (void)
 {
+  /* Register the Global reduction operator */
+  CCTK_RegisterGridArrayReductionOperator(PUGH_ReduceGridArrays);
+  
   /* Register the reduction operators provided by PUGH */
   CCTK_RegisterReductionOperator (PUGH_ReductionAvgGVs,     "average");
   CCTK_RegisterReductionOperator (PUGH_ReductionAvgGVs,     "mean");
diff --git a/src/make.code.defn b/src/make.code.defn
index 7e73b26..1ebc323 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -4,4 +4,4 @@
 # Source files in this directory
 SRCS = ReductionAvg.c ReductionCount.c ReductionMax.c ReductionMin.c          \
        ReductionNorm1.c ReductionNorm2.c ReductionNorm3.c ReductionNorm4.c    \
-       ReductionNormInf.c ReductionSum.c Startup.c Reduction.c
+       ReductionNormInf.c ReductionSum.c Startup.c Reduction.c ReduceGA.c
diff --git a/src/pugh_reductions.h b/src/pugh_reductions.h
index f11d012..7634497 100644
--- a/src/pugh_reductions.h
+++ b/src/pugh_reductions.h
@@ -9,7 +9,7 @@
  @@*/
 
 #ifndef _PUGH_REDUCTIONS_H_
-#define _PUGH_REDUCTIONS_H_
+#define _PUGH_REDUCTIONS_H_ 1
 
 
 #include "cctk.h"
@@ -139,6 +139,7 @@
 #ifdef __cplusplus
 extern "C" {
 #endif
+int PUGH_ReduceGridArrays (REDUCTION_GRID_ARRAY_OPERATOR_REGISTER_ARGLIST);
 
 int PUGH_ReductionAvgGVs     (REDUCTION_OPERATOR_REGISTER_ARGLIST);
 int PUGH_ReductionCountGVs   (REDUCTION_OPERATOR_REGISTER_ARGLIST);