added 2d integration using monte carlo, various communication models (mpi, openmp)

CMakeLists.txt

@@ -25,6 +25,7 @@ set(DWARF_PREFIX 7_montecarlo) # The prefix of the name of the binaries produced
 add_subdirectory(pi)
 add_subdirectory(prng)
 add_subdirectory(integral_basic)
+add_subdirectory(integral_advanced)
 
 # ==================================================================================================
 

integral_advanced/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(integral1_openmp integral1_openmp.c)
+    target_link_libraries(integral1_openmp m)
+    add_executable(integral2_openmp integral2_openmp.c)
+    target_link_libraries(integral2_openmp m)
+    add_executable(integral3_openmp integral3_openmp.c)
+    target_link_libraries(integral3_openmp m)
+else()
+     message("## Skipping 'integral 1/2/3 advanced': 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(integral1_mpi integral1_mpi.c)
+    target_link_libraries(integral1_mpi m ${MPI_LIBRARIES})
+    add_executable(integral2_mpi integral2_mpi.c)
+    target_link_libraries(integral2_mpi m ${MPI_LIBRARIES})
+    add_executable(integral3_mpi integral3_mpi.c)
+    target_link_libraries(integral3_mpi m ${MPI_LIBRARIES})
+else()
+    message("## Skipping 'integral 1/2/3 advanced': no MPI support found")
+    dummy_install(${NAME} "MPI")
+endif()
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${C_FLAGS}")
+
+
+# ==================================================================================================

integral_advanced/integral1_mpi.c

+// calculate integral using master-workers model
+// compile:  mpicc integral1_mpi.c -lm
+// usage: mpirun -n process# ./a.out trial#
+
+#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)
+{
+  return exp( - ( x*x + y*y ) );
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.5 ;
+
+
+   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 ;
+   }
+
+
+}
+
+
+//f is the function to be integrated
+double f(double x, double y)
+{
+
+
+return M_PI *  ( x*x + y*y ) ;
+
+}
+
+
+int main(int argc, char* argv[])
+{
+   double  x, y ;
+   double * x_array ;
+   double * y_array ;
+
+   double sum ;
+   double total_sum ;
+   double integral ;
+   int trials, trials_per_p;
+   int i, j ;
+   int pnumber, myrank;
+
+   MPI_Status status;
+   MPI_Init(&argc,&argv);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+   MPI_Comm_size(MPI_COMM_WORLD, &pnumber);
+
+
+
+
+   x = 0.0 ;
+   y = 0.0 ;
+
+
+   trials = atoi(argv[1]);
+   trials_per_p = trials / (pnumber - 1 );
+
+   x_array = malloc( trials_per_p * sizeof(double)) ;
+   y_array = malloc( trials_per_p * sizeof(double)) ;
+
+
+
+   if(myrank == 0) //master process generates random numbers and sends them to workers
+   {
+      for(j=1 ; j < pnumber ; j++)
+      {
+         for(i=0 ; i < trials_per_p ; i++)
+         {
+            generator(&x, &y);
+               x_array[i] = x ;
+               y_array[i] = y ;
+         }
+
+         MPI_Send(x_array, trials_per_p , MPI_DOUBLE, j, 1234 , MPI_COMM_WORLD) ; 
+         MPI_Send(y_array, trials_per_p , MPI_DOUBLE, j, 5678 , MPI_COMM_WORLD) ;    
+
+      }
+   }
+
+
+   if(myrank != 0) // worker processes receive random numbers from master process and calculate sum
+   {
+
+      MPI_Recv(x_array, trials_per_p, MPI_DOUBLE, 0, 1234, MPI_COMM_WORLD, &status) ;
+      MPI_Recv(y_array, trials_per_p, MPI_DOUBLE, 0, 5678, MPI_COMM_WORLD, &status) ;
+
+      sum = 0.0 ;
+      for(i=0 ; i < trials_per_p ; i++)
+      {
+
+         sum += f(x_array[i], y_array[i]) ;
+      } 
+
+
+   }
+
+   MPI_Reduce(&sum, &total_sum, 1, MPI_DOUBLE , MPI_SUM, 0, MPI_COMM_WORLD) ; // all calculated sums are accumulated in master process
+
+   
+   if(myrank == 0)
+   {
+      integral = total_sum / ( (pnumber -1) * trials_per_p ) ;
+      printf("%f %f\n", integral , M_PI);
+   }
+
+
+   MPI_Finalize();
+
+   return 0;
+
+}

integral_advanced/integral1_openmp.c

+// calculate integral using master-workers model
+// compile: gcc   -fopenmp integral1_openmp.c -lm
+// usage: ./a.out  trial# thread#
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <omp.h>
+
+
+// w is weight function
+double w(double x, double y)
+{
+  return exp( - ( x*x + y*y ) );
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.5 ;
+
+
+   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 ;
+   }
+
+
+}
+
+
+//f is the function to be integrated
+double f(double x, double y)
+{
+
+
+   return M_PI *  ( x*x + y*y ) ;
+
+}
+
+
+int main(int argc, char* argv[])
+{
+   double  x, y ;
+   double * x_array ;
+   double * y_array ;
+
+   double sum = 0.0 ;
+
+   int  count = atoi(argv[1]); // how many random number will generated 
+
+   int tnumber = atoi(argv[2]); // tnumber is total thread number, myrank is rank of each thread
+
+   x_array = malloc(sizeof(double) * count);
+   y_array = malloc(sizeof(double) * count);
+
+   int ready_thread[tnumber] ; //this array holds which thread is ready to run
+
+   #pragma omp parallel num_threads(tnumber) reduction ( +:sum ) 
+   {
+      int myrank = omp_get_thread_num();
+      int i, j ;
+      int number_for_each = count / (tnumber-1) ; // how many random number is used in each thread, master thread is not counted
+
+      
+      if(myrank == 0) //if thread no=0 then generate random numbers, else use generated random numbers
+      {
+
+         x = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1
+         y = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1      
+
+         for(j=1 ; j < tnumber ; j++)
+         {
+            for(i=0 ; i < count ; i++ )
+            {
+               x_array[i] = x ;
+               y_array[i] = y ;
+               generator(&x, &y);
+            }
+            ready_thread[ j ] = 1;
+
+         }
+      }
+      else
+      {
+         int start_index = (myrank - 1) * number_for_each;
+
+         while(ready_thread[myrank] != 1)  
+         {
+            //hold thread until ready bit becomes 1
+         }
+         
+         for(i=0 ; i < number_for_each ; i++ )
+         {
+            sum += f( x_array[start_index+i] , y_array[start_index+i] ) ;
+         }
+
+      }
+      
+
+      #pragma omp barrier
+   } //end pragma
+
+   printf("%f\n", sum / count );
+
+   return 0;
+}

integral_advanced/integral2_mpi.c

+// calculate integral using client-server model
+// compile:  mpicc integral2_mpi.c -lm
+// usage: mpirun -n process# ./a.out trial#
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <mpi.h>
+
+
+#define CHUNK_SIZE 1000
+
+// w is weight function
+double w(double x, double y)
+{
+  return exp( - ( x*x + y*y ) );
+}
+
+//f is the function to be integrated
+double f(double x, double y)
+{
+
+   return M_PI *  ( x*x + y*y ) ;
+
+}
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.5 ;
+
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+   xt = (*x) + delta * (2 * tmp - 1) ; 
+
+   tmp = (double)rand() / (double)RAND_MAX  ;
+   yt = (*y) + delta * (2 * tmp - 1) ; 
+
+   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 sum ;
+   double total_sum ;
+   double integral ;
+   int trials, trials_per_p;
+   int i, j ;
+   int pnumber, myrank;
+   int x_array_tag = 1111 ;
+   int y_array_tag = 2222 ;
+   int demand_tag = 3333 ;
+   int all_done_tag = 4444 ;
+   int sum_tag = 5555 ;
+
+
+   MPI_Init(&argc,&argv);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+   MPI_Comm_size(MPI_COMM_WORLD, &pnumber);
+
+
+
+
+   x = 0.0 ;
+   y = 0.0 ;
+
+
+   trials = atoi(argv[1]);
+
+
+   if(myrank == (pnumber - 1 ))
+   {
+
+      int sent = 0 ;
+      int demander ;
+      int destination = 0;
+      double * x_array ;
+      double * y_array ;
+      int flag=0;
+      MPI_Status status;
+      MPI_Request request;
+      x_array = malloc( trials * sizeof(double)) ;
+      y_array = malloc( trials * sizeof(double)) ;
+
+      int generated = 0 ;
+
+      for(i=0 ;  i < trials ; i++)
+      {
+         generator(&x, &y);
+         x_array[generated] = x ;
+         y_array[generated] = y ;
+         generated++ ;
+      }
+
+
+      while(sent < trials)
+      {
+     
+         MPI_Recv(&demander, 1, MPI_INT, MPI_ANY_SOURCE,  demand_tag, MPI_COMM_WORLD, &status) ;
+
+            MPI_Send(&x_array[sent  ], CHUNK_SIZE , MPI_DOUBLE, demander, x_array_tag, MPI_COMM_WORLD);
+            MPI_Send(&y_array[sent  ], CHUNK_SIZE , MPI_DOUBLE, demander, y_array_tag, MPI_COMM_WORLD);
+            sent += CHUNK_SIZE ;
+
+      }
+
+
+
+
+      for (i=0 ; i < (pnumber-2) ; i++)
+      {
+         MPI_Recv(&demander, 1, MPI_INT, MPI_ANY_SOURCE,  demand_tag, MPI_COMM_WORLD, &status) ;
+
+      }
+
+      double all_done = 1.0 ;
+     
+
+
+      MPI_Send(&all_done, 1 , MPI_DOUBLE, 0, all_done_tag, MPI_COMM_WORLD);
+
+
+   }
+
+
+   if(myrank == 0)
+   {
+      int all_done = 0;
+      double sum = 0.0 ;
+      double temp;
+      MPI_Status status;
+      MPI_Request request;
+
+      while(1)
+      {
+         MPI_Recv(&temp, 1, MPI_DOUBLE, MPI_ANY_SOURCE,  MPI_ANY_TAG , MPI_COMM_WORLD, &status) ;
+         if(status.MPI_SOURCE == (pnumber - 1))
+         {
+            printf("all done received\n");
+            break ;
+         }
+         sum += temp ;
+      }
+
+      printf("%f\n", sum / trials);
+      MPI_Finalize();
+   }
+
+
+
+   
+   if((myrank != 0) && (myrank != pnumber-1) )
+   {
+
+      MPI_Status status;
+      MPI_Request request;
+
+      int flag;
+      double * local_x_array ;
+      double * local_y_array ;
+
+      local_x_array = malloc( CHUNK_SIZE * sizeof(double)) ;
+      local_y_array = malloc( CHUNK_SIZE * sizeof(double)) ;
+
+      while(1)
+      {
+         MPI_Send(&myrank, 1 , MPI_INT, (pnumber-1), demand_tag, MPI_COMM_WORLD);
+         MPI_Recv(local_x_array, CHUNK_SIZE , MPI_DOUBLE, (pnumber - 1 ),  x_array_tag, MPI_COMM_WORLD, &status) ;
+         MPI_Recv(local_y_array, CHUNK_SIZE , MPI_DOUBLE, (pnumber - 1 ),  y_array_tag, MPI_COMM_WORLD, &status) ;
+
+     
+            for(i=0 ; i < CHUNK_SIZE ; i++)
+            {
+               sum += f(local_x_array[i], local_y_array[i]) ;
+
+            }
+         MPI_Send(&sum, 1 , MPI_DOUBLE, 0, sum_tag, MPI_COMM_WORLD);
+         sum = 0.0 ;
+
+       }
+
+         
+
+   }
+
+
+   return 0;
+
+}

integral_advanced/integral2_openmp.c

+// calculate integral using client-server model
+// compile: gcc   -fopenmp integral2_openmp.c -lm
+// usage: ./a.out  trial# thread#
+// note: trial# must be multiple of CHUNK_SIZE
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <omp.h>
+
+#define CHUNK_SIZE 1000
+
+// w is weight function
+double w(double x, double y)
+{
+  return exp( - ( x*x + y*y ) );
+}
+
+
+
+void generator(double * x, double * y)
+{
+   double xt, yt, ratio, tmp ;
+   double delta = 0.5 ;
+
+
+   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 ;
+   }
+
+
+}
+
+
+//f is the function to be integrated
+double f(double x, double y)
+{
+
+
+   return M_PI *  ( x*x + y*y ) ;
+
+}
+
+
+int main(int argc, char* argv[])
+{
+   double  x, y ;
+   double * x_array ;
+   double * y_array ;
+
+   double sum = 0.0 ;
+
+   int  count = atoi(argv[1]); // how many random number will generated 
+
+   int tnumber = atoi(argv[2]); // tnumber is total thread number, myrank is rank of each thread
+
+   x_array = malloc(sizeof(double) * count);
+   y_array = malloc(sizeof(double) * count);
+
+
+   int generated = 0 ;
+   int consumed = 0 ;