local reduce iterator

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

1 files changed, 232 insertions, 1 deletions

diff --git a/src/local_reductions.h b/src/local_reductions.h
index 56a3625..b223d45 100644
--- a/src/local_reductions.h
+++ b/src/local_reductions.h
@@ -1,7 +1,7 @@
  /*@@
    @header    local_reductions.h
    @date      
-   @author    Gabrielle Allen, Yaakoub El Khamra
+   @author    Yaakoub El Khamra
    @desc
               Prototypes for local reduction operators
    @enddesc
@@ -22,6 +22,237 @@ extern "C" {
 #endif
 
 
+#define ITERATE_ON_ARRAY(i,cctk_type, in_data, out_type, out_num, mask_on,input_array_offset, indices, sum_indices, max_iter, iter, flag, actual_indices,input_array_strides, input_array_min_subscripts,product)           \
+{                                                                                 \
+      iter = 0;                                                                   \
+      sum_indices = 0;                                                            \
+  const cctk_type *typed_vdata = (cctk_type *)(in_data);                          \
+  out_type * outval = (out_type *)(out_num);                                      \
+  REDUCTION_INITIAL( *outval)                                                     \
+  if ( mask_on == 1)                                                              \
+  {                                                                               \
+    if ( input_array_offset == 0)                                                 \
+    {                                                                             \
+      while (iter < max_iter)                                                     \
+      {                                                                           \
+        REDUCTION_OPERATION(* outval, typed_vdata[sum_indices])                   \
+        num_points++;                                                             \
+        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                                                                          \
+    {                                                                             \
+      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;                                                              \
+        }                                                                         \
+        num_points++;                                                             \
+        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_offset;                       \
+        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_offset == 0)                                                 \
+    {                                                                             \
+      iter = 0;                                                                   \
+      sum_indices = 0;                                                            \
+      while (iter < max_iter)                                                     \
+      {                                                                           \
+        num_points++;                                                             \
+        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                                                                          \
+    {                                                                             \
+      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;                                                              \
+        }                                                                         \
+        num_points++;                                                             \
+        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_offset;                       \
+        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");                          \
+  }                                                                               \
+}
+  
+
 int LocalReduce_Mean     (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int LocalReduce_Count   (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);
 int LocalReduce_Max  (REDUCTION_LOCAL_ARRAY_OPERATOR_REGISTER_ARGLIST);