vendredi 31 janvier 2020

Random values generator bug with OpenMP tasks

I'm trying to develop a micro code based on the calculation of Pi by the Monte Carlo method while using OpenMP tasks. It seems that I have some troubles in that code. Indeed, my application simply crashes and I got a segmentation fault coming from the gomp_barrier function...

I don't know if I am doing wrong things in my code or not, but I spend a little bit time on time to try to debug it and I found nothing. When I simply remove the computation of x and y with random values generator, the code is doing well... But if I add a call to rand_r or srand48_r, I have this segfault. So maybe, my use of these functions is wrong. Does anyone have an idea please ?

Here is my code :

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <inttypes.h>
#include <time.h>
#include <unistd.h>

#include "omp.h"

#define TRIALS_PER_THREAD 10E6

int main(int argc, char** argv)
{
  uint64_t const n_test = TRIALS_PER_THREAD;
  uint64_t i;
  double pi = 0.;

  int nb_threads = 2;
  #pragma omp parallel shared(nb_threads)
  {
    #pragma omp master
      nb_threads = omp_get_num_threads();

  }
  fprintf(stdout, "Nb threads: %d\n", nb_threads);

  uint64_t* result = (uint64_t*)malloc(sizeof(uint64_t) * nb_threads);
  for(i = 0; i < nb_threads; ++i)
  {
    result[i] = 0;
  }

  int nb_test_per_thread = 20;

  #pragma omp parallel\
  shared(result)\
  firstprivate(n_test,nb_test_per_thread)
  {
    unsigned int seed;
    struct drand48_data randBuffer;

    seed = time(NULL) ^ omp_get_thread_num() ^ getpid();
    srand48_r(seed, &randBuffer);

    for(int k = 0; k < nb_test_per_thread; ++k)
    {
      #pragma omp task shared(result, seed, randBuffer)
      {
        uint64_t local_res = 0;
        double x = 0., y = 0.;
        for(i = 0; i < n_test; ++i)
        {
          drand48_r(&randBuffer, &x);// / (double)RAND_MAX;
          drand48_r(&randBuffer, &y);// / (double)RAND_MAX;

          local_res += (((x * x) + (y * y)) <= 1);
        }
        int tid = omp_get_thread_num();
        result[tid] += local_res;
      }
    }
  }

  for(i = 0; i < nb_threads; ++i)
  {
    pi += result[i];
  }

  fprintf(stdout, "%ld of %ld throws are in the circle !\n", (uint64_t)pi, n_test*nb_test_per_thread*nb_threads);
  pi *= 4;
  pi /= (double)(n_test*nb_test_per_thread*nb_threads);
  fprintf(stdout, "Pi ~= %f\n", pi);

  return 0;
}

Thank you in advance for your help !




Aucun commentaire:

Enregistrer un commentaire