Корректные прагмы OpenMP для pi monte carlo в C с генератором случайных чисел, не зависящим от потока

Question

Мне нужна помощь для распараллеливания вычисления pi методом monte carlo с помощью openmp с помощью генератора случайных чисел, который не является потокобезопасным.

Первый: This SO нитка не помогла мне.

Моя собственная попытка - это следующие операторы #pragma omp. Я думал, что i, x и y vars должны быть инициализированы каждым потоком и должны быть закрытыми. z ist сумма всех попаданий в круге, поэтому ее следует суммировать после предполагаемого барриера после цикла for.

Подумайте, основная проблема ist статическое состояние var генератора случайных чисел. Я сделал критический раздел, где вызываются функции, так что только один поток за время может его выполнить. Но решения Pi не масштабируются с более высокими значениями.

Примечание: Я не должен использовать другой RNG, но его все в порядке, чтобы сделать небольшие изменения на нем.

int main (int argc, char *argv[]) { 

    int i, z = 0, threads = 8, iters = 100000; 
    double x,y, pi; 

    #pragma omp parallel firstprivate(i,x,y) reduction(+:z) num_threads(threads) 
     for (i=0; i<iters; ++i) { 
      #pragma omp critical 
      { 
       x = rng_doub(1.0); 
       y = rng_doub(1.0); 
      } 
      if ((x*x+y*y) <= 1.0) 
       z++; 
     } 

    pi = ((double) z/(double) (iters*threads))*4.0; 
    printf("Pi: %lf\n", pi);; 
    return 0; 
} 

Это ГСЧ фактически включаемый файл, но я не уверен, если я создаю файл заголовка правильно, я включил его в другой файл программы, так что у меня есть только один файл .c.

#define RNG_MOD 741025 

int rng_int(void) { 
    static int state = 0; 

    return (state = (1366 * state + 150889) % RNG_MOD); 
} 

double rng_doub(double range) { 
    return ((double) rng_int())/(double) ((RNG_MOD - 1)/range); 
} 

Я также попытался сделать статическое int state глобальным, но это не меняет моего результата, возможно, я сделал это неправильно. Так что, пожалуйста, не могли бы вы помочь мне внести правильные изменения? Большое спасибо!

Answer 1

Пробуйте использовать приведенный ниже код. Он создает частное состояние для каждого потока. Я сделал что-то подобное с функцией rand_rWhy does calculation with OpenMP take 100x more time than with a single thread?

Редактировать: Я обновил свой код, используя некоторые предложения Христи. Я использовал threadprivate (впервые). Я также использовал лучшую функцию rand, которая дает лучшую оценку pi, но она все еще недостаточно хороша.

Одна странная вещь, которую я должен был определить функцию rng_int после threadprivate иначе я получил ошибку «ошибка:« состояние »объявлено« threadprivate »после первого использования». Вероятно, я должен задать вопрос об этом.

//gcc -O3 -Wall -pedantic -fopenmp main.c 
#include <omp.h> 
#include <stdio.h> 

#define RNG_MOD 0x80000000 
int state; 

int rng_int(void); 
double rng_doub(double range); 

int main() { 
    int i, numIn, n; 
    double x, y, pi; 

    n = 1<<30; 
    numIn = 0; 

    #pragma omp threadprivate(state) 
    #pragma omp parallel private(x, y) reduction(+:numIn) 
    { 

     state = 25234 + 17 * omp_get_thread_num(); 
     #pragma omp for 
     for (i = 0; i <= n; i++) { 
      x = (double)rng_doub(1.0); 
      y = (double)rng_doub(1.0); 
      if (x*x + y*y <= 1) numIn++; 
     } 
    } 
    pi = 4.*numIn/n; 
    printf("asdf pi %f\n", pi); 
    return 0; 
} 

int rng_int(void) { 
    // & 0x7fffffff is equivalent to modulo with RNG_MOD = 2^31 
    return (state = (state * 1103515245 + 12345) & 0x7fffffff); 
} 

double rng_doub(double range) { 
    return ((double)rng_int())/(((double)RNG_MOD)/range); 
} 

Вы можете увидеть результаты (а также редактировать и запускать код) на http://coliru.stacked-crooked.com/a/23c1753a1b7d1b0d

Answer 2

Ваш оригинальный линейный конгруэнтны ПСЧ имеет длину цикла 49400, поэтому вы только получаете 29700 уникальные тестовые точки. Это ужасный генератор, который можно использовать для любого моделирования Монте-Карло. Даже если вы сделаете 100000000 испытаний, вы не приблизитесь к истинному значению Pi, потому что вы просто повторяете одни и те же точки снова и снова, и в результате как конечное значение z, так и iters просто умножаются на одни и те же константа, которая отменит в конце во время деления.

Семена для каждой нити, введенные Z бозоном, немного улучшают ситуацию с количеством уникальных точек, увеличивающихся с общим количеством потоков OpenMP. Увеличение не является линейным, так как если семя одного PRNG попадает в последовательность другого PRNG, оба PRNG производят одну и ту же последовательность, сдвинутую не более чем на 49400 элементов. Учитывая длину цикла, каждый PRNG охватывает 49400/RNG_MOD = 6,7% от общего диапазона выпуска, и это вероятность синхронизации двух PRNG. Всего существует RNG_MOD/49400 = 15 уникальных последовательностей.Это в основном означает, что в наилучшем сценарии посева вы не сможете пройти мимо 30 потоков, поскольку любой другой поток просто повторит результат некоторых других. Множитель 2 исходит из того, что каждая точка использует два элемента из последовательности и, следовательно, можно получить другой набор точек, если вы сдвигаете последовательность на один элемент.

Окончательное решение является полностью уроните ПСЧ и придерживаться что-то вроде Mersenne twister MT19937, который имеет длину цикла 2 - 1 и очень сильный алгоритм высева. Если вы не можете использовать другой ПГСЧ как Вы заявляете в вашем вопросе, по крайней мере, изменить константы LCG, чтобы соответствовать тем, которые используются в rand():

int rng_int(void) { 
    static int state = 1; 

    // & 0x7fffffff is equivalent to modulo with RNG_MOD = 2^31 
    return (state = (state * 1103515245 + 12345) & 0x7fffffff); 
} 

Обратите внимание, что rand() не является хорошей ГСЧ - это еще Плохо. Это немного лучше, чем тот, который используется в вашем коде.