This commit fully enables full weighted reduction operations from PUGHReduce

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

1 files changed, 185 insertions, 71 deletions

diff --git a/src/ReduceGA.c b/src/ReduceGA.c
index 8c90619..b643128 100644
--- a/src/ReduceGA.c
+++ b/src/ReduceGA.c
@@ -17,6 +17,8 @@
 #include "util_Table.h"
 #include "pugh_reductions.h"
 
+#define ABS(x)      ((x) < 0 ? -(x) : (x))
+
 static const char *rcsid = "$Id:";
 
 CCTK_FILEVERSION(CactusPUGH_PUGHReduce_ReduceGA_c);
@@ -143,6 +145,13 @@ static int ReduceGridArrays (const cGH *GH,
   pGA *GA;
   const void *data;
 
+  /* weight variables */
+  int  weight_variable_index;
+  int weight_on = 0;
+  CCTK_REAL weight_sum = 0.0;
+  CCTK_REAL total_weight_sum = 0.0;
+  const void * weight;
+
   CCTK_INT * lower_array_bounds;
   CCTK_INT * min_array_subscript;
   CCTK_INT * array_lsh;
@@ -267,15 +276,29 @@ static int ReduceGridArrays (const cGH *GH,
   /* Set flag to tell local reduction not to divide by num_points */
   ierr = Util_TableSetInt(param_table_handle, 1, "global_calling");
 
+  /* Set the weight variable from the weight variable index */
+  ierr   = Util_TableGetInt(param_table_handle, &weight_variable_index, "weight_variable_index");
+  if (ierr == 1)
+  {
+    weight = CCTK_VarDataPtrI(GH, 0, weight_variable_index);
+    ierr   = Util_TableSetPointerToConst(param_table_handle, weight, "weight");
+  }
+  else
+  {
+    weight = NULL;
+    ierr = Util_TableSetInt(param_table_handle, 0, "weight_on");
+  }
+
   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 */
+  /* Get flag to perform division by number of points or not, number_of_points and weight_sum */
   ierr = Util_TableGetInt(param_table_handle, &perform_division, "perform_division");
-
   ierr = Util_TableGetInt(param_table_handle, &num_points, "num_points");
-  
+  ierr = Util_TableGetReal(param_table_handle, &weight_sum, "weight_sum");
+  ierr = Util_TableGetInt(param_table_handle, &weight_on, "weight_on");
+
   /* Get flags to perform root operations                   */
   ierr = Util_TableGetInt(param_table_handle, &perform_2_root, "perform_2_root");
   ierr = Util_TableGetInt(param_table_handle, &perform_3_root, "perform_3_root");
@@ -290,13 +313,16 @@ static int ReduceGridArrays (const cGH *GH,
   {
     CACTUS_MPI_ERROR (MPI_Allreduce (&num_points, &total_num_points, 1,
                     PUGH_MPI_INT, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Allreduce (&weight_sum, &total_weight_sum, 1,
+                    PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
   }
   else
   {
     CACTUS_MPI_ERROR (MPI_Reduce (&num_points, &total_num_points, 1,
                     PUGH_MPI_INT, MPI_SUM, dest_proc, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Reduce (&weight_sum, &total_weight_sum, 1,
+                    PUGH_MPI_INT, MPI_SUM, dest_proc, pughGH->PUGH_COMM_WORLD));
   }
-
   for (i = 0; i< M_output_values; i++)
   {
     switch (output_value_type_codes[i])
@@ -760,82 +786,170 @@ static int ReduceGridArrays (const cGH *GH,
     }
   }
     num_points = total_num_points;
+    weight_sum = total_weight_sum;
     free (local_outvals);
 #endif
 
   if (perform_division == 0)
   {
-    for (i = 0; i< M_output_values; i++)
+    if (weight_on == 1)
     {
-      switch (output_value_type_codes[i])
+      if (ABS(weight_sum) > 1e-12)
       {
-       /* out values type switches*/
-        case CCTK_VARIABLE_BYTE:
-          *( (CCTK_BYTE *) output_values[i]) = *( (CCTK_BYTE *) 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_REAL16:
-          *( (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:
+        for (i = 0; i< M_output_values; i++)
+        {
+          switch (output_value_type_codes[i])
+          {
+           /* out values type switches*/
+            case CCTK_VARIABLE_BYTE:
+              *( (CCTK_BYTE *) output_values[i]) = *( (CCTK_BYTE *) output_values[i]) / weight_sum; 
+            break;
+            case CCTK_VARIABLE_INT:
+              *( (CCTK_INT *) output_values[i]) = *( (CCTK_INT *) output_values[i]) / weight_sum; 
+            break;
+            #ifdef CCTK_INT1
+            case CCTK_VARIABLE_INT1:
+              *( (CCTK_INT1 *) output_values[i]) = *( (CCTK_INT1 *) output_values[i]) / weight_sum; 
+            break;
+            #endif                                                              
+            #ifdef CCTK_INT2
+            case CCTK_VARIABLE_INT2:
+              *( (CCTK_INT2 *) output_values[i]) = *( (CCTK_INT2 *) output_values[i]) / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_INT4
+            case CCTK_VARIABLE_INT4:
+              *( (CCTK_INT4 *) output_values[i]) = *( (CCTK_INT4 *) output_values[i]) / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_INT8
+            case CCTK_VARIABLE_INT8:
+              *( (CCTK_INT8 *) output_values[i]) = *( (CCTK_INT8 *) output_values[i]) / weight_sum; 
+            break;                                                                
+            #endif
+            case CCTK_VARIABLE_REAL:
+              *( (CCTK_REAL *) output_values[i]) = *( (CCTK_REAL *) output_values[i]) / weight_sum; 
+            break;
+            #ifdef CCTK_REAL4
+            case CCTK_VARIABLE_REAL4:
+              *( (CCTK_REAL4 *) output_values[i]) = *( (CCTK_REAL4 *) output_values[i]) / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_REAL8
+            case CCTK_VARIABLE_REAL8:
+              *( (CCTK_REAL8 *) output_values[i]) = *( (CCTK_REAL8 *) output_values[i]) / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_REAL16
+            case CCTK_VARIABLE_REAL16:
+              *( (CCTK_REAL16 *) output_values[i]) = *( (CCTK_REAL16 *) output_values[i]) / weight_sum; 
+            break;
+            #endif
+            case CCTK_VARIABLE_COMPLEX:
+              (*( (CCTK_COMPLEX *) output_values[i])).Re = (*( (CCTK_COMPLEX *) output_values[i])).Re / weight_sum; 
+              (*( (CCTK_COMPLEX *) output_values[i])).Im = (*( (CCTK_COMPLEX *) output_values[i])).Im / weight_sum; 
+            break;
+            #ifdef CCTK_COMPLEX8
+            case CCTK_VARIABLE_COMPLEX8:
+              (*( (CCTK_COMPLEX8 *) output_values[i])).Re = (*( (CCTK_COMPLEX8 *) output_values[i])).Re / weight_sum; 
+              (*( (CCTK_COMPLEX8 *) output_values[i])).Im = (*( (CCTK_COMPLEX8 *) output_values[i])).Im / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_COMPLEX16
+            case CCTK_VARIABLE_COMPLEX16:
+            (*( (CCTK_COMPLEX16 *) output_values[i])).Re = (*( (CCTK_COMPLEX16 *) output_values[i])).Re / weight_sum; 
+            (*( (CCTK_COMPLEX16 *) output_values[i])).Im = (*( (CCTK_COMPLEX16 *) output_values[i])).Im / weight_sum; 
+            break;
+            #endif
+            #ifdef CCTK_COMPLEX32
+            case CCTK_VARIABLE_COMPLEX32:
+              (*( (CCTK_COMPLEX32 *) output_values[i])).Re = (*( (CCTK_COMPLEX32 *) output_values[i])).Re / weight_sum; 
+              (*( (CCTK_COMPLEX32 *) output_values[i])).Im = (*( (CCTK_COMPLEX32 *) output_values[i])).Im / weight_sum; 
+            break;
+            #endif
+          }
+        }
+      }
+      else
+      {
+        CCTK_WARN (1, "The sum of weights in average reduction is zero");
+        return (-1);
+      }
+    }
+    else
+    {
+      for (i = 0; i< M_output_values; i++)
+      {
+        switch (output_value_type_codes[i])
+        {
+         /* out values type switches*/
+          case CCTK_VARIABLE_BYTE:
+            *( (CCTK_BYTE *) output_values[i]) = *( (CCTK_BYTE *) 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_REAL16:
+            *( (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
+          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
+        }
       }
     }
   }