added mc approximation for pi (serial, mpi, openmp)

CMakeLists.txt

+
+# ==================================================================================================
+# This file is part of the CodeVault project. The project is licensed under Apache Version 2.0.
+# CodeVault is part of the EU-project PRACE-4IP (WP7.3.C).
+#
+# Author(s):
+#   Samet Demir (@itu)
+#
+# ==================================================================================================
+
+# CMake project
+cmake_minimum_required(VERSION 2.8.10 FATAL_ERROR)
+project("7_montecarlo")
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../cmake/common.cmake)
+
+# ==================================================================================================
+
+# Dwarf 7: Monte Carlo
+message("--------------------")
+message("Dwarf 7: Monte Carlo:")
+message("--------------------")
+set(DWARF_PREFIX 7_montecarlo) # The prefix of the name of the binaries produced
+
+# Add the examples
+add_subdirectory(pi)
+
+# ==================================================================================================
+
+
+
+

pi/CMakeLists.txt

+# Packages are optional: if they are not present, certain code samples are not compiled
+find_package(OpenMP)   # Built-in in CMake
+find_package(MPI)      # Built-in in CMake
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../cmake/common.cmake)
+
+# ==================================================================================================
+
+if ("${DWARF_PREFIX}" STREQUAL "")
+  set(DWARF_PREFIX 7_montecarlo)
+endif()
+
+# C compiler settings
+
+find_package(Common)
+
+if (OPENMP_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    add_executable(pi_omp pi_omp.c pi_shared.c)
+else()
+     message("## Skipping 'pi_omp': no OpenMP support found")
+     dummy_install(${NAME} "OpenMP")
+endif()
+
+if (MPI_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    include_directories(${MPI_INCLUDE_PATH})
+    add_executable(pi_mpi pi_mpi.c pi_shared.c)
+    target_link_libraries(pi_mpi ${MPI_LIBRARIES})
+endif()
+
+add_executable(pi_serial pi_serial.c pi_shared.c)
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${C_FLAGS}")
+
+
+# ==================================================================================================

pi/pi_mpi.c

+#include <stdio.h> // printf
+#include <stdlib.h> // EXIT_SUCCESS, srand
+#include <math.h> // M_PI, fabs
+#include <mpi.h>
+
+#include "pi_shared.h"
+
+int main(int argc, char* argv[]) {
+    int seed;
+    int N;
+    int hits;
+    int hits_sum;
+    double my_pi;
+    int number_of_workers;
+    int this_worker;
+    int N_per_worker;
+
+    // initialize MPI sybsystem
+    MPI_Init(&argc,&argv);
+
+    // get number of available processes, and a unique id for this process
+    MPI_Comm_rank(MPI_COMM_WORLD, &this_worker);
+    MPI_Comm_size(MPI_COMM_WORLD, &number_of_workers);
+
+    // get parameters from command line arguments
+    read_arguments(&seed, &N, argc, argv);
+
+    // divide the workload equally among processes
+    N_per_worker = N / number_of_workers;
+    N = N_per_worker * number_of_workers;
+
+    // generate N_per_worker random coordinates, get number of hits
+    // (see count_hits() function for details)
+    // please note that every process will start with a different seed value
+    hits = count_hits(seed + this_worker, N_per_worker);
+
+    // sum hits values from different processes
+    MPI_Reduce(&hits, &hits_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
+
+    // process the combined result (on root process only)
+    if(this_worker == 0) {
+        // estimate PI based on following assumptions
+        //   area of circle                = PI * r^2   ~= hits
+        //   area of square around circle  = (2 * r)^2  ~= N
+        //   r = 1
+        my_pi = 4. * ((double) hits_sum  / (double) N);
+
+        // print seed, N, number_of_workers, result, error
+        printf("%d\t%d\t%d\t%f\t%f\n", seed, N, number_of_workers, my_pi, fabs(my_pi-M_PI));
+
+    }
+
+    // cleanly terminate MPI subsystem
+    MPI_Finalize();
+
+    return EXIT_SUCCESS;
+
+}

pi/pi_omp.c

+#include <stdio.h> // printf
+#include <stdlib.h> // EXIT_SUCCESS
+#include <math.h> // M_PI, fabs
+#include <omp.h>
+
+#include "pi_shared.h"
+
+int main(int argc, char* argv[]) {
+    int seed;
+    int N;
+    int hits;
+    double my_pi;
+    int number_of_workers;
+    int N_per_worker;
+
+    // get parameters from command line arguments
+    read_arguments(&seed, &N, argc, argv);
+
+    // create OMP_NUM_THREADS threads
+    // every thread will have its own hits variable (reduction implies private keyword)
+    // after the parallel section, a single hits variable will remain,
+    // holding the sum of all thread private hits variables
+#pragma omp parallel reduction(+:hits)
+    {
+        // divide the workload equally among threads
+        number_of_workers = omp_get_num_threads();
+        N_per_worker = N / number_of_workers;
+        N = N_per_worker * number_of_workers;
+
+        // generate N_per_worker random coordinates, get number of hits
+        // (see count_hits() function for details)
+        // please note that every thread will start with a different seed value
+        hits = count_hits(seed + omp_get_thread_num(), N_per_worker);
+    } 
+
+    // estimate PI based on following assumptions
+    //   area of circle                = PI * r^2   ~= hits
+    //   area of square around circle  = (2 * r)^2  ~= N
+    //   r = 1
+    my_pi = 4. * ((double) hits  / (double) N);
+
+    // print seed, N, number_of_workers, result, error
+    printf("%d\t%d\t%d\t%f\t%f\n", seed, N, number_of_workers, my_pi, fabs(my_pi-M_PI));
+
+    return EXIT_SUCCESS;
+
+}
+

pi/pi_serial.c

+#include <stdio.h> // printf
+#include <stdlib.h> // EXIT_SUCCESS, srand
+#include <math.h> // M_PI, fabs
+
+#include "pi_shared.h"
+
+int main(int argc, char* argv[]) {
+    int seed;
+    int N;
+    int hits;
+    double my_pi;
+
+    // get parameters from command line arguments
+    read_arguments(&seed, &N, argc, argv);
+
+    // generate N random coordinates, get number of hits
+    // (see count_hits() function for details)
+    hits = count_hits(seed, N);
+
+    // estimate PI based on following assumptions
+    //   area of circle                = PI * r^2   ~= hits
+    //   area of square around circle  = (2 * r)^2  ~= N
+    //   r = 1
+    my_pi = 4. * ((double) hits  / (double) N);
+
+    // print seed, N, number_of_workers, result, error
+    printf("%d\t%d\t%d\t%f\t%f\n", seed, N, 1, my_pi, fabs(my_pi-M_PI));
+
+    return EXIT_SUCCESS;
+
+}
+

pi/pi_shared.c

+#include <stdlib.h> // rand_r, RAND_MAX, EXIT_FAILURE, atoi
+#include <stdio.h> // printf
+#include <time.h> // time
+
+void read_arguments(int *seed, int *N, int argc, char *argv[]) {
+    /*
+     * parse command line arguments to obtain
+     * - seed: seed for the pseudo random number generator
+     * - N: number of random samples
+     */
+
+    // check if enough arguments are given
+    if (argc != 3) { 
+        // print usage information
+        printf("USAGE %s <seed> <N>\n", argv[0]);
+        printf("  seed (integer): seed for the pseudo random number generator\n");
+        printf("                  a negative value indicates current timestamp\n");
+        printf("  N (integer)   : number of random samples to use\n");
+        // do not continue if number of arguments is invalid
+        exit(EXIT_FAILURE);
+    }
+
+    // read seed for the pseudo random numnber generator
+    // WARNING: according to the man page, atoi does not detect errors
+    *seed = atoi(argv[1]);
+    if (*seed < 0) {
+        // use current timestamp as seed
+        *seed = time(NULL);
+    }
+
+    // read number of random samples
+    // WARNING: according to the man page, atoi does not detect errors
+    *N = atoi(argv[2]);
+
+}
+
+
+int count_hits(unsigned int seed, int N) {
+    // create N random ((double) x, (double) y) coordinates, where x:[-1, 1], y:[-1, 1]
+    // return the number of coordinates falling inside a unit circle
+
+    int i;
+    double x;
+    double y;
+    int hits;
+
+    // start with 0 hits
+    hits = 0;
+
+    // for N random samples
+    for (i=0; i<N; i++) {
+
+        // create a random 2D coordinate
+        x = -1. + 2. * (double)rand_r(&seed) / (double)RAND_MAX ;
+        y = -1. + 2. * (double)rand_r(&seed) / (double)RAND_MAX ;
+
+        // if coordinate falls inside unit circle, increase hit counter
+        // ISO C standard guarantees the return value of comparison operators
+        // to be (int) 0 for false and (int) 1 for true results
+        // therefore, we can avoid branching (i.e. if (...) hits++;)
+        hits += ( (x*x + y*y) < 1. );
+
+    }
+
+    // return number of hits
+    return hits;
+
+}

pi/pi_shared.h

+/* pi_shared.c */
+void read_arguments(int *seed, int *N, int argc, char *argv[]);
+int count_hits(unsigned int seed, int N);