Какие функции и операции разрешены в параллельном блоке?

Генератор случайных чисел rand() использует глобальное статически распределенное состояние, разделяемое всеми потоками, и поэтому не является потокобезопасным. Используя его из нескольких потоков, вы столкнетесь с очень серьезным случаем незащищенного доступа к общим переменным, который засоряет кеш и замедляет работу программы. Вместо этого вы должны использовать rand_r() или erand48() - они используют отдельные хранилища состояний, которые вы должны предоставить. Вы должны объявить одно состояние для каждого потока (например, иметь его private), в основном создавая разные PRNG для каждого потока. Затем вы должны раздать их соответствующим образом, иначе вы получите статистически плохие результаты. В принципе, вы можете использовать выходные данные одного генератора rand48() для заполнения других - этого должно быть достаточно, чтобы получить умеренно длинные некоррелированные последовательности.

Вот пример реализации с использованием rand_r() (это не значит, что это очень плохой генератор для моделирования методом Монте-Карло, erand48 лучше, и лучше всего использовать генератор типа "Mersenne Twister" из научной библиотеки GNU, если доступный):

unsigned int prng_state;
#pragma omp threadprivate(prng_state)

double x(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double y(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double z(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}

int d(double x, double y, double z){
    if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
    return 0;
}

double f(double x, double y, double z){
    return 1;
}

...

#pragma omp parallel default(none) \
            private(id,numt,x_,y_,z_) \
            shared(count_of_iterations) \
            reduction(+:total_sum,good_dots)
{
    id = omp_get_thread_num() + 1;
    numt = omp_get_num_threads();

    // Sample PRNG seeding code - DO NOT USE IN PRODUCTION CODE!
    prng_state = 67894 + 1337*id;

    #pragma omp for
    for (j = 1; j <= count_of_iterations;  j++){
        x_=x();
        y_=y();
        z_=z();
        if (d(x_,y_,z_) == 1){
            total_sum += f(x_,y_,z_);
            good_dots += 1;
        }
    }
}

Это просто очень плохая (с точки зрения качества) реализация, но она должна дать вам представление о том, как все работает. Это также то, как вы можете достичь безопасности потоков с минимальными изменениями в исходном коде. Основные моменты:

• состояние PRNG prng_state делается приватным для каждого потока директивой OpenMP threadprivate;
• rand_r() с переменной состояния, зависящей от потока, используется вместо rand() в x(), y() и z();
• состояние PRNG инициализируется в зависимости от потока, например. prng_state = 67894 + 1337*id;, чтобы разные потоки (надеюсь) получали некоррелированные потоки псевдослучайных чисел.

Обратите внимание, что rand() и rand_r() настолько плохого качества, что это просто академический пример. С более длинными последовательностями PRNG вы получите коррелированные потоки в разных потоках, что испортит статистику. Я оставляю на ваше усмотрение переписать код, используя erand48().

Чтобы ответить на ваш первоначальный вопрос, все потокобезопасные вызовы функций разрешены внутри блока parallel. Вы также можете вызывать не потокобезопасные функции, но вы должны защитить вызовы внутри (именованных) конструкций critical, например:

#pragma omp for
for (j = 1; j <= local_coi; j++) {
    #pragma omp critical(prng)
    {
        x_=x();
        y_=y();
        z_=z();
    }
    if (d(x_,y_,z_) == 1) {
        local_sum += f(x_,y_,z_);
        local_good_dots += 1;
    }
}

Это гарантирует, что параллельные вызовы rand() не будут выполняться. Но у вас по-прежнему будет доступ для чтения-изменения-записи к общему состоянию, отсюда и замедление, связанное с кешем.

Кроме того, не пытайтесь повторно реализовать OpenMP reduction или подобные конструкции. Поставщики компиляторов уже прилагают огромные усилия, чтобы убедиться, что они реализованы наилучшим (читай, самым быстрым) способом.