1d integration using monte carlo (serial, mpi, openmp)

CMakeLists.txt

@@ -24,6 +24,7 @@ set(DWARF_PREFIX 7_montecarlo) # The prefix of the name of the binaries produced
 # Add the examples
 add_subdirectory(pi)
 add_subdirectory(prng)
+add_subdirectory(integral1d)
 
 # ==================================================================================================
 

integral1d/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(integral1d_OMP integral1d_OMP.c)
+    target_link_libraries(integral1d_OMP m)
+else()
+     message("## Skipping 'integral1d_OMP': no OpenMP support found")
+     dummy_install(${NAME} "OpenMP")
+endif()
+
+if (MPI_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS}")
+    include_directories(${MPI_INCLUDE_PATH})
+    add_executable(integral1d_mpi integral1d_mpi.c)
+    target_link_libraries(integral1d_mpi m ${MPI_LIBRARIES})
+else()
+    message("## Skipping 'integral1d_mpi': no MPI support found")
+    dummy_install(${NAME} "MPI")
+endif()
+
+add_executable(integral1d integral1d.c)
+target_link_libraries(integral1d m)
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${C_FLAGS}")
+
+
+# ==================================================================================================

integral1d/integral1d.c

+// weight function kullanarak monte carlo integrasyon
+
+/*
+
+COMPILATION: gcc integral1d.c -lm
+USAGE: ./a.out trial#
+
+This program computes integral of ( 1 / (1 + x*x) ) in x=[0,1] by using importance sampling Monte Carlo 
+
+
+*/
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+
+
+// f() is the function to be integrated
+double f(double x)
+{
+
+   return (1.0 / (1.0 + x*x) ) ;
+
+}
+
+
+// random x values are generated according to probability distribution function w() 
+double w(double x)
+{
+
+
+   return ((1.0/3.0) * (4 - 2 * x)) ;
+}
+
+
+// my_rand() function generates random numbers according to p.d.f w()
+// this function uses inverse of the p.d.f w() 
+double my_rand()
+{
+
+   double y = (double)rand() / (double)RAND_MAX ;
+
+   return (2.0 - sqrt(4.0 - 3.0*y) ) ;
+}
+
+
+int main(int argc, char* argv[])
+{
+   int trials;  // number of random points
+   int i;
+   double sum = 0.0;
+   double x;
+
+
+   trials = atoi(argv[1]); // read number of random point from command line argument
+
+
+   for(i=0 ; i < trials ; i++ )
+   {
+      x = my_rand()  ;
+    
+      sum += ( f(x) / w(x)) ;
+
+   }
+
+   double integral = sum / trials ;
+
+   printf("computed value of the integral= %f \n", integral);
+   printf("true value of the integral=     %f \n", atan(1.0));
+   printf("error=                          %f  \n", fabs( integral - atan(1.0)));
+
+
+return 0;
+
+}
+

integral1d/integral1d_OMP.c

+/*
+
+COMPILATION: gcc  -fopenmp integral1d_OMP.c -lm
+USAGE: ./a.out thread#  trial#
+
+This program computes integral of ( 1 / (1 + x*x) ) in x=[0,1] by using importance sampling Monte Carlo 
+
+*/
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <omp.h>
+
+
+// f() is the function to be integrated
+double f(double x)
+{
+
+   return (1.0 / (1.0 + x*x) ) ;
+
+}
+
+// random x values are generated according to probability distribution function w() 
+double w(double x)
+{
+
+
+   return ((1.0/3.0) * (4 - 2 * x)) ;
+}
+
+// my_rand() function generates random numbers according to p.d.f w()
+// this function uses inverse of the p.d.f w() 
+double my_rand()
+{
+
+   double y = (double)rand() / (double)RAND_MAX ;
+
+   return (2.0 - sqrt(4.0 - 3.0*y) ) ;
+}
+
+
+int main(int argc, char* argv[])
+{
+
+
+   int i;
+   double sum = 0;
+   double x;
+
+
+
+   int thread_count = atoi(argv[1]); // number of threads participating computation is read from the command lile arguments
+   int trials = atoi(argv[2]); // trials is total number of random points, it is read from the command line arguments
+
+
+// multiple threads are started after this compiler directive, num_threads() function indicates number of threads to be started
+// reduction(+:sum) function indicates end of the parallel region total_hits variables of each thread are gathered in master thread
+// private(x) indicates variable x is private to each thread instead of being shared by threads
+   #pragma omp parallel num_threads(thread_count) reduction(+:sum) private(x)
+   {
+      // It is desired every thread to generate different random numbers so
+      // each thread are seeded with different number (thread id)
+      srand(omp_get_thread_num());
+
+      //this compiler directive divides the loop iterations between the spawned threads
+      #pragma omp  for       
+      for(i=0 ; i < trials ; i++ )
+      {
+         x = my_rand()  ;
+         sum += ( f(x) / w(x) ) ;
+
+      }
+      
+   }// all threads are joined after this block
+
+
+   double integral = sum / trials ;
+
+   printf("computed value of the integral= %f \n", integral);
+   printf("true value of the integral=     %f \n", atan(1.0));
+   printf("error=                          %f  \n", fabs( integral - atan(1.0)));
+
+
+   return 0;
+
+}
+

integral1d/integral1d_mpi.c

+
+/*
+
+COMPILATION: mpicc integral1d_mpi.c -lm
+USAGE: mpirun -n processor# ./a.out trial#
+
+This program computes integral of ( 1 / (1 + x*x) ) in x=[0,1] by using importance sampling Monte Carlo 
+
+
+*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <mpi.h>
+
+// f() is the function to be integrated
+double f(double x)
+{
+
+   return (1.0 / (1.0 + x*x) ) ;
+
+}
+
+// random x values are generated according to probability distribution function w() 
+double w(double x)
+{
+
+
+   return ((1.0/3.0) * (4 - 2 * x)) ;
+}
+
+
+// my_rand() function generates random numbers according to p.d.f w()
+// this function uses inverse of the p.d.f w() 
+double my_rand()
+{
+
+
+   double y = (double)rand() / (double)RAND_MAX ;
+
+   return (2.0 - sqrt(4.0 - 3.0*y) ) ;
+}
+
+
+int main(int argc, char* argv[])
+{
+   int trials, trials_per_p; // trials is total number of random points and trials_per_p number of random points for each processor 
+   int i;
+   double partial_sum = 0; // each processor uses its variable partial_sum to sum function value
+   double total_sum; // at the end of the computation all processors' partial_sum gathered in total_sum
+   double x;
+   int pnumber, myrank; // pnumber is the number of processor in the computation and myrank is the rank of the processor
+
+   MPI_Init(&argc,&argv);  // Initialize the MPI execution environment
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank); // Each processor learns its rank number
+   MPI_Comm_size(MPI_COMM_WORLD, &pnumber); // Each processor learns how many processor participate to the calculation
+
+   trials = atoi(argv[1]); // command line argument is readed for total random number
+
+   trials_per_p = trials / pnumber ; // number of points for each processor is calculated
+
+
+   srand(myrank); // It is desired every processor to generate different random numbers so
+                  //      each processor seeded with different number (rank   oprocessor)
+
+
+   for(i=0 ; i < trials_per_p ; i++ )
+   {
+      x = my_rand()  ;
+    
+      partial_sum += ( f(x) / w(x) ) ;
+
+   }
+
+   // Assume the processor which has rank=0 is master
+   int master_p = 0;
+   
+   //total sum and total random number are gathered to the master processor   
+   MPI_Reduce(&partial_sum,  &total_sum, 1,  MPI_DOUBLE,  MPI_SUM, master_p, MPI_COMM_WORLD);
+   MPI_Reduce(&trials,  &trials_per_p, 1,  MPI_INT,  MPI_SUM, master_p, MPI_COMM_WORLD);
+
+   //Master processor prints the result
+   if(myrank == master_p)
+   {
+      double integral = total_sum / trials ;
+
+      printf("computed value of the integral= %f \n", integral);
+      printf("true value of the integral=     %f \n", atan(1.0));
+      printf("error=                          %f  \n", fabs( integral - atan(1.0)));
+
+   }
+
+
+   MPI_Finalize(); //Terminates MPI execution environment
+
+
+return 0;
+
+}
+