diff --git a/.vscode/settings.json b/.vscode/settings.json index cd1d7b1..04036d3 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -43,7 +43,9 @@ "opt_vm.h": "c", "stack": "c", "ios": "c", - "stdint.h": "c" + "stdint.h": "c", + "bitmap_alloc.h": "c", + "buddy.h": "c" }, "C_Cpp.errorSquiggles": "disabled" } \ No newline at end of file diff --git a/allocator/bitmap/HugePage/bitmap_alloc.c b/allocator/bitmap/HugePage/bitmap_alloc.c index 4d2129a..ea98898 100644 --- a/allocator/bitmap/HugePage/bitmap_alloc.c +++ b/allocator/bitmap/HugePage/bitmap_alloc.c @@ -82,7 +82,7 @@ void *HugePageAlloc(void) { if (ftruncate(fd, sz) < 0) { perror("ftruncate"); - close(fd); + close(fd); exit(EXIT_FAILURE); } // if (error != 0 && errno == ENOMEM) diff --git a/benchmarks/benchmarks/XSbench/GridInit.c b/benchmarks/benchmarks/XSbench/GridInit.c new file mode 100644 index 0000000..08a0f27 --- /dev/null +++ b/benchmarks/benchmarks/XSbench/GridInit.c @@ -0,0 +1,230 @@ +#include "XSbench_header.h" + +SimulationData grid_init_do_not_profile( Inputs in, int mype ) +{ + // Structure to hold all allocated simuluation data arrays + SimulationData SD; + + // Keep track of how much data we're allocating + size_t nbytes = 0; + + // Set the initial seed value + uint64_t seed = 42; + + // loop variable + long e = 0; + + //////////////////////////////////////////////////////////////////// + // Initialize Nuclide Grids + //////////////////////////////////////////////////////////////////// + + if(mype == 0) printf("Intializing nuclide grids...\n"); + + // First, we need to initialize our nuclide grid. This comes in the form + // of a flattened 2D array that hold all the information we need to define + // the cross sections for all isotopes in the simulation. + // The grid is composed of "NuclideGridPoint" structures, which hold the + // energy level of the grid point and all associated XS data at that level. + // An array of structures (AOS) is used instead of + // a structure of arrays, as the grid points themselves are accessed in + // a random order, but all cross section interaction channels and the + // energy level are read whenever the gridpoint is accessed, meaning the + // AOS is more cache efficient. + + // Initialize Nuclide Grid + SD.length_nuclide_grid = in.n_isotopes * in.n_gridpoints; + SD.nuclide_grid = (NuclideGridPoint *) malloc( SD.length_nuclide_grid * sizeof(NuclideGridPoint)); + assert(SD.nuclide_grid != NULL); + nbytes += SD.length_nuclide_grid * sizeof(NuclideGridPoint); + for( int i = 0; i < SD.length_nuclide_grid; i++ ) + { + SD.nuclide_grid[i].energy = LCG_random_double(&seed); + SD.nuclide_grid[i].total_xs = LCG_random_double(&seed); + SD.nuclide_grid[i].elastic_xs = LCG_random_double(&seed); + SD.nuclide_grid[i].absorbtion_xs = LCG_random_double(&seed); + SD.nuclide_grid[i].fission_xs = LCG_random_double(&seed); + SD.nuclide_grid[i].nu_fission_xs = LCG_random_double(&seed); + } + + // Sort so that each nuclide has data stored in ascending energy order. + for( int i = 0; i < in.n_isotopes; i++ ) + qsort( &SD.nuclide_grid[i*in.n_gridpoints], in.n_gridpoints, sizeof(NuclideGridPoint), NGP_compare); + + // error debug check + /* + for( int i = 0; i < in.n_isotopes; i++ ) + { + printf("NUCLIDE %d ==============================\n", i); + for( int j = 0; j < in.n_gridpoints; j++ ) + printf("E%d = %lf\n", j, SD.nuclide_grid[i * in.n_gridpoints + j].energy); + } + */ + + + //////////////////////////////////////////////////////////////////// + // Initialize Acceleration Structure + //////////////////////////////////////////////////////////////////// + + if( in.grid_type == NUCLIDE ) + { + SD.length_unionized_energy_array = 0; + SD.length_index_grid = 0; + } + + if( in.grid_type == UNIONIZED ) + { + if(mype == 0) printf("Intializing unionized grid...\n"); + + // Allocate space to hold the union of all nuclide energy data + SD.length_unionized_energy_array = in.n_isotopes * in.n_gridpoints; + SD.unionized_energy_array = (double *) malloc( SD.length_unionized_energy_array * sizeof(double)); + assert(SD.unionized_energy_array != NULL ); + nbytes += SD.length_unionized_energy_array * sizeof(double); + + // Copy energy data over from the nuclide energy grid + for( int i = 0; i < SD.length_unionized_energy_array; i++ ) + SD.unionized_energy_array[i] = SD.nuclide_grid[i].energy; + + // Sort unionized energy array + qsort( SD.unionized_energy_array, SD.length_unionized_energy_array, sizeof(double), double_compare); + + // Allocate space to hold the acceleration grid indices + SD.length_index_grid = SD.length_unionized_energy_array * in.n_isotopes; + SD.index_grid = (int *) malloc( SD.length_index_grid * sizeof(int)); + assert(SD.index_grid != NULL); + nbytes += SD.length_index_grid * sizeof(int); + + // Generates the double indexing grid + int * idx_low = (int *) calloc( in.n_isotopes, sizeof(int)); + assert(idx_low != NULL ); + double * energy_high = (double *) malloc( in.n_isotopes * sizeof(double)); + assert(energy_high != NULL ); + + for( int i = 0; i < in.n_isotopes; i++ ) + energy_high[i] = SD.nuclide_grid[i * in.n_gridpoints + 1].energy; + + for( long e = 0; e < SD.length_unionized_energy_array; e++ ) + { + double unionized_energy = SD.unionized_energy_array[e]; + for( long i = 0; i < in.n_isotopes; i++ ) + { + if( unionized_energy < energy_high[i] ) + SD.index_grid[e * in.n_isotopes + i] = idx_low[i]; + else if( idx_low[i] == in.n_gridpoints - 2 ) + SD.index_grid[e * in.n_isotopes + i] = idx_low[i]; + else + { + idx_low[i]++; + SD.index_grid[e * in.n_isotopes + i] = idx_low[i]; + energy_high[i] = SD.nuclide_grid[i * in.n_gridpoints + idx_low[i] + 1].energy; + } + } + } + + free(idx_low); + free(energy_high); + } + + if( in.grid_type == HASH ) + { + if(mype == 0) printf("Intializing hash grid...\n"); + SD.length_unionized_energy_array = 0; + SD.length_index_grid = in.hash_bins * in.n_isotopes; + SD.index_grid = (int *) malloc( SD.length_index_grid * sizeof(int)); + assert(SD.index_grid != NULL); + nbytes += SD.length_index_grid * sizeof(int); + + double du = 1.0 / in.hash_bins; + + // For each energy level in the hash table + #pragma omp parallel for + for( e = 0; e < in.hash_bins; e++ ) + { + double energy = e * du; + + // We need to determine the bounding energy levels for all isotopes + for( long i = 0; i < in.n_isotopes; i++ ) + { + SD.index_grid[e * in.n_isotopes + i] = grid_search_nuclide( in.n_gridpoints, energy, SD.nuclide_grid + i * in.n_gridpoints, 0, in.n_gridpoints-1); + } + } + } + + //////////////////////////////////////////////////////////////////// + // Initialize Materials and Concentrations + //////////////////////////////////////////////////////////////////// + if(mype == 0) printf("Intializing material data...\n"); + + // Set the number of nuclides in each material + SD.num_nucs = load_num_nucs(in.n_isotopes); + SD.length_num_nucs = 12; // There are always 12 materials in XSBench + + // Intialize the flattened 2D grid of material data. The grid holds + // a list of nuclide indices for each of the 12 material types. The + // grid is allocated as a full square grid, even though not all + // materials have the same number of nuclides. + SD.mats = load_mats(SD.num_nucs, in.n_isotopes, &SD.max_num_nucs); + SD.length_mats = SD.length_num_nucs * SD.max_num_nucs; + + // Intialize the flattened 2D grid of nuclide concentration data. The grid holds + // a list of nuclide concentrations for each of the 12 material types. The + // grid is allocated as a full square grid, even though not all + // materials have the same number of nuclides. + SD.concs = load_concs(SD.num_nucs, SD.max_num_nucs); + SD.length_concs = SD.length_mats; + + // Allocate and initialize replicas +#ifdef AML + // num_nucs + aml_replicaset_hwloc_create(&(SD.num_nucs_replica), + SD.length_num_nucs * sizeof(*(SD.num_nucs)), + HWLOC_OBJ_CORE, + HWLOC_DISTANCES_KIND_FROM_OS | + HWLOC_DISTANCES_KIND_MEANS_LATENCY); + nbytes += (SD.num_nucs_replica)->n * (SD.num_nucs_replica)->size; + aml_replicaset_init(SD.num_nucs_replica, SD.num_nucs); + + // concs + aml_replicaset_hwloc_create(&(SD.concs_replica), + SD.length_concs * sizeof(*(SD.concs)), + HWLOC_OBJ_CORE, + HWLOC_DISTANCES_KIND_FROM_OS | + HWLOC_DISTANCES_KIND_MEANS_LATENCY); + nbytes += (SD.concs_replica)->n * (SD.concs_replica)->size; + aml_replicaset_init(SD.concs_replica, SD.concs); + + // unionized_energy_array + if( in.grid_type == UNIONIZED ){ + aml_replicaset_hwloc_create(&(SD.unionized_energy_array_replica), + SD.length_unionized_energy_array * sizeof(*(SD.unionized_energy_array)), + HWLOC_OBJ_CORE, + HWLOC_DISTANCES_KIND_FROM_OS | + HWLOC_DISTANCES_KIND_MEANS_LATENCY); + nbytes += (SD.unionized_energy_array_replica)->n * (SD.unionized_energy_array_replica)->size; + aml_replicaset_init(SD.unionized_energy_array_replica, SD.unionized_energy_array); + } + + // index grid + if( in.grid_type == UNIONIZED || in.grid_type == HASH ){ + aml_replicaset_hwloc_create(&(SD.index_grid_replica), + SD.length_index_grid * sizeof(*(SD.index_grid)), + HWLOC_OBJ_CORE, + HWLOC_DISTANCES_KIND_FROM_OS | + HWLOC_DISTANCES_KIND_MEANS_LATENCY); + nbytes += (SD.index_grid_replica)->n * (SD.index_grid_replica)->size; + aml_replicaset_init(SD.index_grid_replica, SD.index_grid); + } + + // nuclide grid + aml_replicaset_hwloc_create(&(SD.nuclide_grid_replica), + SD.length_nuclide_grid * sizeof(*(SD.nuclide_grid)), + HWLOC_OBJ_CORE, + HWLOC_DISTANCES_KIND_FROM_OS | + HWLOC_DISTANCES_KIND_MEANS_LATENCY); + nbytes += (SD.nuclide_grid_replica)->n * (SD.nuclide_grid_replica)->size; + aml_replicaset_init(SD.nuclide_grid_replica, SD.nuclide_grid); +#endif + + if(mype == 0) printf("Intialization complete. Allocated %.0lf MB of data.\n", nbytes/1024.0/1024.0 ); + return SD; +} \ No newline at end of file diff --git a/benchmarks/benchmarks/XSbench/Main.c b/benchmarks/benchmarks/XSbench/Main.c index 489d3e7..6498bec 100644 --- a/benchmarks/benchmarks/XSbench/Main.c +++ b/benchmarks/benchmarks/XSbench/Main.c @@ -120,4 +120,4 @@ int main( int argc, char* argv[] ) #endif return is_invalid_result; -} +} \ No newline at end of file diff --git a/benchmarks/benchmarks/XSbench/Makefile b/benchmarks/benchmarks/XSbench/Makefile index 39febd1..4fcfab8 100644 --- a/benchmarks/benchmarks/XSbench/Makefile +++ b/benchmarks/benchmarks/XSbench/Makefile @@ -20,7 +20,7 @@ program = XSBench source = \ Main.c \ io.c \ -Simulation.c \ +Simulations.c \ GridInit.c \ XSutils.c \ Materials.c diff --git a/benchmarks/benchmarks/XSbench/Materials.c b/benchmarks/benchmarks/XSbench/Materials.c index e695d34..5606aeb 100644 --- a/benchmarks/benchmarks/XSbench/Materials.c +++ b/benchmarks/benchmarks/XSbench/Materials.c @@ -1,3 +1,4 @@ + // Material data is hard coded into the functions in this file. // Note that there are 12 materials present in H-M (large or small) @@ -114,4 +115,3 @@ double * load_concs( int * num_nucs, int max_num_nucs ) return concs; } - diff --git a/benchmarks/benchmarks/XSbench/Simulations.c b/benchmarks/benchmarks/XSbench/Simulations.c index e69de29..0676e48 100644 --- a/benchmarks/benchmarks/XSbench/Simulations.c +++ b/benchmarks/benchmarks/XSbench/Simulations.c @@ -0,0 +1,871 @@ +#include "XSbench_header.h" + +//////////////////////////////////////////////////////////////////////////////////// +// BASELINE FUNCTIONS +//////////////////////////////////////////////////////////////////////////////////// +// All "baseline" code is at the top of this file. The baseline code is a simple +// implementation of the algorithm, with only minor CPU optimizations in place. +// Following these functions are a number of optimized variants, +// which each deploy a different combination of optimizations strategies. By +// default, XSBench will only run the baseline implementation. Optimized variants +// must be specifically selected using the "-k " command +// line argument. +//////////////////////////////////////////////////////////////////////////////////// + +unsigned long long run_event_based_simulation(Inputs in, SimulationData SD, int mype) +{ + if( mype == 0) + printf("Beginning event based simulation...\n"); + + //////////////////////////////////////////////////////////////////////////////// + // SUMMARY: Simulation Data Structure Manifest for "SD" Object + // Here we list all heap arrays (and lengths) in SD that would need to be + // offloaded manually if using an accelerator with a seperate memory space + //////////////////////////////////////////////////////////////////////////////// + // int * num_nucs; // Length = length_num_nucs; + // double * concs; // Length = length_concs + // int * mats; // Length = length_mats + // double * unionized_energy_array; // Length = length_unionized_energy_array + // int * index_grid; // Length = length_index_grid + // NuclideGridPoint * nuclide_grid; // Length = length_nuclide_grid + // + // Note: "unionized_energy_array" and "index_grid" can be of zero length + // depending on lookup method. + // + // Note: "Lengths" are given as the number of objects in the array, not the + // number of bytes. + //////////////////////////////////////////////////////////////////////////////// + + + //////////////////////////////////////////////////////////////////////////////// + // Begin Actual Simulation Loop + //////////////////////////////////////////////////////////////////////////////// + unsigned long long verification = 0; + int i = 0; + #pragma omp parallel for schedule(dynamic,100) reduction(+:verification) + for( i = 0; i < in.lookups; i++ ) + { + #ifdef AML + int * num_nucs = aml_replicaset_hwloc_local_replica(SD.num_nucs_replica); + double * concs = aml_replicaset_hwloc_local_replica(SD.concs_replica); + double * unionized_energy_array = aml_replicaset_hwloc_local_replica(SD.unionized_energy_array_replica); + int * index_grid = aml_replicaset_hwloc_local_replica(SD.index_grid_replica); + NuclideGridPoint * nuclide_grid = aml_replicaset_hwloc_local_replica(SD.nuclide_grid_replica); + #else + int * num_nucs = SD.num_nucs; + double * concs = SD.concs; + double * unionized_energy_array = SD.unionized_energy_array; + int * index_grid = SD.index_grid; + NuclideGridPoint * nuclide_grid = SD.nuclide_grid; + #endif + + // Set the initial seed value + uint64_t seed = STARTING_SEED; + + // Forward seed to lookup index (we need 2 samples per lookup) + seed = fast_forward_LCG(seed, 2*i); + + // Randomly pick an energy and material for the particle + double p_energy = LCG_random_double(&seed); + int mat = pick_mat(&seed); + + double macro_xs_vector[5] = {0}; + + // Perform macroscopic Cross Section Lookup + calculate_macro_xs( + p_energy, // Sampled neutron energy (in lethargy) + mat, // Sampled material type index neutron is in + in.n_isotopes, // Total number of isotopes in simulation + in.n_gridpoints, // Number of gridpoints per isotope in simulation + num_nucs, // 1-D array with number of nuclides per material + concs, // Flattened 2-D array with concentration of each nuclide in each material + unionized_energy_array, // 1-D Unionized energy array + index_grid, // Flattened 2-D grid holding indices into nuclide grid for each unionized energy level + nuclide_grid, // Flattened 2-D grid holding energy levels and XS_data for all nuclides in simulation + SD.mats, // Flattened 2-D array with nuclide indices defining composition of each type of material + macro_xs_vector, // 1-D array with result of the macroscopic cross section (5 different reaction channels) + in.grid_type, // Lookup type (nuclide, hash, or unionized) + in.hash_bins, // Number of hash bins used (if using hash lookup type) + SD.max_num_nucs // Maximum number of nuclides present in any material + ); + + // For verification, and to prevent the compiler from optimizing + // all work out, we interrogate the returned macro_xs_vector array + // to find its maximum value index, then increment the verification + // value by that index. In this implementation, we prevent thread + // contention by using an OMP reduction on the verification value. + // For accelerators, a different approach might be required + // (e.g., atomics, reduction of thread-specific values in large + // array via CUDA thrust, etc). + double max = -1.0; + int max_idx = 0; + for(int j = 0; j < 5; j++ ) + { + if( macro_xs_vector[j] > max ) + { + max = macro_xs_vector[j]; + max_idx = j; + } + } + verification += max_idx+1; + } + + return verification; +} + +unsigned long long run_history_based_simulation(Inputs in, SimulationData SD, int mype) +{ + if( mype == 0) + printf("Beginning history based simulation...\n"); + + + //////////////////////////////////////////////////////////////////////////////// + // SUMMARY: Simulation Data Structure Manifest for "SD" Object + // Here we list all heap arrays (and lengths) in SD that would need to be + // offloaded manually if using an accelerator with a seperate memory space + //////////////////////////////////////////////////////////////////////////////// + // int * num_nucs; // Length = length_num_nucs; + // double * concs; // Length = length_concs + // int * mats; // Length = length_mats + // double * unionized_energy_array; // Length = length_unionized_energy_array + // int * index_grid; // Length = length_index_grid + // NuclideGridPoint * nuclide_grid; // Length = length_nuclide_grid + // + // Note: "unionized_energy_array" and "index_grid" can be of zero length + // depending on lookup method. + // + // Note: "Lengths" are given as the number of objects in the array, not the + // number of bytes. + //////////////////////////////////////////////////////////////////////////////// + + unsigned long long verification = 0; + + // Begin outer lookup loop over particles. This loop is independent. + int p = 0; + #pragma omp parallel for schedule(dynamic, 100) reduction(+:verification) + for( p = 0; p < in.particles; p++ ) + { + #ifdef AML + int * num_nucs = aml_replicaset_hwloc_local_replica(SD.num_nucs_replica); + double * concs = aml_replicaset_hwloc_local_replica(SD.concs_replica); + double * unionized_energy_array = aml_replicaset_hwloc_local_replica(SD.unionized_energy_array_replica); + int * index_grid = aml_replicaset_hwloc_local_replica(SD.index_grid_replica); + NuclideGridPoint * nuclide_grid = aml_replicaset_hwloc_local_replica(SD.nuclide_grid_replica); + #else + int * num_nucs = SD.num_nucs; + double * concs = SD.concs; + double * unionized_energy_array = SD.unionized_energy_array; + int * index_grid = SD.index_grid; + NuclideGridPoint * nuclide_grid = SD.nuclide_grid; + #endif + + // Set the initial seed value + uint64_t seed = STARTING_SEED; + + // Forward seed to lookup index (we need 2 samples per lookup, and + // we may fast forward up to 5 times after each lookup) + seed = fast_forward_LCG(seed, p*in.lookups*2*5); + + // Randomly pick an energy and material for the particle + double p_energy = LCG_random_double(&seed); + int mat = pick_mat(&seed); + + // Inner XS Lookup Loop + // This loop is dependent! + // i.e., Next iteration uses data computed in previous iter. + for( int i = 0; i < in.lookups; i++ ) + { + double macro_xs_vector[5] = {0}; + + // Perform macroscopic Cross Section Lookup + calculate_macro_xs( + p_energy, // Sampled neutron energy (in lethargy) + mat, // Sampled material type neutron is in + in.n_isotopes, // Total number of isotopes in simulation + in.n_gridpoints, // Number of gridpoints per isotope in simulation + num_nucs, // 1-D array with number of nuclides per material + concs, // Flattened 2-D array with concentration of each nuclide in each material + unionized_energy_array, // 1-D Unionized energy array + index_grid, // Flattened 2-D grid holding indices into nuclide grid for each unionized energy level + nuclide_grid, // Flattened 2-D grid holding energy levels and XS_data for all nuclides in simulation + SD.mats, // Flattened 2-D array with nuclide indices for each type of material + macro_xs_vector, // 1-D array with result of the macroscopic cross section (5 different reaction channels) + in.grid_type, // Lookup type (nuclide, hash, or unionized) + in.hash_bins, // Number of hash bins used (if using hash lookups) + SD.max_num_nucs // Maximum number of nuclides present in any material + ); + + + // For verification, and to prevent the compiler from optimizing + // all work out, we interrogate the returned macro_xs_vector array + // to find its maximum value index, then increment the verification + // value by that index. In this implementation, we prevent thread + // contention by using an OMP reduction on it. For other accelerators, + // a different approach might be required (e.g., atomics, reduction + // of thread-specific values in large array via CUDA thrust, etc) + double max = -1.0; + int max_idx = 0; + for(int j = 0; j < 5; j++ ) + { + if( macro_xs_vector[j] > max ) + { + max = macro_xs_vector[j]; + max_idx = j; + } + } + verification += max_idx+1; + + // Randomly pick next energy and material for the particle + // Also incorporates results from macro_xs lookup to + // enforce loop dependency. + // In a real MC app, this dependency is expressed in terms + // of branching physics sampling, whereas here we are just + // artificially enforcing this dependence based on fast + // forwarding the LCG state + uint64_t n_forward = 0; + for( int j = 0; j < 5; j++ ) + if( macro_xs_vector[j] > 1.0 ) + n_forward++; + if( n_forward > 0 ) + seed = fast_forward_LCG(seed, n_forward); + + p_energy = LCG_random_double(&seed); + mat = pick_mat(&seed); + } + + } + return verification; +} + +// Calculates the microscopic cross section for a given nuclide & energy +void calculate_micro_xs( double p_energy, int nuc, long n_isotopes, + long n_gridpoints, + double * restrict egrid, int * restrict index_data, + NuclideGridPoint * restrict nuclide_grids, + long idx, double * restrict xs_vector, int grid_type, int hash_bins ){ + // Variables + double f; + NuclideGridPoint * low, * high; + + // If using only the nuclide grid, we must perform a binary search + // to find the energy location in this particular nuclide's grid. + if( grid_type == NUCLIDE ) + { + // Perform binary search on the Nuclide Grid to find the index + idx = grid_search_nuclide( n_gridpoints, p_energy, &nuclide_grids[nuc*n_gridpoints], 0, n_gridpoints-1); + + // pull ptr from nuclide grid and check to ensure that + // we're not reading off the end of the nuclide's grid + if( idx == n_gridpoints - 1 ) + low = &nuclide_grids[nuc*n_gridpoints + idx - 1]; + else + low = &nuclide_grids[nuc*n_gridpoints + idx]; + } + else if( grid_type == UNIONIZED) // Unionized Energy Grid - we already know the index, no binary search needed. + { + // pull ptr from energy grid and check to ensure that + // we're not reading off the end of the nuclide's grid + if( index_data[idx * n_isotopes + nuc] == n_gridpoints - 1 ) + low = &nuclide_grids[nuc*n_gridpoints + index_data[idx * n_isotopes + nuc] - 1]; + else + low = &nuclide_grids[nuc*n_gridpoints + index_data[idx * n_isotopes + nuc]]; + } + else // Hash grid + { + // load lower bounding index + int u_low = index_data[idx * n_isotopes + nuc]; + + // Determine higher bounding index + int u_high; + if( idx == hash_bins - 1 ) + u_high = n_gridpoints - 1; + else + u_high = index_data[(idx+1)*n_isotopes + nuc] + 1; + + // Check edge cases to make sure energy is actually between these + // Then, if things look good, search for gridpoint in the nuclide grid + // within the lower and higher limits we've calculated. + double e_low = nuclide_grids[nuc*n_gridpoints + u_low].energy; + double e_high = nuclide_grids[nuc*n_gridpoints + u_high].energy; + int lower; + if( p_energy <= e_low ) + lower = 0; + else if( p_energy >= e_high ) + lower = n_gridpoints - 1; + else + lower = grid_search_nuclide( n_gridpoints, p_energy, &nuclide_grids[nuc*n_gridpoints], u_low, u_high); + + if( lower == n_gridpoints - 1 ) + low = &nuclide_grids[nuc*n_gridpoints + lower - 1]; + else + low = &nuclide_grids[nuc*n_gridpoints + lower]; + } + + high = low + 1; + + // calculate the re-useable interpolation factor + f = (high->energy - p_energy) / (high->energy - low->energy); + + // Total XS + xs_vector[0] = high->total_xs - f * (high->total_xs - low->total_xs); + + // Elastic XS + xs_vector[1] = high->elastic_xs - f * (high->elastic_xs - low->elastic_xs); + + // Absorbtion XS + xs_vector[2] = high->absorbtion_xs - f * (high->absorbtion_xs - low->absorbtion_xs); + + // Fission XS + xs_vector[3] = high->fission_xs - f * (high->fission_xs - low->fission_xs); + + // Nu Fission XS + xs_vector[4] = high->nu_fission_xs - f * (high->nu_fission_xs - low->nu_fission_xs); +} + +// Calculates macroscopic cross section based on a given material & energy +void calculate_macro_xs( double p_energy, int mat, long n_isotopes, + long n_gridpoints, int * restrict num_nucs, + double * restrict concs, + double * restrict egrid, int * restrict index_data, + NuclideGridPoint * restrict nuclide_grids, + int * restrict mats, + double * restrict macro_xs_vector, int grid_type, int hash_bins, int max_num_nucs ){ + int p_nuc; // the nuclide we are looking up + long idx = -1; + double conc; // the concentration of the nuclide in the material + + // cleans out macro_xs_vector + for( int k = 0; k < 5; k++ ) + macro_xs_vector[k] = 0; + + // If we are using the unionized energy grid (UEG), we only + // need to perform 1 binary search per macroscopic lookup. + // If we are using the nuclide grid search, it will have to be + // done inside of the "calculate_micro_xs" function for each different + // nuclide in the material. + if( grid_type == UNIONIZED ) + idx = grid_search( n_isotopes * n_gridpoints, p_energy, egrid); + else if( grid_type == HASH ) + { + double du = 1.0 / hash_bins; + idx = p_energy / du; + } + + // Once we find the pointer array on the UEG, we can pull the data + // from the respective nuclide grids, as well as the nuclide + // concentration data for the material + // Each nuclide from the material needs to have its micro-XS array + // looked up & interpolatied (via calculate_micro_xs). Then, the + // micro XS is multiplied by the concentration of that nuclide + // in the material, and added to the total macro XS array. + // (Independent -- though if parallelizing, must use atomic operations + // or otherwise control access to the xs_vector and macro_xs_vector to + // avoid simulataneous writing to the same data structure) + for( int j = 0; j < num_nucs[mat]; j++ ) + { + double xs_vector[5]; + p_nuc = mats[mat*max_num_nucs + j]; + conc = concs[mat*max_num_nucs + j]; + calculate_micro_xs( p_energy, p_nuc, n_isotopes, + n_gridpoints, egrid, index_data, + nuclide_grids, idx, xs_vector, grid_type, hash_bins ); + for( int k = 0; k < 5; k++ ) + macro_xs_vector[k] += xs_vector[k] * conc; + } +} + + +// binary search for energy on unionized energy grid +// returns lower index +long grid_search( long n, double quarry, double * restrict A) +{ + long lowerLimit = 0; + long upperLimit = n-1; + long examinationPoint; + long length = upperLimit - lowerLimit; + + while( length > 1 ) + { + examinationPoint = lowerLimit + ( length / 2 ); + + if( A[examinationPoint] > quarry ) + upperLimit = examinationPoint; + else + lowerLimit = examinationPoint; + + length = upperLimit - lowerLimit; + } + + return lowerLimit; +} + +// binary search for energy on nuclide energy grid +long grid_search_nuclide( long n, double quarry, NuclideGridPoint * A, long low, long high) +{ + long lowerLimit = low; + long upperLimit = high; + long examinationPoint; + long length = upperLimit - lowerLimit; + + while( length > 1 ) + { + examinationPoint = lowerLimit + ( length / 2 ); + + if( A[examinationPoint].energy > quarry ) + upperLimit = examinationPoint; + else + lowerLimit = examinationPoint; + + length = upperLimit - lowerLimit; + } + + return lowerLimit; +} + +// picks a material based on a probabilistic distribution +int pick_mat( uint64_t * seed ) +{ + // I have a nice spreadsheet supporting these numbers. They are + // the fractions (by volume) of material in the core. Not a + // *perfect* approximation of where XS lookups are going to occur, + // but this will do a good job of biasing the system nonetheless. + + double dist[12]; + dist[0] = 0.140; // fuel + dist[1] = 0.052; // cladding + dist[2] = 0.275; // cold, borated water + dist[3] = 0.134; // hot, borated water + dist[4] = 0.154; // RPV + dist[5] = 0.064; // Lower, radial reflector + dist[6] = 0.066; // Upper reflector / top plate + dist[7] = 0.055; // bottom plate + dist[8] = 0.008; // bottom nozzle + dist[9] = 0.015; // top nozzle + dist[10] = 0.025; // top of fuel assemblies + dist[11] = 0.013; // bottom of fuel assemblies + + double roll = LCG_random_double(seed); + + // makes a pick based on the distro + for( int i = 0; i < 12; i++ ) + { + double running = 0; + for( int j = i; j > 0; j-- ) + running += dist[j]; + if( roll < running ) + return i; + } + + return 0; +} + +double LCG_random_double(uint64_t * seed) +{ + // LCG parameters + const uint64_t m = 9223372036854775808ULL; // 2^63 + const uint64_t a = 2806196910506780709ULL; + const uint64_t c = 1ULL; + *seed = (a * (*seed) + c) % m; + return (double) (*seed) / (double) m; +} + +uint64_t fast_forward_LCG(uint64_t seed, uint64_t n) +{ + // LCG parameters + const uint64_t m = 9223372036854775808ULL; // 2^63 + uint64_t a = 2806196910506780709ULL; + uint64_t c = 1ULL; + + n = n % m; + + uint64_t a_new = 1; + uint64_t c_new = 0; + + while(n > 0) + { + if(n & 1) + { + a_new *= a; + c_new = c_new * a + c; + } + c *= (a + 1); + a *= a; + + n >>= 1; + } + + return (a_new * seed + c_new) % m; + +} + +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// +// OPTIMIZED VARIANT FUNCTIONS +//////////////////////////////////////////////////////////////////////////////////// +// This section contains a number of optimized variants of some of the above +// functions, which each deploy a different combination of optimizations strategies. +// By default, XSBench will not run any of these variants. They +// must be specifically selected using the "-k " command +// line argument. +// +// As fast parallel sorting will be required for these optimizations, we will +// first define a set of key-value parallel quicksort routines. +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////// + + +//////////////////////////////////////////////////////////////////////////////////// +// Parallel Quicksort Key-Value Sorting Algorithms +//////////////////////////////////////////////////////////////////////////////////// +// +// These algorithms are based on the parallel quicksort implementation by +// Eduard Lopez published at https://github.com/eduardlopez/quicksort-parallel +// +// Eduard's original version was for an integer type quicksort, but I have modified +// it to form two different versions that can sort key-value pairs together without +// having to bundle them into a separate object. Additionally, I have modified the +// optimal chunk sizes and restricted the number of threads for the array sizing +// that XSBench will be using by default. +// +// Eduard's original implementation carries the following license, which applies to +// the following functions only: +// +// void quickSort_parallel_internal_i_d(int* key,double * value, int left, int right, int cutoff) +// void quickSort_parallel_i_d(int* key,double * value, int lenArray, int numThreads) +// void quickSort_parallel_internal_d_i(double* key,int * value, int left, int right, int cutoff) +// void quickSort_parallel_d_i(double* key,int * value, int lenArray, int numThreads) +// +// The MIT License (MIT) +// +// Copyright (c) 2016 Eduard López +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//////////////////////////////////////////////////////////////////////////////////// +void quickSort_parallel_internal_i_d(int* key,double * value, int left, int right, int cutoff) +{ + int i = left, j = right; + int tmp; + int pivot = key[(left + right) / 2]; + + { + while (i <= j) { + while (key[i] < pivot) + i++; + while (key[j] > pivot) + j--; + if (i <= j) { + tmp = key[i]; + key[i] = key[j]; + key[j] = tmp; + double tmp_v = value[i]; + value[i] = value[j]; + value[j] = tmp_v; + i++; + j--; + } + } + + } + + if ( ((right-left) 16 ) + numThreads = 16; + + #pragma omp parallel num_threads(numThreads) + { + #pragma omp single nowait + { + quickSort_parallel_internal_i_d(key,value, 0, lenArray-1, cutoff); + } + } + +} + +void quickSort_parallel_internal_d_i(double* key,int * value, int left, int right, int cutoff) +{ + int i = left, j = right; + double tmp; + double pivot = key[(left + right) / 2]; + + { + while (i <= j) { + while (key[i] < pivot) + i++; + while (key[j] > pivot) + j--; + if (i <= j) { + tmp = key[i]; + key[i] = key[j]; + key[j] = tmp; + int tmp_v = value[i]; + value[i] = value[j]; + value[j] = tmp_v; + i++; + j--; + } + } + + } + + if ( ((right-left) 16 ) + numThreads = 16; + + #pragma omp parallel num_threads(numThreads) + { + #pragma omp single nowait + { + quickSort_parallel_internal_d_i(key,value, 0, lenArray-1, cutoff); + } + } + +} + +//////////////////////////////////////////////////////////////////////////////////// +// Optimization 1 -- Event-based Sample/XS Lookup kernel splitting + Sorting +// lookups by material and energy +//////////////////////////////////////////////////////////////////////////////////// +// This kernel separates out the sampling and lookup regions of the event-based +// model, and then sorts the lookups by material type and energy. The goal of this +// optimization is to allow for greatly improved cache locality, and XS indices +// loaded from memory may be re-used for multiple lookups. +// +// As efficienct sorting is key for performance, we also must implement an +// efficient key-value parallel sorting algorithm. We also experimented with using +// the C++ version of thrust for these purposes, but found that our own implemtation +// was slightly faster than the thrust library version, so for speed and +// simplicity we will do not add the thrust dependency. +//////////////////////////////////////////////////////////////////////////////////// + + +unsigned long long run_event_based_simulation_optimization_1(Inputs in, SimulationData SD, int mype) +{ + char * optimization_name = "Optimization 1 - Kernel splitting + full material & energy sort"; + + if( mype == 0) printf("Simulation Kernel:\"%s\"\n", optimization_name); + + //////////////////////////////////////////////////////////////////////////////// + // Allocate Additional Data Structures Needed by Optimized Kernel + //////////////////////////////////////////////////////////////////////////////// + if( mype == 0) printf("Allocating additional data required by optimized kernel...\n"); + size_t sz; + size_t total_sz = 0; + double start, stop; + + // loop variables + int i = 0; + int m = 0; + + sz = in.lookups * sizeof(double); + SD.p_energy_samples = (double *) malloc(sz); + total_sz += sz; + SD.length_p_energy_samples = in.lookups; + + sz = in.lookups * sizeof(int); + SD.mat_samples = (int *) malloc(sz); + total_sz += sz; + SD.length_mat_samples = in.lookups; + + if( mype == 0) printf("Allocated an additional %.0lf MB of data on GPU.\n", total_sz/1024.0/1024.0); + + //////////////////////////////////////////////////////////////////////////////// + // Begin Actual Simulation + //////////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////////////////////////////////////////// + // Sample Materials and Energies + //////////////////////////////////////////////////////////////////////////////// + #pragma omp parallel for schedule(dynamic, 100) + for( i = 0; i < in.lookups; i++ ) + { + // Set the initial seed value + uint64_t seed = STARTING_SEED; + + // Forward seed to lookup index (we need 2 samples per lookup) + seed = fast_forward_LCG(seed, 2*i); + + // Randomly pick an energy and material for the particle + double p_energy = LCG_random_double(&seed); + int mat = pick_mat(&seed); + + SD.p_energy_samples[i] = p_energy; + SD.mat_samples[i] = mat; + } + if(mype == 0) printf("finished sampling...\n"); + + //////////////////////////////////////////////////////////////////////////////// + // Sort by Material + //////////////////////////////////////////////////////////////////////////////// + + start = get_time(); + + quickSort_parallel_i_d(SD.mat_samples, SD.p_energy_samples, in.lookups, in.nthreads); + + stop = get_time(); + + if(mype == 0) printf("Material sort took %.3lf seconds\n", stop-start); + + //////////////////////////////////////////////////////////////////////////////// + // Sort by Energy + //////////////////////////////////////////////////////////////////////////////// + + start = get_time(); + + // Count up number of each type of sample. + int num_samples_per_mat[12] = {0}; + for( int l = 0; l < in.lookups; l++ ) + num_samples_per_mat[ SD.mat_samples[l] ]++; + + // Determine offsets + int offsets[12] = {0}; + for( int m = 1; m < 12; m++ ) + offsets[m] = offsets[m-1] + num_samples_per_mat[m-1]; + + stop = get_time(); + if(mype == 0) printf("Counting samples and offsets took %.3lf seconds\n", stop-start); + start = stop; + + // Sort each material type by energy level + int offset = 0; + for( int m = 0; m < 12; m++ ) + quickSort_parallel_d_i(SD.p_energy_samples + offsets[m],SD.mat_samples + offsets[m], num_samples_per_mat[m], in.nthreads); + + stop = get_time(); + if(mype == 0) printf("Energy Sorts took %.3lf seconds\n", stop-start); + + //////////////////////////////////////////////////////////////////////////////// + // Perform lookups for each material separately + //////////////////////////////////////////////////////////////////////////////// + start = get_time(); + + unsigned long long verification = 0; + + // Individual Materials + offset = 0; + for( m = 0; m < 12; m++ ) + { + #pragma omp parallel for schedule(dynamic,100) reduction(+:verification) + for( i = offset; i < offset + num_samples_per_mat[m]; i++) + { + #ifdef AML + int * num_nucs = aml_replicaset_hwloc_local_replica(SD.num_nucs_replica); + double * concs = aml_replicaset_hwloc_local_replica(SD.concs_replica); + double * unionized_energy_array = aml_replicaset_hwloc_local_replica(SD.unionized_energy_array_replica); + int * index_grid = aml_replicaset_hwloc_local_replica(SD.index_grid_replica); + NuclideGridPoint * nuclide_grid = aml_replicaset_hwloc_local_replica(SD.nuclide_grid_replica); + #else + int * num_nucs = SD.num_nucs; + double * concs = SD.concs; + double * unionized_energy_array = SD.unionized_energy_array; + int * index_grid = SD.index_grid; + NuclideGridPoint * nuclide_grid = SD.nuclide_grid; + #endif + + // load pre-sampled energy and material for the particle + double p_energy = SD.p_energy_samples[i]; + int mat = SD.mat_samples[i]; + + double macro_xs_vector[5] = {0}; + + // Perform macroscopic Cross Section Lookup + calculate_macro_xs( + p_energy, // Sampled neutron energy (in lethargy) + mat, // Sampled material type index neutron is in + in.n_isotopes, // Total number of isotopes in simulation + in.n_gridpoints, // Number of gridpoints per isotope in simulation + num_nucs, // 1-D array with number of nuclides per material + concs, // Flattened 2-D array with concentration of each nuclide in each material + unionized_energy_array, // 1-D Unionized energy array + index_grid, // Flattened 2-D grid holding indices into nuclide grid for each unionized energy level + nuclide_grid, // Flattened 2-D grid holding energy levels and XS_data for all nuclides in simulation + SD.mats, // Flattened 2-D array with nuclide indices defining composition of each type of material + macro_xs_vector, // 1-D array with result of the macroscopic cross section (5 different reaction channels) + in.grid_type, // Lookup type (nuclide, hash, or unionized) + in.hash_bins, // Number of hash bins used (if using hash lookup type) + SD.max_num_nucs // Maximum number of nuclides present in any material + ); + + // For verification, and to prevent the compiler from optimizing + // all work out, we interrogate the returned macro_xs_vector array + // to find its maximum value index, then increment the verification + // value by that index. In this implementation, we prevent thread + // contention by using an OMP reduction on the verification value. + // For accelerators, a different approach might be required + // (e.g., atomics, reduction of thread-specific values in large + // array via CUDA thrust, etc). + double max = -1.0; + int max_idx = 0; + for(int j = 0; j < 5; j++ ) + { + if( macro_xs_vector[j] > max ) + { + max = macro_xs_vector[j]; + max_idx = j; + } + } + verification += max_idx+1; + } + offset += num_samples_per_mat[m]; + } + + stop = get_time(); + if(mype == 0) printf("XS Lookups took %.3lf seconds\n", stop-start); + return verification; +} \ No newline at end of file diff --git a/benchmarks/benchmarks/XSbench/XSbench_header.h b/benchmarks/benchmarks/XSbench/XSbench_header.h index 53157a6..a0229ed 100644 --- a/benchmarks/benchmarks/XSbench/XSbench_header.h +++ b/benchmarks/benchmarks/XSbench/XSbench_header.h @@ -148,4 +148,4 @@ double get_time(void); int * load_num_nucs(long n_isotopes); int * load_mats( int * num_nucs, long n_isotopes, int * max_num_nucs ); double * load_concs( int * num_nucs, int max_num_nucs ); -#endif +#endif \ No newline at end of file diff --git a/benchmarks/benchmarks/XSbench/XSutils.c b/benchmarks/benchmarks/XSbench/XSutils.c index 9e54220..72984c1 100644 --- a/benchmarks/benchmarks/XSbench/XSutils.c +++ b/benchmarks/benchmarks/XSbench/XSutils.c @@ -60,4 +60,4 @@ double get_time(void) double time = (double) ms / 1000.0; return time; -} +} \ No newline at end of file diff --git a/benchmarks/benchmarks/XSbench/io.c b/benchmarks/benchmarks/XSbench/io.c index 2f60188..f8b85c3 100644 --- a/benchmarks/benchmarks/XSbench/io.c +++ b/benchmarks/benchmarks/XSbench/io.c @@ -505,4 +505,4 @@ SimulationData binary_read( Inputs in ) fclose(fp); return SD; -} +} \ No newline at end of file diff --git a/benchmarks/benchmarks/kmeans/bitmap.h b/benchmarks/benchmarks/kmeans/bitmap.h new file mode 100644 index 0000000..30a3900 --- /dev/null +++ b/benchmarks/benchmarks/kmeans/bitmap.h @@ -0,0 +1,395 @@ +/* + * Copyright (c) 2011 Bharath Ramesh + * + * Distributed under the terms of GNU LGPL, version 2.1 + */ + +#include +#include +// #include "brmalloc.h" + +// ------------ +// User defined +#define MIN_ALLOC_SIZE ((size_t)8U) +#define MAX_ALLOC_SIZE ((size_t)4096U) +#define MAX_NUM_ALLOCS 8192 +#define NUM_ITERATIONS 32UNIMPLEMENTED + +// ------------ + +#define BITS_PER_LONG 64 +#define LOG_BITS_PER_BYTE 3 +#define LOG_BITS_PER_LONG 6 +#define LOG_BYTES_PER_LONG 3 +#define MAX_LONG 0xffffffffffffffff +#define K ((size_t)1024U) +#define M ((size_t)1024U * K) +#define G ((size_t)1024U * M) +#define MALLOC_OVERHEAD ((size_t)16U) + +#ifndef ABORT +#define ABORT() abort() +#endif + +#ifndef EXIT_ERR +#define EXIT_ERR() exit(EXIT_FAILURE) +#endif + +#ifndef DEBUG +#include +#define DEBUG(fmt, args...) printf("%d, %s: " fmt, __LINE__, __FUNCTION__, ##args) +#endif + +#ifndef MALLOC_QUANTA +#define MALLOC_QUANTA ((size_t)16U) +#endif + +#ifndef MALLOC_THRESHOLD +#define MALLOC_THRESHOLD ((size_t)1U * M) +#endif + +#ifndef MALLOC_ZONE_SIZE +#define MALLOC_ZONE_SIZE ((size_t)16U * M) +#endif + +#ifndef MMAP +#include +#define MMAP_FLAGS (MAP_ANONYMOUS | MAP_PRIVATE) +#define MMAP_PROT (PROT_READ | PROT_WRITE) +#define MMAP(s) mmap (0, (s), MMAP_PROT, MMAP_FLAGS, -1, 0) +#define MUNMAP(a, s) munmap ((a), (s)) +#endif + +struct __malloc_zone_type { + unsigned char *bmap; + size_t bmap_longs; + size_t free_size; + uint64_t id; + void *region; + uint64_t start_byte; + struct __malloc_zone_type *next; +}; + +typedef struct __malloc_zone_type malloc_zone_t; + +static malloc_zone_t *mz_head = NULL, *mz_tail = NULL; +static uint8_t mqbits = 0; + +static inline void +free_region (malloc_zone_t *mz, uint64_t loc, size_t nbytes) +{ + unsigned char *bmap; + size_t bmap_longs; + uint64_t i, i_b, j_b, lbits, nbits, pattern, *wordptr; + + lbits = nbits = nbytes >> mqbits; + i_b = loc >> LOG_BITS_PER_LONG; + j_b = loc - (i_b << LOG_BITS_PER_LONG); + bmap = mz->bmap; + bmap_longs = mz->bmap_longs; + for (i = i_b; i < bmap_longs; ++i) { + if (lbits == 0) + break; + + wordptr = (uint64_t *) bmap + i; + if (i == i_b) { + if ((j_b + nbits) < BITS_PER_LONG) { + pattern = ((uint64_t) 1 << nbits) - 1; + pattern = ~(pattern << j_b); + *wordptr &= pattern; + return; + } + + pattern = ~(MAX_LONG << j_b); + *wordptr &= pattern; + lbits -= (BITS_PER_LONG -j_b); + continue; + } + + if (lbits >= BITS_PER_LONG) { + *wordptr = 0; + lbits -= BITS_PER_LONG; + continue; + } + + pattern = ~(((uint64_t) 1 << lbits) - 1); + *wordptr &= pattern; + return; + } + + return; +} + +static inline uint64_t +get_region (malloc_zone_t *mz, size_t nbytes) +{ + unsigned char *bmap; + size_t bmap_longs; + uint64_t i = 0, i_b, j, j_b, l, nbits, rbits, pattern, word, *wordptr; + int64_t loc = -1; + + nbits = rbits = nbytes >> mqbits; + bmap = mz->bmap; + bmap_longs = mz->bmap_longs; + l = mz->start_byte; + while (i < bmap_longs) { + if (rbits == 0) + break; + + word = *((uint64_t *) bmap + l); + if (word == MAX_LONG) { + loc = -1; + rbits = nbits; + ++i; + ++l; + if (l == bmap_longs) + l = 0; + + continue; + } + + if (word == 0) { + if (rbits >= BITS_PER_LONG) { + rbits -= BITS_PER_LONG; + if (loc == -1) + loc = l << LOG_BITS_PER_LONG; + + ++i; + ++l; + if (l == bmap_longs) { + l = 0; + loc = -1; + rbits = nbits; + } + + continue; + } + + rbits = 0; + if (loc == -1) + loc = l << LOG_BITS_PER_LONG; + + break; + } + + if (rbits >= BITS_PER_LONG) { + loc = -1; + rbits = nbits; + ++i; + ++l; + if (l == bmap_longs) + l = 0; + + continue; + } + + for (j = 0; j < BITS_PER_LONG; ++j) { + if (rbits == 0) + break; + + if ((word >> j) & 1) { + loc = -1; + rbits = nbits; + ++i; + ++l; + if (l == bmap_longs) + l = 0; + + break; + } + + --rbits; + if (loc == -1) + loc = (l << LOG_BITS_PER_LONG) + j; + } + } + + if ((rbits != 0) || (loc == -1)) + return -1; + + mz->start_byte = (loc + nbits) >> LOG_BITS_PER_LONG; + i_b = loc >> LOG_BITS_PER_LONG; + j_b = loc - (i_b << LOG_BITS_PER_LONG); + rbits = nbits; + for (i = i_b; i < bmap_longs; ++i) { + if (rbits == 0) + break; + + wordptr = (uint64_t *) bmap + i; + if (i == i_b) { + if ((j_b + nbits) < BITS_PER_LONG) { + pattern = ((uint64_t) 1 << nbits) - 1; + pattern = pattern << j_b; + *wordptr |= pattern; + return loc; + } + + pattern = MAX_LONG << j_b; + *wordptr |= pattern; + rbits -= (BITS_PER_LONG - j_b); + continue; + } + + if (rbits >= BITS_PER_LONG) { + *wordptr = MAX_LONG; + rbits -= BITS_PER_LONG; + continue; + } + + pattern = ((uint64_t) 1 << rbits) - 1; + *wordptr |= pattern; + return loc; + } + + + return loc; +} + +static inline malloc_zone_t * +new_malloc_zone (void) +{ + size_t bmap_size, bmap_longs; + malloc_zone_t *mz; + void *region; + + region = MMAP (MALLOC_ZONE_SIZE); + if (region == MAP_FAILED) { + DEBUG ("ERROR: MMAP failed.\n"); + return NULL; + } + + bmap_size = (MALLOC_ZONE_SIZE / MALLOC_QUANTA) >> LOG_BITS_PER_BYTE; + bmap_longs = bmap_size >> LOG_BYTES_PER_LONG; + mz = (malloc_zone_t *) malloc (sizeof (malloc_zone_t)); + mz->bmap = (unsigned char *) malloc (bmap_size); + if (mz->bmap == NULL) { + DEBUG ("ERROR: Unable to allocate bitmap.\n"); + free (mz); + MUNMAP (region, MALLOC_ZONE_SIZE); + return NULL; + } + + memset (mz->bmap, 0, bmap_size); + mz->bmap_longs = bmap_longs; + mz->free_size = MALLOC_ZONE_SIZE; + mz->id = (uint64_t) mz; + mz->region = region; + mz->start_byte = 0; + mz->next = NULL; + if (mz_tail != NULL) + mz_tail->next = mz; + + mz_tail = mz; + if (mz_head == NULL) + mz_head = mz; + + return mz; +} + +static void __attribute__ ((constructor)) +brm_init (void) +{ + uint8_t i, no_ones; + size_t size; + + no_ones = 0; + size = MALLOC_QUANTA; + for (i = 0; i < BITS_PER_LONG; ++i) { + if ((size >> i) & 1) { + ++no_ones; + if (no_ones > 1) { + DEBUG ("ERROR: MALLOC_QUANTA not power of " + "2.\n"); + EXIT_ERR (); + } + + mqbits = i; + } + } + + if (no_ones == 0) { + DEBUG ("ERROR: MALLOC_QUANTA set to 0 (zero).\n"); + EXIT_ERR (); + } + + if (new_malloc_zone () == NULL) { + DEBUG ("ERROR: new_malloc_zone failed.\n"); + EXIT_ERR (); + } + + return; +} + +void +brm_free (void *ptr) +{ + uint64_t *base; + malloc_zone_t *mz; + int64_t offset; + size_t size; + + base = (uint64_t *) ptr - 2; + mz = (malloc_zone_t *) *((uint64_t *) base); + size = (size_t) *((uint64_t *) base + 1); + if (mz->id != (uint64_t) mz) { + DEBUG ("ptr: %p, base: %p\n", ptr, base); + DEBUG ("ERROR: data corruption.\n"); + ABORT (); + } + + offset = ((void *) base - mz->region) >> mqbits; + free_region (mz, offset, size); + mz->free_size += size; + + return; +} + +void * +brm_malloc (size_t size) +{ + uint64_t addr; + malloc_zone_t *mz; + size_t nbytes; + int64_t offset; + + if (size == 0) + return NULL; + + nbytes = (size + MALLOC_OVERHEAD + MALLOC_QUANTA - 1) & + ~(MALLOC_QUANTA - 1); + // printf(nbytes); + // if (nbytes >= MALLOC_THRESHOLD) { + // DEBUG ("UNIMPLEMENTED\n"); + // EXIT_ERR(); + // } + + mz = mz_head; + while (1) { + if (nbytes > mz->free_size) { + DEBUG ("UNIMPLEMENTED\n"); + EXIT_ERR (); + } + + offset = get_region (mz, nbytes); + if (offset == -1) { + mz = mz->next; + if (mz == NULL) { + if ((mz = new_malloc_zone ()) == NULL) { + DEBUG ("ERROR: new_malloc zone " + "failed.\n"); + EXIT_ERR (); + } + } + + continue; + } + + mz->free_size -= nbytes; + addr = (uint64_t) mz->region + ((size_t) offset << mqbits); + *((uint64_t *) addr) = (uint64_t) mz; + *((uint64_t *) addr + 1) = (uint64_t) nbytes; + break; + } + + return (void *) (addr + 2 * sizeof (uint64_t)); +} \ No newline at end of file diff --git a/benchmarks/benchmarks/kmeans/kmeans-pthread.c b/benchmarks/benchmarks/kmeans/kmeans-pthread.c index 0c0f927..7d08df8 100644 --- a/benchmarks/benchmarks/kmeans/kmeans-pthread.c +++ b/benchmarks/benchmarks/kmeans/kmeans-pthread.c @@ -39,8 +39,8 @@ #include "coz.h" -// #define malloc MALLOCCHERI -// #define free FREECHERI +#define malloc MALLOC +#define free FREE #define DEF_NUM_POINTS 150000 #define DEF_NUM_MEANS 100 @@ -271,6 +271,7 @@ int main(int argc, char **argv) // printf("Initial alloc called\n"); //INITAlloc(); //INITREGULARALLOC(); + init_malloc(); int num_procs, curr_point; int i; diff --git a/benchmarks/benchmarks/kmeans/simple.h b/benchmarks/benchmarks/kmeans/simple.h new file mode 100644 index 0000000..9de2a96 --- /dev/null +++ b/benchmarks/benchmarks/kmeans/simple.h @@ -0,0 +1,177 @@ +/* Copyright (C) 2023. Shivashish Das. Licensed under the MIT License.*/ +#include +#include +#include +#include +#include + +// source: https://www.reddit.com/r/C_Programming/comments/1bt8dyz/github_dasshivamalloc_a_simple_memory_allocator/ + +// #include "alloc.h" +#ifndef _MSC_VER +#include +#endif +/* This is a simple memory allocator meant for use in single threaded applications. + First we get memory from the system allocator which is defined by the pool size + Larger the pool size, more is the amount you can allocate before running out of memory. + On linux, mmap() is used for memory allocation while on windows good old calloc() is used as the system allocator + + Some basic definitions: + Block - A memory region always of size 16 bytes. This is the basic unit of allocation + All allocations are made in multiples of blocks. If any allocation request is not a multiple of 16 bytes, we return memory of a size + that is the closest multiple to 16 and greater than the user requested size. + + Metadata blocks - For every allocation we allocate two extra blocks. These two blocks hold data about the allocation itself + and serve to prevent buffer overflows too. Check the comment in alloc() to find out more. For example suppose that if the user asks for 80 bytes + i.e 80 / 16 = 5 blocks, we will allocate 7 blocks but the pointer passed to the user will point to the seconf block so the user is unable to access + these blocks. They do sound like a waste of some bytes but help provide protection from buffer overflows + + Posioning - This means that user code has overflown the buffer it was allocated. All memory allocated by alloc() is now invalid + However this may also be caused by the user simply passing an invalid pointer to us i.e a pointer allocated by some other allocator etc. + In this case the user can clear the poisoned state by calling clear_posion() but be absolutely sure as a user about this before doing so +*/ +static uint8_t* mem = 0; +static uint8_t* bitmap = 0; +static uint64_t blocks = 0; +static uint8_t poison = 0; + +// Always call this before anything else. +void alloc_init(uint64_t pool) { + if (pool % 16 != 0) { + int rem = pool % 16; + pool += (16 - rem); + } + +#ifndef _MSC_VER + mem = mmap(0, pool, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANON, -1, 0); + if (mem == MAP_FAILED) { + fprintf(stderr, "alloc_init(): could not allocate heap\n"); + perror("mmap"); + return; + } +#else + mem = calloc(pool, 1); + if (!mem) { + fprintf(stderr, "alloc_init(): could not allocate heap\n"); + exit(1); + } +#endif + // Each bitmap entry can represent 8 blocks and each block is 16 bytes + // So space representable in one uint8_t is 16 * 8 = 128 bytes + uint64_t sz = pool / 128; + if (sz == 0) + sz = 1; // allocate at least one to keep track of small pools + +#ifndef _MSC_VER + bitmap = mmap(0, sz , PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANON, -1, 0); + if (bitmap == MAP_FAILED) { + fprintf(stderr, "alloc_init(): could not allocate bitmap"); + munmap(mem, pool); + } +#else + bitmap = calloc(sz, 1); + if (!bitmap) { + fprintf(stderr, "alloc_init(): could not allocate bitmap\n"); + exit(1); + } +#endif + // Zero the entire bitmap + memset(bitmap, 0, sz); + blocks = pool / 16; +} + +#define IS_FREE(blkid) (bitmap[blkid / 8] & (((uint8_t)1) << blkid)) == 0 +#define MARK(blkid) (bitmap[blkid / 8] ^= ((((uint8_t)1) << blkid) - 1)) + +// Allocate sz bytes of memory. Caution: May allocate upto 15 bytes more than sz +void* alloc(uint64_t sz) { + if (sz % 16 != 0) { + int rem = sz % 16; + sz += (16 - rem); + } + // Allocate two extra blocks + // First will be allocated at just behind the first user accessible block + // This block will have the number of blocks allocated and a randomly generated magic number each 8 bytes long + // The last block has the "magic number" present in the first block + // If this magic number gets modified then when free() tries to free the memory + // Buffer overruns will be caught and this allocator gets poisoned i.e it can no longer allocate memory + // This is because all blocks are laid out sequentially and if the user overruns the blocks allocated + // Then the user may have overwritten the contents of other blocks and it is not possible to estimate the damage caused + // and data corrupted. All pointers to blocks allocated immediately become invalid and free() posions the allocator + // This helps catch buffer overflows early on + uint64_t blk = (sz / 16) + 2; + + // if we are posioned, all allocation requests will fail + if (poison) + return 0; + // Loop through the entire bitmap. If a free block is found, check if there are at least blk free blocks after it. + // If such a contigious group of blocks is found, take appropriate actions and return to user + // Otherwise we have ran out of memory so inform the user about it + for (uint64_t i = 0; i < blocks; i++) { + if (IS_FREE(i)) { + // Check for contigious free blocks + for (uint64_t j = i; j < (i + blk); j++) { + if (!IS_FREE(j)) + goto next; + } + + // Mark all free blocks + for (uint64_t j = i; j < (i + blk + 1); j++) { + MARK(j); + } + uint64_t* ptr = mem + (i * 16); + *ptr = blk; + + // I needed a number which was large enough to occupy 8 bytes so rand() is not enough as in most cases RAND_MAX is only USHORT_MAX + // Instead use time() which returns a 64 bit value and is almost guaranteed to be unique on every call to alloc() + uint64_t magic = time(0); + *(ptr + 1) = magic; + + // Store a magic number in the last block. For the reason see free_mem() + ptr = mem + (i * 16) + ((blk - 1) * 16); + *ptr = magic; + *(ptr + 1) = magic; + // Return the user a pointer which points to the region just above our metadata block + return mem + ((i + 1) * 16); + } +next: + } + fprintf(stderr, "Pool has been exhausted...Cannot allocate more memory"); + return 0; +} + +// Frees memory allocated by alloc() +void free_mem(void* data) { + // First get the number of blocks allocated and magic from the metadata block (i.e the block right behind what alloc() returned) + uint64_t* ptr = data; + ptr -= 2; + uint64_t blk = *ptr; + ptr++; + uint64_t magic = *ptr; + + // The magic is stored in the last block of the allocation + // Compare the two magic values + // If they are equal, this memory block was allocated by us and we can free this + // Otherwise the buffer has been overflown which has overwritten the magic number or this was not allocated by alloc() and is not ours to deal with + ptr = data + (blk - 2) * 16; + if (magic != *ptr) { + // If the buffer has overflown then mark this allocator posioned. + // You may change the poison back to 0 in your code but be careful and do this only if you know that the buffer was not overrun + fprintf(stderr, "Invalid pointer or buffer overrun detected..Poisoning ourself"); + poison = 1; + return; + } + uint64_t offset = ((uint8_t*) data) - mem; + offset -= 16; + offset /= 16; + + // Clear all bits representing this block so next call to alloc() can use this + for (uint64_t j = offset; j < offset + blk + 1; j++) { + MARK(j); + } +} + +// Do not call this unless you are absolutely sure about the cause of poisoning +void clear_posion() { + poison = 0; +} \ No newline at end of file diff --git a/benchmarks/benchmarks/kmeans/stddefines.h b/benchmarks/benchmarks/kmeans/stddefines.h index 595bdcb..966bf69 100644 --- a/benchmarks/benchmarks/kmeans/stddefines.h +++ b/benchmarks/benchmarks/kmeans/stddefines.h @@ -48,9 +48,37 @@ #include #include +// #include "bitmap.h" +#include "simple.h" + #define MAXPAGESIZES 2 +init_malloc(void) { + // Init malloc implementation + //INITREGULARALLOC(); + + // init_alloc(100,100000); +// brm_init(); + alloc_init(104857632); +} + +void* MALLOC(size_t sz) { + // malloc implementation + //return MALLOCCHERI(sz); +// return malloc_buddy(sz); + return alloc(sz); + // return (void *)alloc_chunk(); +} + +void FREE(void *ptr) { + // free implementation + //FREECHERI(ptr); + // free_chunk(ptr); + free_mem(ptr); +} + + //#define TIMING /* Debug printf */ @@ -64,12 +92,12 @@ assert ((a) == 0); \ } -static inline void *MALLOC(size_t size) -{ - void * temp = malloc(size); - assert(temp); - return temp; -} +// static inline void *MALLOC(size_t size) +// { +// void * temp = malloc(size); +// assert(temp); +// return temp; +// } static inline void *CALLOC(size_t num, size_t size) { @@ -142,43 +170,43 @@ int MallocCounter; size_t sizeUsed; -INITAlloc(void) { +// INITAlloc(void) { - size_t sz; - // Pre Allocate 600 MB - sz = 100000000; +// size_t sz; +// // Pre Allocate 600 MB +// sz = 100000000; - int fd = open(FILENAME, O_RDWR, 0600); +// int fd = open(FILENAME, O_RDWR, 0600); - if (fd < 0) { - perror("open"); - exit(EXIT_FAILURE); - } +// if (fd < 0) { +// perror("open"); +// exit(EXIT_FAILURE); +// } - off_t offset = 0; // offset to seek to. +// off_t offset = 0; // offset to seek to. - if (ftruncate(fd, sz) < 0) { - perror("ftruncate"); - close(fd); - exit(EXIT_FAILURE); - } +// if (ftruncate(fd, sz) < 0) { +// perror("ftruncate"); +// close(fd); +// exit(EXIT_FAILURE); +// } - // ptr = mmap(NULL, sz, - // PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANON,-1,0); +// // ptr = mmap(NULL, sz, +// // PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANON,-1,0); - ptr = mmap(NULL, sz, - PROT_READ|PROT_WRITE, MAP_SHARED,fd,0); +// ptr = mmap(NULL, sz, +// PROT_READ|PROT_WRITE, MAP_SHARED,fd,0); - // Added error handling - if(ptr == MAP_FAILED) - { - perror("mmap"); - exit(EXIT_FAILURE); - } +// // Added error handling +// if(ptr == MAP_FAILED) +// { +// perror("mmap"); +// exit(EXIT_FAILURE); +// } - MallocCounter = (int)sz; +// MallocCounter = (int)sz; -} +// } // Quick malloc implementation with mmap void* MALLOCCHERI(size_t sz) @@ -285,14 +313,170 @@ INITREGULARALLOC(void) { MallocCounter = (int)sz; } -// Standard Alloc -// void* MALLOCREGULAR(size_t sz) { - -// } -// void* CLEARALLOC(void) { -// / -// } +// ------------------------------ bitmap allocator ---------------------------- -#endif // STDDEFINES_H_ +#define BITS_PER_BYTE 8 + +char *buffer = NULL; /* allocation buffer */ +unsigned char *bitmap = NULL; /* bitmap for the buffer */ + +int buffer_size = 0; /* size of buffer (in bytes) */ +int bitmap_size = 0; /* size of bitmap (in bytes) */ +int bytes_per_chunk = 0; /* size of single chunk (in bytes) */ + +void init_alloc(int num_chunks, int chunk_size) +{ + int i = 0; + + /* we need to increase the num_chunks + * so every bit in bitmap will be used + */ + int adjusted_num_chunks = (num_chunks % BITS_PER_BYTE == 0) + ? num_chunks + : (num_chunks + (BITS_PER_BYTE - (num_chunks % BITS_PER_BYTE))); + + /* we need to increase the chunk_size + * so chunks will be CHERI aligned + * (i.e. 16 bytes for RISC-V 64-bit arch) + */ + int adjusted_chunk_size = + (chunk_size % (sizeof(void *)) == 0) + ? chunk_size + : (chunk_size + (sizeof(void *)) - (chunk_size % (sizeof(void *)))); + + /* check this chunk size is small enough so we can represent + * bounds precisely with CHERI compressed representation + */ + adjusted_chunk_size = cheri_representable_length(adjusted_chunk_size); + + /* request memory for our allocation buffer */ + char *res = mmap(NULL, adjusted_num_chunks * adjusted_chunk_size, PROT_READ | PROT_WRITE, + MAP_ANON | MAP_PRIVATE, -1, 0); + /* request memory for our bitmap */ + bitmap = (void *) mmap(NULL, adjusted_num_chunks / BITS_PER_BYTE, + PROT_READ | PROT_WRITE, MAP_ANON | MAP_PRIVATE, -1, 0); + + if (res == MAP_FAILED || bitmap == MAP_FAILED) + { + perror("error in initial mem allocation"); + exit(-1); + } + + /* NB mmap min bounds for capability is 1 page (4K) */ + buffer = res; + /* check buffer is aligned */ + assert((uintptr_t) buffer % sizeof(void *) == 0); + /* check bitmap is aligned */ + assert((uintptr_t) bitmap % sizeof(void *) == 0); + + bytes_per_chunk = adjusted_chunk_size; + buffer_size = adjusted_num_chunks * adjusted_chunk_size; + bitmap_size = adjusted_num_chunks / BITS_PER_BYTE; + + /* zero bitmap, since all chunks are free initially */ + for (i = 0; i < bitmap_size; i++) + { + bitmap[i] = 0; + } + + // set exact bounds for buffer and bitmap? + buffer = cheri_setbounds(buffer, buffer_size); + bitmap = cheri_setbounds(bitmap, bitmap_size); + return; +} + +/* + * allocate fixed size chunk with bitmap allocator + * this is our simplistic `malloc` function + */ +char *alloc_chunk() +{ + unsigned char updated_byte = 0; + int chunk_index = 0; + char *chunk = NULL; + // iterate over all bits in bitmap, looking for a 0 + // when we find a 0, set it to 1 and + // return the corresponding chunk + // (setting its capability bounds) + int i = 0; + while (bitmap[i] == (unsigned char) 0xff) + { + i++; + if (i >= bitmap_size) + break; + } + // do we have a 0? + if (i < bitmap_size && bitmap[i] != (unsigned char) 0xff) + { + // find the lowest 0 ... + int j = 0; + // right shift until bottom bit is 0 + for (j = 0; j < BITS_PER_BYTE; j++) + { + int bit = (bitmap[i] >> j) & 1; + if (bit == 0) + { + break; + } + } + // now i is the word index, j is the bit index + // set this bit to 1 ... + // and work out the chunk to allocate + updated_byte = bitmap[i] + (unsigned char) (1 << j); + bitmap[i] = updated_byte; + + chunk_index = i * BITS_PER_BYTE + j; + chunk = buffer + (chunk_index * bytes_per_chunk); + + /* restrict capability range before returning ptr */ + chunk = cheri_setbounds(chunk, bytes_per_chunk); + } + + return chunk; +} + +void free_chunk(void *chunk) +{ + vaddr_t base = cheri_getbase(chunk); + vaddr_t buff_base = cheri_getbase(buffer); + /* calculate chunk index in buffer */ + int chunk_index = (base - buff_base) / bytes_per_chunk; + assert(chunk_index >= 0); + /* calculate corresponding bitmap index */ + int bitmap_index = chunk_index / BITS_PER_BYTE; + assert(bitmap_index < bitmap_size); + int bitmap_offset = chunk_index % BITS_PER_BYTE; + /* set this bitmap entry to 0 */ + unsigned char updated_byte = bitmap[bitmap_index] & (unsigned char) (~(1 << bitmap_offset)); + bitmap[bitmap_index] = updated_byte; + return; +} + +int num_used_chunks() +{ + int i = 0; + int used_chunks = 0; + + while (i < bitmap_size) + { + unsigned char x = bitmap[i]; + if (x != 0) + { + /* some used chunks here */ + unsigned char j; + for (j = 1; j <= x; j = j << 1) + { + if (x & j) + { + used_chunks++; + } + } + } + i++; + } + return used_chunks; +} + +#endif // STDDEFINES_H_ \ No newline at end of file