diff options
Diffstat (limited to 'src/ReductionAvg.c')
| -rw-r--r-- | src/ReductionAvg.c | 535 |
1 files changed, 422 insertions, 113 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); } + + |