-
Notifications
You must be signed in to change notification settings - Fork 58
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
sw: Add sample Monte Carlo application
- Loading branch information
Showing
6 changed files
with
172 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,59 @@ | ||
// Copyright 2023 ETH Zurich and University of Bologna. | ||
// Licensed under the Apache License, Version 2.0, see LICENSE for details. | ||
// SPDX-License-Identifier: Apache-2.0 | ||
// | ||
// Fabio Cappellini <[email protected]> | ||
// Hakim Filali <[email protected]> | ||
// Luca Colagrande <[email protected]> | ||
|
||
#include "lcg.h" | ||
#include "math.h" | ||
#include "pi_estimation.h" | ||
#include "snrt.h" | ||
|
||
#define N_SAMPLES 1024 | ||
|
||
__thread uint32_t seed0, seed1, Ap, Cp; | ||
|
||
double pi_estimate; | ||
|
||
inline void mc_init() { | ||
// Double the sequences as each core produces two random numbers | ||
unsigned int num_sequences = 2 * snrt_cluster_compute_core_num(); | ||
init_2d_lcg_params(num_sequences, 0, LCG_A, LCG_C, &seed0, &seed1, &Ap, | ||
&Cp); | ||
} | ||
|
||
int main() { | ||
// Initialize the PRNGs for parallel Monte Carlo | ||
if (snrt_is_compute_core()) mc_init(); | ||
|
||
// Store partial sum array at first free address in TCDM | ||
uint32_t* reduction_array = (uint32_t*)snrt_l1_next(); | ||
|
||
// Calculate partial sums | ||
uint32_t n_samples_per_core = N_SAMPLES / snrt_cluster_compute_core_num(); | ||
reduction_array[snrt_cluster_core_idx()] = | ||
calculate_partial_sum(seed0, seed1, Ap, Cp, n_samples_per_core); | ||
|
||
// Synchronize cores | ||
snrt_cluster_hw_barrier(); | ||
|
||
// First core in cluster performs the final calculation | ||
if (snrt_cluster_core_idx() == 0) { | ||
// Reduce partial sums | ||
uint32_t sum = 0; | ||
for (int i = 0; i < snrt_cluster_compute_core_num(); i++) { | ||
sum += reduction_array[i]; | ||
} | ||
|
||
// Estimate pi | ||
pi_estimate = estimate_pi(sum, N_SAMPLES); | ||
|
||
// Check result | ||
double err = fabs(pi_estimate - M_PI); | ||
if (err > 0.05) return 1; | ||
} | ||
|
||
return 0; | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,34 @@ | ||
// Copyright 2023 ETH Zurich and University of Bologna. | ||
// Licensed under the Apache License, Version 2.0, see LICENSE for details. | ||
// SPDX-License-Identifier: Apache-2.0 | ||
// | ||
// Fabio Cappellini <[email protected]> | ||
// Hakim Filali <[email protected]> | ||
// Luca Colagrande <[email protected]> | ||
|
||
__thread double one = 1.0; | ||
|
||
inline uint32_t calculate_partial_sum(uint32_t seed0, uint32_t seed1, | ||
uint32_t Ap, uint32_t Cp, | ||
unsigned int n_samples) { | ||
uint32_t x1 = seed0; | ||
uint32_t x2 = seed1; | ||
double u1, u2; | ||
unsigned int result = 0; | ||
|
||
for (unsigned int i = 0; i < n_samples; i++) { | ||
u1 = normalize(x1); | ||
u2 = normalize(x2); | ||
x1 = lcg(Ap, Cp, x1); | ||
x2 = lcg(Ap, Cp, x2); | ||
|
||
if ((u1 * u1 + u2 * u2) < one) { | ||
result++; | ||
} | ||
} | ||
return result; | ||
} | ||
|
||
inline double estimate_pi(uint32_t sum, uint32_t n_samples) { | ||
return (double)(4 * sum) / (double)n_samples; | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,63 @@ | ||
// Copyright 2023 ETH Zurich and University of Bologna. | ||
// Licensed under the Apache License, Version 2.0, see LICENSE for details. | ||
// SPDX-License-Identifier: Apache-2.0 | ||
// | ||
// Fabio Cappellini <[email protected]> | ||
// Hakim Filali <[email protected]> | ||
// Luca Colagrande <[email protected]> | ||
|
||
#include "snrt.h" | ||
|
||
// Numerical Recipes from the "quick and dirty generators" list, Chapter 7.1, | ||
// Eq. 7.1.6 parameters from Knuth and H. W. Lewis | ||
#define MAX_UINT_PLUS1 4294967296.0 | ||
#define LCG_A 1664525 | ||
#define LCG_C 1013904223 | ||
|
||
__thread double max_uint_plus_1_inverse = (double)1.0 / (double)MAX_UINT_PLUS1; | ||
|
||
// Calculate A' and C' constants | ||
inline void leapfrog_constants(unsigned int num_sequences, uint32_t a, | ||
uint32_t c, uint32_t* Ap, uint32_t* Cp) { | ||
// Ap = a^k | ||
// Cp = (a^(k-1) + a^(k-2) + ... + a^1 + a^0)*c | ||
uint32_t Ap_tmp = 1; | ||
uint32_t Cp_tmp = 0; | ||
for (unsigned int p = 0; p < num_sequences; p++) { | ||
Cp_tmp += Ap_tmp; | ||
Ap_tmp *= a; | ||
} | ||
Cp_tmp *= c; | ||
|
||
// Store temporary variables to outputs | ||
*Ap = Ap_tmp; | ||
*Cp = Cp_tmp; | ||
} | ||
|
||
// Calculate seed for leapfrog method's right sequence of index `sequence_idx` | ||
inline uint32_t right_seed(uint32_t left_seed, uint32_t a, uint32_t c, | ||
unsigned int sequence_idx) { | ||
uint32_t seed = left_seed; | ||
for (unsigned int p = 1; p <= sequence_idx; p++) { | ||
seed = seed * a + c; | ||
} | ||
return seed; | ||
} | ||
|
||
// Generate next PRN from LCG recurrence equation | ||
inline uint32_t lcg(uint32_t a, uint32_t c, uint32_t previous) { | ||
return previous * a + c; | ||
} | ||
|
||
// Normalize LCG PRN to [0, 1) range | ||
inline double normalize(uint32_t x) { | ||
return (double)x * max_uint_plus_1_inverse; | ||
} | ||
|
||
inline void init_2d_lcg_params(unsigned int num_sequences, uint32_t seed, | ||
uint32_t a, uint32_t c, uint32_t* seed0, | ||
uint32_t* seed1, uint32_t* Ap, uint32_t* Cp) { | ||
*seed0 = right_seed(seed, a, c, 2 * snrt_global_compute_core_idx()); | ||
*seed1 = right_seed(seed, a, c, 2 * snrt_global_compute_core_idx() + 1); | ||
leapfrog_constants(num_sequences, a, c, Ap, Cp); | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
14 changes: 14 additions & 0 deletions
14
target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,14 @@ | ||
# Copyright 2023 ETH Zurich and University of Bologna. | ||
# Licensed under the Apache License, Version 2.0, see LICENSE for details. | ||
# SPDX-License-Identifier: Apache-2.0 | ||
# | ||
# Luca Colagrande <[email protected]> | ||
|
||
PI_ESTIMATION_DIR = ../../../../../../sw/apps/montecarlo/pi_estimation | ||
PRNG_DIR = ../../../../../../sw/apps/prng | ||
|
||
APP ?= pi_estimation | ||
SRCS ?= $(PI_ESTIMATION_DIR)/main.c | ||
INCDIRS += $(PI_ESTIMATION_DIR) $(PRNG_DIR) | ||
|
||
include ../../common.mk |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters