diff --git a/answers/01/pi_monte_carlo_parallel.cpp b/answers/01/pi_monte_carlo_parallel.cpp new file mode 100644 index 0000000..447ed24 --- /dev/null +++ b/answers/01/pi_monte_carlo_parallel.cpp @@ -0,0 +1,52 @@ +#include +#include +#include + +using namespace std; + +int main() { + unsigned int seed = 0; + default_random_engine re{seed}; + uniform_real_distribution zero_to_one{0.0, 1.0}; + + int n = 100000000; // number of points to generate + int counter = 0; // counter for points lying in the first quadrant of a unit circle + auto start_time = omp_get_wtime(); // omp_get_wtime() is an OpenMP library routine + +#pragma omp parallel + { + int num_threads = omp_get_num_threads(); + int thread_id = omp_get_thread_num(); + + int thread_shift = n / num_threads; // the remainder will just get handled by the last thread + int thread_start = thread_shift*thread_id; + int thread_end = thread_id == num_threads-1 ? n : thread_shift*(thread_id+1); +#pragma omp critical + { // Just a pretty lil debug print out to make sure my division logic is logicing + std::cout << "Thread " << thread_id << " of " << num_threads << ": Computing " << thread_start << " to " << thread_end << std::endl; + } + + int thread_counter = 0; + + // compute n points and test if they lie within the first quadrant of a unit circle + for (int i = thread_start; i < thread_end; i++) { + auto x = zero_to_one(re); // generate random number between 0.0 and 1.0 + auto y = zero_to_one(re); // generate random number between 0.0 and 1.0 + if (x * x + y * y <= 1.0) { // if the point lies in the first quadrant of a unit circle + thread_counter++; + } + } + +#pragma omp critical + { + counter = counter+thread_counter; + } + } + + auto run_time = omp_get_wtime() - start_time; + auto pi = 4 * (double(counter) / n); + + cout << "pi: " << pi << endl; + cout << "run_time: " << run_time << " s" << endl; + cout << "n: " << n << endl; +}