started putting stuff back

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

3 files changed, 861 insertions, 338 deletions

diff --git a/src/ReductionAvg.c b/src/ReductionAvg.c
index 976efd1..fb466a2 100644
--- a/src/ReductionAvg.c
+++ b/src/ReductionAvg.c
@@ -19,18 +19,14 @@ static const char *rcsid = "$Id$";
 CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionAvg_c);
 
 /* local function prototypes */
-static int ReductionAvg (int num_dims,
-                         const int from[/* dim */],
-                         const int to[/* dim */],
-                         int iterator[/* dim */],
-                         const int points_per_dim[/* dim */],
-                         int num_points,
-                         int have_local_points,
-                         int num_inarrays,
-                         const int intypes[/* num_inarrays */],
-                         const void *const inarrays[/* num_inarrays */],
-                         int num_outvals,
-                         CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+static int ReductionAvg (int N_dims, int operator_handle, 
+                        int param_table_handle,   int N_input_arrays,
+                        const CCTK_INT input_array_dims[], 
+                        const CCTK_INT input_array_type_codes[],
+                        const void *const input_arrays[],
+                        int M_output_numbers,
+                        const CCTK_INT output_number_type_codes[],
+                        void * const output_numbers[]);
 
 
 /*@@
@@ -99,7 +95,7 @@ int LocalReduce_Mean (int N_dims, int operator_handle,
                       const void *const input_arrays[],
                       int M_output_numbers,
                       const CCTK_INT output_number_type_codes[],
-                      void *const output_numbers[])
+                      void * const output_numbers[])
 {
   return (LocalReduce_Reduce (N_dims, operator_handle, 
                           param_table_handle, N_input_arrays,
@@ -109,6 +105,8 @@ int LocalReduce_Mean (int N_dims, int operator_handle,
                           ReductionAvg));
 }
 
+
+
 /*****************************************************************************/
 /*                             local functions                               */
 /*****************************************************************************/
@@ -121,128 +119,439 @@ int LocalReduce_Mean (int N_dims, int operator_handle,
                (num_outvals == 1) or on the results of the local reductions.
    @enddesc
 @@*/
-static int ReductionAvg (int num_dims,
-                         const int from[/* dim */],
-                         const int to[/* dim */],
-                         int iterator[/* dim */],
-                         const int points_per_dim[/* dim */],
-                         int num_points,
-                         int have_local_points,
-                         int num_inarrays,
-                         const int intypes[/* num_inarrays */],
-                         const void *const inarrays[/* num_inarrays */],
-                         int num_outvals,
-                         CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+static int ReductionAvg (int N_dims, int operator_handle, 
+                        int param_table_handle,   int N_input_arrays,
+                        const CCTK_INT input_array_dims[], 
+                        const CCTK_INT input_array_type_codes[],
+                        const void *const input_arrays[],
+                        int M_output_numbers,
+                        const CCTK_INT output_number_type_codes[],
+                        void * const output_numbers[])
 {
-  int i, total_outvals;
-  const char *vtypename;
+  /* utility variables */
+  int i, j, k, flag, product, num_points;
+  int ierr;
+  int * iters_per_dim;
+  void * data_pointer;
+  void * output_buffer;
+
+  output_buffer = (CCTK_REAL *) malloc (M_output_numbers*sizeof(CCTK_REAL));
+
+  /* indices to hold the temp indices of size N_dims  and iteration indices*/
+  int * indices;
+  int * actual_indices;
+  int * actual_iters_per_dim;
+  int index = 0;
+  int max_iter = 0;;
+  int iter = 0;
+  int sum_indices = 0;
+  int max_index = 1;
 
+  /* data pointer offset and strides declared here */
+  CCTK_INT * input_array_offsets;
+  CCTK_INT * input_array_strides;
+  CCTK_INT * input_array_min_subscripts;
+  CCTK_INT * input_array_max_subscripts;
 
-  /* avoid compiler warnings about unused parameters */
-  (void) (num_points + 0);
+  /* excesion variables declared here */
+  int mask_on = 1;          /* mask is by default off=1 */
+  void * mask_array;        /* same dimensions/indexing as input arrays */
+  CCTK_INT mask_type_code;  /* one of the CCTK_VARIABLE_* codes */
+  CCTK_INT mask_offset;
+  CCTK_INT mask_time_level;
 
-/* macros to complete the ITERATE_ARRAY macro */
-#define INITIAL_REDUCTION_VALUE(array)     0
-#define REDUCTION_OPERATION(avg, scalar)   avg += scalar
+  /* set the number of points */
+  num_points = 0;
 
-  for (i = total_outvals = 0; i < num_inarrays; i++)
+  /* allocate memory for iters_per_dim */
+  iters_per_dim = (int *)malloc(N_dims * sizeof(int));
+
+  /* allocate then initialize the values of the strides and subscripts */
+  indices = (int *)malloc (N_dims * sizeof(int));
+  actual_indices = (int *)malloc (N_dims * sizeof(int));
+  actual_iters_per_dim = (int *)malloc (N_dims * sizeof(int));
+
+  /* allocate then initialize the values of the strides and subscripts */
+  input_array_offsets = (CCTK_INT *)malloc (N_input_arrays * sizeof(CCTK_INT));
+  input_array_strides = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+  input_array_min_subscripts = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+  input_array_max_subscripts = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+
+  for (i = 0; i<N_input_arrays; i++)
   {
-    switch (intypes[i])
-    {
-      case CCTK_VARIABLE_CHAR:
-        ITERATE_ARRAY (CCTK_BYTE, num_dims, inarrays[i],
-                       from, to, iterator, points_per_dim,
-                       CCTK_BYTE, outvals, num_outvals, total_outvals);
-        break;
+    input_array_offsets[i] = 0;
+  }
 
-      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;
+  for (i = 0; i<N_dims; i++)
+  {
+    input_array_strides[i]        = 1;
+    input_array_min_subscripts[i] = 0;
+    input_array_max_subscripts[i] = input_array_dims[i];
+    max_index *= input_array_max_subscripts[i];
+  }
 
-#ifdef CCTK_INT1
-      case CCTK_VARIABLE_INT1:
-        ITERATE_ARRAY (CCTK_INT1, num_dims, inarrays[i],
-                       from, to, iterator, points_per_dim,
-                       CCTK_INT1, outvals, num_outvals, total_outvals);
-        break;
-#endif
+  /* for strides and subscripts get values from param table (it they exist) */
+  if ( Util_TableQueryNKeys(param_table_handle) != 0)
+  {
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_input_arrays, input_array_offsets, "input_array_offsets");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_strides, "input_array_strides");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_min_subscripts, "input_array_min_subscripts");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_min_subscripts, "input_array_min_subscripts");  
+  }
 
-#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
+  /* for masks get values from param table (it they exist) */
+  if ( Util_TableQueryNKeys(param_table_handle) != 0)
+  {
+    ierr = 0;
+    ierr = Util_TableGetInt(param_table_handle, mask_type_code, "mask_type_code"); 
 
-#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
+     /*   mask_valid_min, mask_valid_max; 
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_range, "mask_valid_min");
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_range, "mask_valid_max");
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_array, "mask_array");  */ 
+  }
+  /* reduction maps an array to a single value of the same type */
+  if (M_output_numbers != N_input_arrays)
+  {
+    CCTK_WARN (1, "Average reduction returns a single value\n                 \
+                   for each input array\n");
+    return (-1);
+  }
 
-#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);
+  /* set the indices to their minimum values */
+  max_iter = 1;
+  for (j = 0; j <N_dims; j++)
+  {
+    indices       [j]        = 0;
+    actual_indices[j]        = input_array_min_subscripts[j];
+    actual_iters_per_dim [j] = (int) (input_array_max_subscripts[j] - input_array_min_subscripts[j]);
+    iters_per_dim [j]        = (int) ((input_array_max_subscripts[j] - input_array_min_subscripts[j])/input_array_strides[j]);
+    max_iter *= iters_per_dim [j];
+  }
+                                                                                
+  for (i = 0; i < N_input_arrays; i++)
+  { 
+/*
+      switch (output_number_type_codes[i])
+      {
+        case CCTK_VARIABLE_CHAR:
+        output_buffer = (CCTK_BYTE *) output_buffer;
+        output_buffer = (CCTK_BYTE *) malloc (M_output_numbers*sizeof(CCTK_BYTE));
         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);
+        case CCTK_VARIABLE_INT:
+        output_buffer = (CCTK_INT *) output_buffer;
+        output_buffer = (CCTK_INT *) malloc (M_output_numbers*sizeof(CCTK_INT));
         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);
+        #ifdef CCTK_INT1
+        case CCTK_VARIABLE_INT1:
+        output_buffer = (CCTK_INT1 *) output_buffer;
+        output_buffer = (CCTK_INT1 *) malloc (M_output_numbers*sizeof(CCTK_INT1));
         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);
+        #endif                                                              
+        #ifdef CCTK_INT2
+        case CCTK_VARIABLE_INT2:
+        output_buffer = (CCTK_INT2 *) output_buffer;
+        output_buffer = (CCTK_INT2 *) malloc (M_output_numbers*sizeof(CCTK_INT2));
         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);
+        #endif
+        #ifdef CCTK_INT4
+        case CCTK_VARIABLE_INT4:
+        output_buffer = (CCTK_INT4 *) output_buffer;
+        output_buffer = (CCTK_INT4 *) malloc (M_output_numbers*sizeof(CCTK_INT4));
+        break;
+        #endif
+        #ifdef CCTK_INT8
+        case CCTK_VARIABLE_INT8:
+        output_buffer = (CCTK_INT8 *) output_buffer;
+        output_buffer = (CCTK_INT8 *) malloc (M_output_numbers*sizeof(CCTK_INT8));
+        break;                                                                
+        #endif
+        case CCTK_VARIABLE_REAL:
+        #define type  CCTK_REAL
+        output_buffer = (CCTK_REAL *) output_buffer;
+        output_buffer = (CCTK_REAL *) malloc (M_output_numbers*sizeof(CCTK_REAL));
+        break;
+        #ifdef CCTK_REAL4
+        case CCTK_VARIABLE_REAL4:
+        output_buffer = (CCTK_REAL4 *) output_buffer;
+        output_buffer = (CCTK_REAL4 *) malloc (M_output_numbers*sizeof(CCTK_REAL4));
         break;
-#endif
+        #endif
+        #ifdef CCTK_REAL8
+        case CCTK_VARIABLE_REAL8:
+        output_buffer = (CCTK_REAL8 *) output_buffer;
+        output_buffer = (CCTK_REAL8 *) malloc (M_output_numbers*sizeof(CCTK_REAL8));
+        break;
+        #endif
+        #ifdef CCTK_REAL6
+        case CCTK_VARIABLE_REAL6:
+        output_buffer = (CCTK_REAL6 *) output_buffer;
+        output_buffer = (CCTK_REAL6 *) malloc (M_output_numbers*sizeof(CCTK_REAL6));
+        break;
+        #endif
+        default:   
+        CCTK_WARN (1, "Average reduction: Unknown variable type");
+        return (-1);
+      }*/
+    if ( mask_on == 1)
+    {
+      if ( input_array_offsets[i] == 0)
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
 
-      default:
-        vtypename = CCTK_VarTypeName (intypes[i]);
-        if (vtypename && strncmp (vtypename, "CCTK_VARIABLE_COMPLEX", 21) == 0)
+        while (iter < max_iter)
         {
-          CCTK_WARN (1, "ReductionAvg: Don't know how to compute "
-                        "the average of complex variables !!!");
+          
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
         }
-        else
-        {
-          CCTK_WARN (1, "ReductionAvg: Unknown variable type");
+      }
+      else
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        { 
+          /* prevent offset from giving segfaults */
+          if (sum_indices >= max_iter)
+          {
+            CCTK_WARN(1,"offsets and strides access unallocated memory");
+            return -1;
+          }
+
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0]+input_array_offsets[i];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
         }
-        return (-1);
+      }             
     }
-  }
+    else if (mask_on == 0)
+    {
+      if ( input_array_offsets[i] == 0)
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
 
-  /* finally assign the return value */
-  for (i = 0; i < total_outvals; i++)
-  {
-    outvals[i] /= num_points;
-  }
+        while (iter < max_iter)
+        {
+          
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }
+      else
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        { 
+          /* prevent offset from giving segfaults */
+          if (sum_indices >= max_iter)
+          {
+            CCTK_WARN(1,"offsets and strides access unallocated memory");
+            return -1;
+          }
+
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
 
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0]+input_array_offsets[i];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }         
+    }  
+    else
+    {
+      CCTK_WARN(1, "mask_on is not set to a valid value");
+    }
+  }
   return (0);
 }
+
+
diff --git a/src/ReductionCount.c b/src/ReductionCount.c
index 0964b04..55585d4 100644
--- a/src/ReductionCount.c
+++ b/src/ReductionCount.c
@@ -19,144 +19,477 @@ static const char *rcsid = "$Id$";
 CCTK_FILEVERSION(CCTDevelopment_LocalReduce_ReductionCount_c);
 
 /* local function prototypes */
-static int ReductionCount (int num_dims,
-                           const int from[/* dim */],
-                           const int to[/* dim */],
-                           int iterator[/* dim */],
-                           const int points_per_dim[/* dim */],
-                           int num_points,
-                           int have_local_points,
-                           int num_inarrays,
-                           const int intypes[/* num_inarrays */],
-                           const void *const inarrays[/* num_inarrays */],
-                           int num_outvals,
-                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+static int ReductionCount (int N_dims, int operator_handle, 
+                        int param_table_handle,   int N_input_arrays,
+                        const CCTK_INT input_array_dims[], 
+                        const CCTK_INT input_array_type_codes[],
+                        const void *const input_arrays[],
+                        int M_output_numbers,
+                        const CCTK_INT output_number_type_codes[],
+                        void * const output_numbers[]);
 
 
 /*@@
-  @routine LocalReduce_Count
+  @routine LocalReduce_Mean
   @author  Thomas Radke
   @date    19 Aug 1999
   @desc
-           Registered reduction routine for computing the counts
-           of a set of arrays.
-           The arrays are described by their dimensions and variable type.
-           For the number of output values only 1 is accepted.
-           Type casting of the result is provided by specifying the
-           requested output datatype. The intermediate reduction value
-           is always computed as a CCTK_REAL value internally.
   @enddesc
-  @var     GH
-  @vdesc   Pointer to CCTK grid hierarchy
-  @vtype   const cGH *
+  @history
+  @endhistory
+  @var     N_dims
+  @vdesc   number of dimensions in the *reduction*
+  @vtype   int
   @vio     in
   @endvar
-  @var     proc
-  @vdesc   processor that should receive the result of operation
-           (negative value means all processors receive the result)
+  @var     operator_handle
+  @vdesc   operator handle specificies the type of reduction we will perform
   @vtype   int
   @vio     in
   @endvar
-  @var     nDims
-  @vdesc   number of dimensions of input arrays
+  @var     param_table_handle
+  @vdesc   handle to "parameter table", a key-value table
   @vtype   int
   @vio     in
   @endvar
-  @var     dims
-  @vdesc   dimensions of input arrays
-  @vtype   const int *
+  @var     N_input_arrays
+  @vdesc   number of input arrays
+  @vtype   int 
   @vio     in
   @endvar
-  @var     nArrays
-  @vdesc   number of input arrays
-  @vtype   int
+  @var     input_array_dims
+  @vdesc   array of input array dimensions (common to all input arrays)
+  @vtype   const CCTK_INT
   @vio     in
   @endvar
-  @var     inarrays
-  @vdesc   field of input arrays
-  @vtype   const void *const
+  @var     input_array_type_codes
+  @vdesc   array of CCTK_VARIABLE_* codes giving data types of input arrays
+  @vtype   const CCTK_INT
   @vio     in
   @endvar
-  @var     inType
-  @vdesc   (common) variable type of input arrays
-  @vtype   int
+  @var     input_arrays
+  @vdesc   array of pointers to input arrays
+  @vtype   const void *const
   @vio     in
   @endvar
-  @var     num_outvals
-  @vdesc   number of values per output array
+  @var     M_output_numbers
+  @vdesc   
   @vtype   int
   @vio     in
   @endvar
-  @var     outvals
-  @vdesc   pointer to buffer holding the output values
-  @vtype   void *
+  @var     output_number_type_codes
+  @vdesc   array of CCTK_VARIABLE_* codes giving data types of output numbers
+  @vtype   const CCTK_INT
   @vio     in
   @endvar
-  @var     outtype
-  @vdesc   (common) variable type of output arrays
-  @vtype   int
+  @var     output_numbers
+  @vdesc   array[M_output_numbers] of pointers to output numbers[M_reduce_numbers]
+  @vtype   void *const
   @vio     in
   @endvar
 @@*/
 int LocalReduce_Count (int N_dims, int operator_handle, 
-                       int param_table_handle,   int N_input_arrays,
-                       const CCTK_INT input_array_dims[], 
-                       const CCTK_INT input_array_type_codes[],
-                       const void *const input_arrays[],
-                       int M_output_numbers,
-                       const CCTK_INT output_number_type_codes[],
-                       void *const output_numbers[])
+                      int param_table_handle,   int N_input_arrays,
+                      const CCTK_INT input_array_dims[], 
+                      const CCTK_INT input_array_type_codes[],
+                      const void *const input_arrays[],
+                      int M_output_numbers,
+                      const CCTK_INT output_number_type_codes[],
+                      void * const output_numbers[])
 {
   return (LocalReduce_Reduce (N_dims, operator_handle, 
-                           param_table_handle, N_input_arrays,
-                           input_array_dims, input_array_type_codes,
-                           input_arrays, M_output_numbers,
-                           output_number_type_codes, output_numbers,
-                           ReductionCount));
+                          param_table_handle, N_input_arrays,
+                          input_array_dims, input_array_type_codes,
+                          input_arrays, M_output_numbers,
+                          output_number_type_codes, output_numbers,
+                          ReductionAvg));
 }
 
+
+
 /*****************************************************************************/
 /*                             local functions                               */
 /*****************************************************************************/
 /*@@
-   @routine    ReductionCount
+   @routine    ReductionAvg
    @date       Aug 19 1999
    @author     Thomas Radke
-   @desc
-               Returns the number of grid points of a distributed array
+   @desc       Returns the average 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 ReductionCount (int num_dims,
-                           const int from[/* dim */],
-                           const int to[/* dim */],
-                           int iterator[/* dim */],
-                           const int points_per_dim[/* dim */],
-                           int num_points,
-                           int have_local_points,
-                           int num_inarrays,
-                           const int intypes[/* num_inarrays */],
-                           const void *const inarrays[/* num_inarrays */],
-                           int num_outvals,
-                           CCTK_REAL outvals[/*num_inarrays*num_outvals*/])
+static int ReductionCount (int N_dims, int operator_handle, 
+                        int param_table_handle,   int N_input_arrays,
+                        const CCTK_INT input_array_dims[], 
+                        const CCTK_INT input_array_type_codes[],
+                        const void *const input_arrays[],
+                        int M_output_numbers,
+                        const CCTK_INT output_number_type_codes[],
+                        void * const output_numbers[])
 {
-  int i;
-
-  /* avoid compiler warnings about unused parameters */
-  (void) (num_dims + 0);
-  (void) (from + 0);
-  (void) (to + 0);
-  (void) (iterator + 0);
-  (void) (points_per_dim + 0);
-  (void) (have_local_points + 0);
-  (void) (intypes + 0);
-  (void) (inarrays + 0);
-  (void) (num_outvals + 0);
-
-  /* assign the return value */
-  for (i = 0; i < num_inarrays; i++)
+  /* utility variables */
+  int i, j, k, flag, product, num_points;
+  int ierr;
+  int * iters_per_dim;
+  void * data_pointer;
+  void * output_buffer;
+
+  output_buffer = (CCTK_REAL *) malloc (M_output_numbers*sizeof(CCTK_REAL));
+
+  /* indices to hold the temp indices of size N_dims  and iteration indices*/
+  int * indices;
+  int * actual_indices;
+  int * actual_iters_per_dim;
+  int index = 0;
+  int max_iter = 0;;
+  int iter = 0;
+  int sum_indices = 0;
+  int max_index = 1;
+
+  /* data pointer offset and strides declared here */
+  CCTK_INT * input_array_offsets;
+  CCTK_INT * input_array_strides;
+  CCTK_INT * input_array_min_subscripts;
+  CCTK_INT * input_array_max_subscripts;
+
+  /* excesion variables declared here */
+  int mask_on = 1;          /* mask is by default off=1 */
+  void * mask_array;        /* same dimensions/indexing as input arrays */
+  CCTK_INT mask_type_code;  /* one of the CCTK_VARIABLE_* codes */
+  CCTK_INT mask_offset;
+  CCTK_INT mask_time_level;
+
+  /* set the number of points */
+  num_points = 0;
+
+  /* allocate memory for iters_per_dim */
+  iters_per_dim = (int *)malloc(N_dims * sizeof(int));
+
+  /* allocate then initialize the values of the strides and subscripts */
+  indices = (int *)malloc (N_dims * sizeof(int));
+  actual_indices = (int *)malloc (N_dims * sizeof(int));
+  actual_iters_per_dim = (int *)malloc (N_dims * sizeof(int));
+
+  /* allocate then initialize the values of the strides and subscripts */
+  input_array_offsets = (CCTK_INT *)malloc (N_input_arrays * sizeof(CCTK_INT));
+  input_array_strides = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+  input_array_min_subscripts = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+  input_array_max_subscripts = (CCTK_INT *) malloc (N_dims * sizeof(CCTK_INT));
+
+  for (i = 0; i<N_input_arrays; i++)
+  {
+    input_array_offsets[i] = 0;
+  }
+
+  for (i = 0; i<N_dims; i++)
+  {
+    input_array_strides[i]        = 1;
+    input_array_min_subscripts[i] = 0;
+    input_array_max_subscripts[i] = input_array_dims[i];
+    max_index *= input_array_max_subscripts[i];
+  }
+
+  /* for strides and subscripts get values from param table (it they exist) */
+  if ( Util_TableQueryNKeys(param_table_handle) != 0)
   {
-    outvals[i] = num_points;
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_input_arrays, input_array_offsets, "input_array_offsets");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_strides, "input_array_strides");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_min_subscripts, "input_array_min_subscripts");  
+     ierr = Util_TableGetGenericArray(param_table_handle, CCTK_VARIABLE_INT, 
+              N_dims, input_array_min_subscripts, "input_array_min_subscripts");  
   }
 
+  /* for masks get values from param table (it they exist) */
+  if ( Util_TableQueryNKeys(param_table_handle) != 0)
+  {
+    ierr = 0;
+    ierr = Util_TableGetInt(param_table_handle, mask_type_code, "mask_type_code"); 
+
+     /*   mask_valid_min, mask_valid_max; 
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_range, "mask_valid_min");
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_range, "mask_valid_max");
+    ierr = Util_TableGetGeneric(param_table_handle, mask_type_code, 
+              mask_array, "mask_array");  */ 
+  }
+  /* reduction maps an array to a single value of the same type */
+  if (M_output_numbers != N_input_arrays)
+  {
+    CCTK_WARN (1, "Average reduction returns a single value\n                 \
+                   for each input array\n");
+    return (-1);
+  }
+
+  /* set the indices to their minimum values */
+  max_iter = 1;
+  for (j = 0; j <N_dims; j++)
+  {
+    indices       [j]        = 0;
+    actual_indices[j]        = input_array_min_subscripts[j];
+    actual_iters_per_dim [j] = (int) (input_array_max_subscripts[j] - input_array_min_subscripts[j]);
+    iters_per_dim [j]        = (int) ((input_array_max_subscripts[j] - input_array_min_subscripts[j])/input_array_strides[j]);
+    max_iter *= iters_per_dim [j];
+  }
+                                                                                
+  for (i = 0; i < N_input_arrays; i++)
+  { 
+    if ( mask_on == 1)
+    {
+      if ( input_array_offsets[i] == 0)
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        {
+          
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }
+      else
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        { 
+          /* prevent offset from giving segfaults */
+          if (sum_indices >= max_iter)
+          {
+            CCTK_WARN(1,"offsets and strides access unallocated memory");
+            return -1;
+          }
+
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0]+input_array_offsets[i];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }             
+    }
+    else if (mask_on == 0)
+    {
+      if ( input_array_offsets[i] == 0)
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        {
+          
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }
+      else
+      {
+        output_buffer = 0;
+        iter = 0;
+        sum_indices = 0;
+
+        while (iter < max_iter)
+        { 
+          /* prevent offset from giving segfaults */
+          if (sum_indices >= max_iter)
+          {
+            CCTK_WARN(1,"offsets and strides access unallocated memory");
+            return -1;
+          }
+
+          ((CCTK_REAL *)output_buffer)[i] += ((CCTK_REAL* )data_pointer)[sum_indices];
+          num_points++;
+
+          printf("\nvalue is:%f",((CCTK_REAL* )data_pointer)[sum_indices]);
+
+          flag = 0;
+          for (k=0;k<N_dims;k++)
+          {
+            if( indices[k] < iters_per_dim[k]-1)
+            {
+              if (flag==1)
+              {
+                actual_indices[k] += input_array_strides[k-1];
+                indices[k]++;
+                iter++;
+                flag = 0;
+                break;
+              }
+              indices[k]++;
+              actual_indices[k] += input_array_strides[k];
+              iter++;
+              break;
+            }
+            else if (indices[k] == iters_per_dim[k]-1)
+            { 
+              indices[k] = 0;
+              actual_indices[k] = input_array_min_subscripts[k];
+              flag = 1;
+              continue;
+            }
+            else
+            {
+              CCTK_WARN(1,"indices out of bounds, this should not happen");
+              return -1;
+            }
+          }
+          sum_indices = actual_indices[0]+input_array_offsets[i];
+          for (k=N_dims-1;k>0;k--)
+          {
+            product = 1;           
+            for (j=k-1;j>=0;j--)
+            {
+              product *= actual_iters_per_dim[j];
+            }
+            sum_indices += actual_indices[k]*product;
+          }          
+        }
+      }         
+    }  
+    else
+    {
+      CCTK_WARN(1, "mask_on is not set to a valid value");
+    }
+  }
   return (0);
 }
+
+
diff --git a/src/local_reductions.h b/src/local_reductions.h
index 07e51d7..56a3625 100644
--- a/src/local_reductions.h
+++ b/src/local_reductions.h
@@ -16,126 +16,6 @@
 #include "cctk_Reduction.h"
 
 
-/***
-   Macro to iterate over every element of an arbitrary sized array
-
-   The array is given by
-     - its CCTK type, cctk_type,
-     - the number of dimensions, vdim,
-     - an (untyped) pointer to its data, vdata,
-     - one-dimensional vectors, from[] and to[], of length vdim,
-       specifying the start- and endpoint to iterate over for each dimension
-
-   The result of the reduction is specified by
-     - the array where to store the results, outvals[]
-     - the type of the intermediate result scalar, outvals_type
-     - the number of scalars to reduce from this input array, num_outvals
-     - the total number of scalars for all input arrays, total_outvals
-       (incremented after each iteration)
-
-   An iteration vector, iterator[], is passed in to hold the current indices
-   for each dimension of the array.
-   Another vector of size vdim, points_per_dim[], is used to compute the
-   linear index out of the iterator[] elements.
-
-   Before the iteration loop, a local reduction variable, typed_outval,
-   is defined and initialized by the macro INITIAL_REDUCTION_VALUE(typed_vdata).
-
-   For each iteration the macro REDUCTION_OPERATION(typed_outval, typed_inval)
-   is called, with the local reduction variable and the pointer to the
-   current data element.
-
-   Finally, the local reduction variable, typed_outval, is type-casted
-   to CCTK_REAL and assigned to the current output array value outvals[].
-   The total output scalar counter, total_outvals, is incremented.
-
-
-   NOTE: The macro assumes from[i] > to[i] for all 1 <= i < vdim such that
-         at least one iteration will be performed.
-         Make sure this is true when you set up the from[] and to[] vectors.
-
-         For efficiency reasons, the innermost loop (the one for dimension 0)
-         is coded as a for() loop. Thus the overhead to iterate over a
-         one-dimensional array should be kept to a minimum.
-***/
-#define ITERATE_ARRAY(cctk_type, vdim, vdata, from, to, iterator,             \
-                      points_per_dim,                                         \
-                      outvals_type, outvals, num_outvals, total_outvals)      \
-    {                                                                         \
-      int _i, _j, _dim, _vindex;                                              \
-      const cctk_type *typed_vdata = (vdata);                                 \
-      outvals_type typed_outval;                                              \
-                                                                              \
-      if (have_local_points)                                                  \
-      {                                                                       \
-        for (_j = 0; _j < num_outvals; _j++, typed_vdata++)                   \
-        {                                                                     \
-          /* get the linear index of the element to start with */             \
-          _vindex = from[0];                                                  \
-          for (_i = 1; _i < vdim; _i++)                                       \
-            _vindex += from [_i] * points_per_dim [_i];                       \
-  CCTK_VInfo(CCTK_THORNSTRING, "initial reduction value from and to are: %d ,%d ", from[0], to[0]);  \
-          typed_outval = INITIAL_REDUCTION_VALUE(typed_vdata + _vindex);      \
-  CCTK_VInfo(CCTK_THORNSTRING, "initial reduction vars value  are: %d ,%d ", typed_vdata, _vindex);  \
-                                                                              \
-          /* set iterator to local startpoint */                              \
-          memcpy (iterator, from, vdim * sizeof (int));                       \
-                                                                              \
-          /* do the nested loops starting with the second-innermost */        \
-          _dim = 1;                                                           \
-          while (1)                                                           \
-          {                                                                   \
-            /* get the linear index */                                        \
-            _vindex = 0;                                                      \
-            for (_i = 1; _i < vdim; _i++)                                     \
-              _vindex += iterator [_i] * points_per_dim [_i];                 \
-                                                                              \
-            /* do the reduction for the innermost loop (lowest dimension) */  \
-            for (_i = from [0]; _i < to [0]; _i++)                            \
-            {                                                                 \
-              REDUCTION_OPERATION (typed_outval, typed_vdata [_i + _vindex]); \
-            }                                                                 \
-                                                                              \
-            if (vdim > 1)                                                     \
-            {                                                                 \
-              /* increment current looper and check for end */                \
-              if (++iterator [_dim] >= to [_dim])                             \
-              {                                                               \
-                /* increment outermost loopers */                             \
-                for (_dim++; _dim < vdim; _dim++)                             \
-                {                                                             \
-                  if (++iterator [_dim] < to [_dim])                          \
-                    break;                                                    \
-                }                                                             \
-                                                                              \
-                /* done if beyond outermost loop */                           \
-                if (_dim >= vdim)                                             \
-                  break;                                                      \
-                                                                              \
-                /* reset innermost loopers */                                 \
-                for (_dim--; _dim >= 0; _dim--)                               \
-                  iterator [_dim] = from [_dim];                              \
-                _dim = 1;                                                     \
-              }                                                               \
-            }                                                                 \
-            else                                                              \
-            {                                                                 \
-              /* exit loop if array is one-dimensional */                     \
-              break;                                                          \
-            }                                                                 \
-                                                                              \
-          } /* end of nested loops over all dimensions */                     \
-                                                                              \
-          outvals [total_outvals++] = (CCTK_REAL) typed_outval;               \
-                                                                              \
-        } /* end of loop over num_outvals */                                  \
-      }                                                                       \
-      else                                                                    \
-      {                                                                       \
-        total_outvals += num_outvals;                                         \
-      }                                                                       \
-    }
-
 
 #ifdef __cplusplus
 extern "C" {
@@ -152,19 +32,19 @@ int LocalReduce_L3   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int LocalReduce_L4   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int LocalReduce_LInf (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int LocalReduce_Sum     (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
-
-typedef int (*reduction_fn_t) (int num_dims,
-                               const int from[/* dim */],
-                               const int to[/* dim */],
-                               int iterator[/* dim */],
-                               const int points_per_dim[/* dim */],
-                               int num_points,
-                               int have_local_points,
-                               int num_inarrays,
-                               const int intypes[/* num_inarrays */],
-                               const void *const inarrays[/* num_inarrays */],
-                               int num_outvals,
-                               CCTK_REAL outvals[/*num_inarrays*num_outvals*/]);
+int LocalReduce_CmplxMax1 (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int LocalReduce_CmplxMax2 (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int LocalReduce_CmplxMin1 (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+int LocalReduce_CmplxMin2 (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
+
+typedef int (*reduction_fn_t) (int N_dims, int operator_handle, 
+                        int param_table_handle,   int N_input_arrays,
+                        const CCTK_INT input_array_dims[], 
+                        const CCTK_INT input_array_type_codes[],
+                        const void *const input_arrays[],
+                        int M_output_numbers,
+                        const CCTK_INT output_number_type_codes[],
+                        void * const output_numbers[]);
 
 int LocalReduce_Reduce (int N_dims, int operator_handle, 
                           int param_table_handle,   int N_input_arrays,
@@ -173,10 +53,11 @@ int LocalReduce_Reduce (int N_dims, int operator_handle,
                           const void *const input_arrays[],
                           int M_output_numbers,
                           const CCTK_INT output_number_type_codes[],
-                          void *const output_numbers[],reduction_fn_t reduction_fn);
+                          void * const output_numbers[],reduction_fn_t reduction_fn);
 
 #ifdef __cplusplus
 }
 #endif
 
+
 #endif