diff options
Diffstat (limited to 'src/local_reductions.h')
| -rw-r--r-- | src/local_reductions.h | 119 |
1 files changed, 115 insertions, 4 deletions
diff --git a/src/local_reductions.h b/src/local_reductions.h index f7451f1..6316063 100644 --- a/src/local_reductions.h +++ b/src/local_reductions.h @@ -16,14 +16,14 @@ #include "util_Table.h" #include "cctk_Reduction.h" - +#define ABS(x) ((x) < 0 ? -(x) : (x)) #ifdef __cplusplus 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,input_array_dims,product) \ +#define ITERATE_ON_ARRAY(i,cctk_type, in_data, out_type, out_num, weight_on, weight, input_array_offset, indices, sum_indices, max_iter, iter, flag, actual_indices,input_array_strides, input_array_min_subscripts,input_array_dims,product) \ { \ const cctk_type * typed_vdata = (const cctk_type *)(in_data); \ out_type inval; \ @@ -40,7 +40,117 @@ extern "C" { } \ else \ { \ - if ( mask_on == 1) \ + if ( weight_on == 1) \ + { \ + if ( input_array_offset == 0) \ + { \ + while (iter < max_iter) \ + { \ + sum_indices = actual_indices[0]; \ + for (k=N_dims-1;k>0;k--) \ + { \ + product = 1; \ + for (j=k-1;j>=0;j--) \ + { \ + product *= input_array_dims[j]; \ + } \ + sum_indices += actual_indices[k]*product; \ + } \ + weight_value = ((CCTK_REAL *) weight)[sum_indices]; \ + \ + REDUCTION_PREOP_CAST(inval, typed_vdata,sum_indices, out_type); \ + WEIGHTED_REDUCTION_OPERATION(*outval,inval,weight_value); \ + num_points++; \ + weight_sum += weight_value; \ + iter++; \ + 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]; \ + indices[k]++; \ + flag = 0; \ + break; \ + } \ + indices[k]++; \ + actual_indices[k] += input_array_strides[k]; \ + 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; \ + } \ + } \ + } \ + } \ + else \ + { \ + while (iter < max_iter) \ + { \ + 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; \ + } \ + /* prevent offset from giving segfaults */ \ + if (sum_indices >= max_iter) \ + { \ + CCTK_WARN(1,"offsets and strides access unallocated memory"); \ + return -1; \ + } \ + REDUCTION_PREOP_CAST(inval, typed_vdata,sum_indices, out_type); \ + WEIGHTED_REDUCTION_OPERATION(*outval,inval,weight_value); \ + num_points++; \ + iter++; \ + 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]++; \ + flag = 0; \ + break; \ + } \ + indices[k]++; \ + actual_indices[k] += input_array_strides[k]; \ + 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; \ + } \ + } \ + } \ + } \ + } \ + else if (weight_on ==0) \ { \ if ( input_array_offset == 0) \ { \ @@ -59,6 +169,7 @@ extern "C" { REDUCTION_PREOP_CAST(inval, typed_vdata,sum_indices, out_type); \ REDUCTION_OPERATION(*outval,inval); \ num_points++; \ + weight_sum += weight_value; \ iter++; \ flag = 0; \ for (k=0;k<N_dims;k++) \ @@ -149,7 +260,7 @@ extern "C" { } \ else \ { \ - CCTK_WARN(1, "mask_on is not set to a valid value"); \ + CCTK_WARN(1, "weight_on is not set to a valid value"); \ } \ } \ EXTRA_STEP(*outval, (out_type)num_points) \ |