From 60b8f6d237bf99c9ec0504f32d7fe0d14332a5b6 Mon Sep 17 00:00:00 2001 From: Luca Colagrande Date: Tue, 26 Sep 2023 18:19:13 +0200 Subject: [PATCH] sw: Add sample Monte Carlo application --- sw/apps/montecarlo/pi_estimation/main.c | 59 +++++++++++++++++ .../montecarlo/pi_estimation/pi_estimation.h | 34 ++++++++++ sw/apps/prng/lcg.h | 63 +++++++++++++++++++ target/snitch_cluster/sw/apps/Makefile | 1 + .../sw/apps/montecarlo/pi_estimation/Makefile | 14 +++++ target/snitch_cluster/sw/run.yaml | 1 + 6 files changed, 172 insertions(+) create mode 100644 sw/apps/montecarlo/pi_estimation/main.c create mode 100644 sw/apps/montecarlo/pi_estimation/pi_estimation.h create mode 100644 sw/apps/prng/lcg.h create mode 100644 target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile diff --git a/sw/apps/montecarlo/pi_estimation/main.c b/sw/apps/montecarlo/pi_estimation/main.c new file mode 100644 index 000000000..f5a279bec --- /dev/null +++ b/sw/apps/montecarlo/pi_estimation/main.c @@ -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 +// Hakim Filali +// Luca Colagrande + +#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; +} diff --git a/sw/apps/montecarlo/pi_estimation/pi_estimation.h b/sw/apps/montecarlo/pi_estimation/pi_estimation.h new file mode 100644 index 000000000..99da3b532 --- /dev/null +++ b/sw/apps/montecarlo/pi_estimation/pi_estimation.h @@ -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 +// Hakim Filali +// Luca Colagrande + +__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; +} diff --git a/sw/apps/prng/lcg.h b/sw/apps/prng/lcg.h new file mode 100644 index 000000000..19ffb2cbe --- /dev/null +++ b/sw/apps/prng/lcg.h @@ -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 +// Hakim Filali +// Luca Colagrande + +#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); +} diff --git a/target/snitch_cluster/sw/apps/Makefile b/target/snitch_cluster/sw/apps/Makefile index d0c410aa5..d86145c85 100644 --- a/target/snitch_cluster/sw/apps/Makefile +++ b/target/snitch_cluster/sw/apps/Makefile @@ -17,6 +17,7 @@ SUBDIRS += dnn/layernorm SUBDIRS += dnn/linear SUBDIRS += dnn/maxpool SUBDIRS += dnn/softmax +SUBDIRS += montecarlo/pi_estimation .PHONY: all clean $(SUBDIRS) diff --git a/target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile b/target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile new file mode 100644 index 000000000..ed455b3b5 --- /dev/null +++ b/target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile @@ -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 + +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 diff --git a/target/snitch_cluster/sw/run.yaml b/target/snitch_cluster/sw/run.yaml index 8e80cf35c..f25ea7641 100644 --- a/target/snitch_cluster/sw/run.yaml +++ b/target/snitch_cluster/sw/run.yaml @@ -84,3 +84,4 @@ runs: # throws illegal instruction on FDIV in simulation # - elf: apps/dnn/softmax/build/softmax.elf # throws illegal instruction on FDIV in simulation + - elf: apps/montecarlo/pi_estimation/build/pi_estimation.elf