Cheri-research/FAT-Allocator
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Ported XSBench

Browse Source
This commit is contained in:
Akilan Selvacoumar
2024-12-02 14:32:04 +00:00
parent 1ea61d9c54
commit 81397da55e
14 changed files with 1911 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

230
benchmarks/benchmarks/XSbench/GridInit.c Normal file
View File

@@ -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;
}

2
benchmarks/benchmarks/XSbench/Main.c
View File

@@ -120,4 +120,4 @@ int main( int argc, char* argv[] )
#endif
return is_invalid_result;
}
}

2
benchmarks/benchmarks/XSbench/Makefile
View File

@@ -20,7 +20,7 @@ program = XSBench
source = \
Main.c \
io.c \
Simulation.c \
Simulations.c \
GridInit.c \
XSutils.c \
Materials.c

2
benchmarks/benchmarks/XSbench/Materials.c
View File

@@ -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;
}

871
benchmarks/benchmarks/XSbench/Simulations.c
View File

@@ -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 <optimized variant ID>" 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 <optimized variant ID>" 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)<cutoff) ){
if (left < j){ quickSort_parallel_internal_i_d(key, value, left, j, cutoff); }
if (i < right){ quickSort_parallel_internal_i_d(key, value, i, right, cutoff); }
}else{
#pragma omp task
{ quickSort_parallel_internal_i_d(key, value, left, j, cutoff); }
#pragma omp task
{ quickSort_parallel_internal_i_d(key, value, i, right, cutoff); }
}
}
void quickSort_parallel_i_d(int* key,double * value, int lenArray, int numThreads){
// Set minumum problem size to still spawn threads for
int cutoff = 10000;
// For this problem size, more than 16 threads on CPU is not helpful
if( numThreads > 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)<cutoff) ){
if (left < j){ quickSort_parallel_internal_d_i(key, value, left, j, cutoff); }
if (i < right){ quickSort_parallel_internal_d_i(key, value, i, right, cutoff); }
}else{
#pragma omp task
{ quickSort_parallel_internal_d_i(key, value, left, j, cutoff); }
#pragma omp task
{ quickSort_parallel_internal_d_i(key, value, i, right, cutoff); }
}
}
void quickSort_parallel_d_i(double* key,int * value, int lenArray, int numThreads){
// Set minumum problem size to still spawn threads for
int cutoff = 10000;
// For this problem size, more than 16 threads on CPU is not helpful
if( numThreads > 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;
}

2
benchmarks/benchmarks/XSbench/XSbench_header.h
View File

@@ -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

2
benchmarks/benchmarks/XSbench/XSutils.c
View File

@@ -60,4 +60,4 @@ double get_time(void)
double time = (double) ms / 1000.0;
return time;
}
}

2
benchmarks/benchmarks/XSbench/io.c
View File

@@ -505,4 +505,4 @@ SimulationData binary_read( Inputs in )
fclose(fp);
return SD;
}
}

395
benchmarks/benchmarks/kmeans/bitmap.h Normal file
View File

@@ -0,0 +1,395 @@
/*
* Copyright (c) 2011 Bharath Ramesh <bramesh.dev@gmail.com>
*
* Distributed under the terms of GNU LGPL, version 2.1
*/
#include <stdint.h>
#include <string.h>
// #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 <stdio.h>
#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 <sys/mman.h>
#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));
}

5
benchmarks/benchmarks/kmeans/kmeans-pthread.c
View File

@@ -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;

177
benchmarks/benchmarks/kmeans/simple.h Normal file
View File

@@ -0,0 +1,177 @@
/* Copyright (C) 2023. Shivashish Das. Licensed under the MIT License.*/
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
// source: https://www.reddit.com/r/C_Programming/comments/1bt8dyz/github_dasshivamalloc_a_simple_memory_allocator/
// #include "alloc.h"
#ifndef _MSC_VER
#include <sys/mman.h>
#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;
}

266
benchmarks/benchmarks/kmeans/stddefines.h
View File

@@ -48,9 +48,37 @@
#include <sys/stat.h>
#include <fcntl.h>
// #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_