added pseudo random number generation (serial, mpi, openmp)

CMakeLists.txt

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

prng/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(random_openmp random_openmp.c)
+    target_link_libraries(random_openmp m)
+else()
+     message("## Skipping 'random_openmp': 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(random_mpi random_mpi.c)
+    target_link_libraries(random_mpi m ${MPI_LIBRARIES})
+else()
+    message("## Skipping 'random_mpi': no MPI support found")
+    dummy_install(${NAME} "MPI")
+endif()
+
+add_executable(random random.c)
+target_link_libraries(random m)
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${C_FLAGS}")
+
+
+# ==================================================================================================

prng/random.c

+// metropolis yöntemi ile random sayı üretir, monte darlo integrasyon
+// compile: gcc random.c -lm
+// usage: ./a.out trial#
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+
+
+// w is gaussian weight function
+double w(double x, double y)
+{
+   double sigma ;
+// sigma is standard deviation
+   sigma = 0.05 ;
+
+   return exp(-(x *x + y*y) / (2 * sigma * sigma) ) ;
+
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.2 ;
+
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update x by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   xt = (*x) + delta * (2 * tmp - 1) ; 
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update y by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   yt = (*y) + delta * (2 * tmp - 1) ; 
+
+
+// compare updated x,y values with old x,y values, accept or reject updated values as new values according to weight function
+   ratio = w(xt, yt ) / w(*x, *y) ;
+   tmp = (double)rand() / (double)RAND_MAX  ;
+
+   if(ratio > tmp )
+   {
+
+      *x = xt ;
+      *y = yt ;
+   }
+
+
+}
+
+
+int main(int argc, char* argv[])
+{
+
+   int i;
+   double x,y ;
+
+// how many random number will generated 
+   int  count = atoi(argv[1]);
+
+//assign initial random numbers between 0-1 to x,y
+
+   x = (double) rand() / RAND_MAX ;
+   y = (double) rand() / RAND_MAX ;
+
+   for(i=0 ; i < count ; i++)
+   {
+
+      printf("%f %f\n", x, y) ;
+      generator(&x,&y);
+   }
+
+
+   return 0;
+
+}

prng/random_mpi.c

+// metropolis yöntemi ile random sayı üretir it generates random numbers by using metropolis method
+// compile: mpicc metrop.c -lm
+// usage: mpirun -n process# ./a.out trial#
+// note: process# must be square of an integer
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <mpi.h>
+
+
+// w is gaussian weight function
+double w(double x, double y)
+{
+   double sigma ;
+// sigma is standard deviation
+   sigma = 0.05 ;
+
+   return exp(-(x *x + y*y) / (2 * sigma * sigma) ) ;
+
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.05 ;
+
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update x by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   xt = (*x) + delta * (2 * tmp - 1) ; 
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update y by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   yt = (*y) + delta * (2 * tmp - 1) ; 
+
+
+// compare updated x,y values with old x,y values, accept or reject updated values as new values according to weight function
+   ratio = w(xt, yt ) / w(*x, *y) ;
+   tmp = (double)rand() / (double)RAND_MAX  ;
+
+   if(ratio > tmp )
+   {
+
+      *x = xt ;
+      *y = yt ;
+   }
+
+
+}
+
+
+int main(int argc, char* argv[])
+{
+   double  x, y ;
+   double * x_array ;
+   double * y_array ;
+
+   double * local_x_array ;
+   double * local_y_array ;
+
+
+   int  count = atoi(argv[1]); // how many random number will generated 
+
+   int i ;
+   int pnumber, myrank; // pnumber is total process number, myrank is rank of each process
+
+   MPI_Status status;
+   MPI_Init(&argc,&argv);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+   MPI_Comm_size(MPI_COMM_WORLD, &pnumber);
+
+
+
+
+
+
+   local_x_array = malloc( count * sizeof(double)) ;
+   local_y_array = malloc( count * sizeof(double)) ;
+
+
+
+// define boundries for each process
+   double x_min, x_max, y_min, y_max ;
+
+   double interval = 2.0 / sqrt(pnumber) ; 
+   x_min = -1.0 + interval * (myrank % ((int) sqrt(pnumber))) ;
+   x_max = x_min + interval ;
+   y_min = -1.0 + interval * (myrank / (int) sqrt(pnumber)) ;
+   y_max = y_min + interval ;
+
+   int counter = 0 ; // counter holds number of random points in local_x_array and local_y_array
+
+   srand(myrank) ; // seed each process a different number to generate different random number
+
+
+
+
+
+   while(counter != 1)
+   {
+
+      x = x_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to x
+      y = y_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to y
+
+
+
+      if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
+      { 
+         local_x_array[0] = x ;
+         local_y_array[0] = y ;
+         counter++ ;
+      }
+
+   }
+
+
+   for(i=1 ; i < count ; i++ )
+   {
+      generator(&x, &y);
+      if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
+      { 
+         local_x_array[counter] = x ;
+         local_y_array[counter] = y ;
+         counter++ ;
+      }
+
+   }
+
+
+
+// compute total generated numbers
+
+   int total_number[pnumber];
+   int master = 0 ;
+
+//all processes will know how many number generated in each process
+   MPI_Allgather(&counter, 1, MPI_INT, total_number, 1, MPI_INT,  MPI_COMM_WORLD);
+
+
+
+//compute displacements from initial element in receive buffer for each send buffer
+   int displs[pnumber] ;
+   displs[0] = 0 ;
+   int sum = 0;
+   for(i=1 ; i < pnumber  ; i++)
+   {
+      sum = sum + total_number[i-1] ;
+      displs[i] = sum ;
+   }
+
+
+   sum = sum + total_number[pnumber-1] ;
+
+   x_array = malloc(sizeof(double) * sum );
+   y_array = malloc(sizeof(double) * sum );
+
+
+
+
+
+   MPI_Gatherv(local_x_array, counter, MPI_DOUBLE, x_array, total_number, displs,  MPI_DOUBLE, master, MPI_COMM_WORLD);
+   MPI_Gatherv(local_y_array, counter, MPI_DOUBLE, y_array, total_number, displs,  MPI_DOUBLE, master, MPI_COMM_WORLD);
+
+
+
+   MPI_Finalize();
+
+   return 0;
+
+}

prng/random_openmp.c

+// metropolis yöntemi ile random sayı üretir it generates random numbers by using metropolis method
+// compile: gcc  -lm -fopenmp random_openmp.c
+// usage: ./a.out  trial# thread#
+// note: process# must be square of an integer
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <omp.h>
+
+
+// w is gaussian weight function
+double w(double x, double y)
+{
+   double sigma ;
+// sigma is standard deviation
+   sigma = 0.05 ;
+
+   return exp(-(x *x + y*y) / (2 * sigma * sigma) ) ;
+
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.05 ;
+
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update x by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   xt = (*x) + delta * (2 * tmp - 1) ; 
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+// update y by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0
+   yt = (*y) + delta * (2 * tmp - 1) ; 
+
+
+// compare updated x,y values with old x,y values, accept or reject updated values as new values according to weight function
+   ratio = w(xt, yt ) / w(*x, *y) ;
+   tmp = (double)rand() / (double)RAND_MAX  ;
+
+   if(ratio > tmp )
+   {
+
+      *x = xt ;
+      *y = yt ;
+   }
+
+
+}
+
+
+int main(int argc, char* argv[])
+{
+   double  x, y ;
+   double * x_array ;
+   double * y_array ;
+
+
+
+   int  count = atoi(argv[1]); // how many random number will generated 
+  
+
+   int tnumber; // tnumber is total thread number, myrank is rank of each thread
+   tnumber = atoi(argv[2]);
+
+   int total_number[tnumber];  // this array holds how many number generated in each thread, and this is shared among all threads
+   
+
+   #pragma omp parallel num_threads(tnumber) 
+   {
+      double * local_x_array ;
+      double * local_y_array ;
+      int myrank = omp_get_thread_num();
+      int i ;
+
+      local_x_array = malloc( count * sizeof(double)) ;
+      local_y_array = malloc( count * sizeof(double)) ;
+
+
+
+   // define boundries for each thread
+      double x_min, x_max, y_min, y_max ;
+
+      double interval = 2.0 / sqrt(tnumber) ; 
+      x_min = -1.0 + interval * (myrank % ((int) sqrt(tnumber))) ;
+      x_max = x_min + interval ;
+      y_min = -1.0 + interval * (myrank / (int) sqrt(tnumber)) ;
+      y_max = y_min + interval ;
+
+      int counter = 0 ; // counter holds number of random points in local_x_array and local_y_array
+
+      srand(myrank) ; // seed each process a different number to generate different random number
+
+
+      while(counter != 1)
+      {
+
+         x = x_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to x
+         y = y_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to y
+
+
+
+         if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
+         { 
+            local_x_array[0] = x ;
+            local_y_array[0] = y ;
+            counter++ ;
+         }
+
+      }
+
+
+      for(i=1 ; i < count ; i++ )
+      {
+         generator(&x, &y);
+         if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
+         { 
+            local_x_array[counter] = x ;
+            local_y_array[counter] = y ;
+            counter++ ;
+         }
+
+      }
+
+
+
+      total_number[myrank] = counter ;
+
+// to wait until total_number array is populated place a barrier 
+
+#pragma omp barrier
+
+
+      if (myrank == 0)
+      {
+         int sum = 0 ; // sum holds sum of generated numbers by each thead 
+         for(i=0 ; i < tnumber ; i++)
+         {
+            sum = sum + total_number[i] ;
+         }
+      printf("sum is %d\n", sum);
+         x_array = malloc(sizeof(double) * sum);
+         y_array = malloc(sizeof(double) * sum);
+
+      }
+
+// to prevent threads to write before thread 0 create x_array, y_array place another barrier
+#pragma omp barrier
+
+   // now each array populates global x_array, y_array by using private array local_x_array, local_y_array   
+   // firstly, each thread computes displacement from initial element of global array 
+      int displs = 0 ;
+      
+      for(i=0 ; i < myrank ; i++)
+      {
+         displs = displs + total_number[i];
+      }
+
+
+      for(i=0 ; i < counter ; i++ )
+      {
+
+         x_array[displs + i] = local_x_array[i] ;
+         y_array[displs + i] = local_y_array[i] ;
+
+      }
+
+   } //end pragma
+
+
+   return 0;
+
+}