Вычисление пи с помощью Монте-Карло с использованием OpenMP

Я новичок в C и изучаю параллелизм с помощью C. Я наткнулся на упражнение в книга с просьбой найти приблизительное значение Pi, используя технику Монте-Карло с OpenMP. Я придумал следующее:

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <math.h>
#include <omp.h>

int in_circle = 0;
int total = 10000;
int thread_numb = 4;

void calculate_in_circle_random_count() {
  srand((unsigned int)time(NULL));
  for (int i = 0; i < total / thread_numb; i++) {
    float x = (float)rand()/(float)(RAND_MAX);
    float y = (float)rand()/(float)(RAND_MAX);

    float val = (x * x) + (y * y);

    if (val < 1.0) {
      #pragma omp critical
      {
        in_circle++;
      }
    }
  }
}

int main(int argc, char const *argv[]) {
  #pragma omp parallel num_threads(thread_numb)
  {
    printf("Hello from process: %dn", omp_get_thread_num());
    calculate_in_circle_random_count();
  }

  float pi_approx = 4 * (float)in_circle / (float)total;

  printf("In circle: %dn", in_circle);
  printf("Total: %dn", total);
  printf("Pi approximation: %.6fn", pi_approx);

  return 0;
}

У меня есть выборка 10000 для расчета приближения, и я хочу выделить 4 потока, которые будут параллельны для расчета. Так как мне нужно 4 потока, я убеждаюсь, что цикл for выполняется до тех пор, пока total / 4 для каждого потока. И всякий раз, когда (val < 1.0) является true Я увеличиваю in_circle переменная в критическом сечении.

Имеет ли смысл такой подход? Если нет, то как это можно улучшить?

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

Чтобы исправить это, вам в значительной степени нужно использовать другой генератор случайных чисел. Ваша реализация может предоставить один (например, erand48 довольно часто). Если вам действительно нужно, чтобы ваш код был переносимым, вы можете написать свой собственный, что-нибудь в таком порядке:

#include <time.h>

typedef unsigned long long rand_state;

// multiplier/modulus taken from Knuth Volume 2, Table 1
static const int multiplier = 314159269;
static const int addend = 1;
static const int modulus = 0xffffffff;

// note that this works differently from srand, returning a seed rather than setting
// a hidden seed.
rand_state omp_srand() {
    rand_state state = time(NULL);
    state ^= (unsigned long long)omp_get_thread_num() << 32;
    return state;
}

int omp_rand(rand_state *state) {
    *state = *state * multiplier + addend;
    return *state & modulus;
}

Обратите внимание: поскольку вполне вероятно, что все потоки запускаются в одну и ту же секунду, это объединяет текущее время с идентификатором потока для заполнения генератора каждого потока. Хотя это не на 100% гарантировано, это дает довольно высокую вероятность того, что генератор каждого потока будет запускаться с уникальным начальным значением. С другой стороны, это также означает, что мы должны запустить поток, а затем заполнить генератор внутри этого потока.

Чтобы использовать это, наш код будет примерно таким:

double calculate_in_circle_random_count(void) {

    static const unsigned total = 1000000000;
    unsigned in_circle = 0;

#pragma omp parallel reduction(+:in_circle)
    {
        // do per-thread initialization
        rand_state seed = omp_srand();
        int count = total / omp_get_num_threads();
        int i;

        // then do this thread's portion of the computation:
        for (i = 0;  i < count;  ++i) {

            double x = (double)omp_rand(&seed) / OMP_RAND_MAX;
            double y = (double)omp_rand(&seed) / OMP_RAND_MAX;

            double val = x * x + y * y;
            in_circle +=  val < 1.0;
        }
    }
    return 4.0 * in_circle / total;
}

int main(int argc, char const *argv[])
{
    int num_threads = 1;

    if (argc > 1) 
        num_threads = atoi(argv[1]);

    omp_set_num_threads(num_threads);

    float pi_approx = calculate_in_circle_random_count();

    printf("Pi approximation: %.6fn", pi_approx);

    return 0;
}

С этой модификацией код как минимум способный масштабирования. Чтобы он хорошо масштабировался, вам нужно изменить total к гораздо большему числу — с таким маленьким, как вы указали, для запуска нескольких потоков требуется больше времени, чем это экономит при расчетах. Но хотя бы этот код может хорошо масштабируется, поэтому если мы сделаем total совсем немного больше, он действительно будет работать быстрее. Для тайминга я добавил еще несколько нулей к total, и переписал main немного, чтобы создать цикл для использования 1, 2, 3 и 4 потоков и распечатать время каждой итерации. На моей машине это произвело следующее:

With 1 threads: Pi approximation: 3.140259, time: 3,652,683 microseconds
With 2 threads: Pi approximation: 3.139389, time: 1,892,415 microseconds
With 3 threads: Pi approximation: 3.138885, time: 1,268,917 microseconds
With 4 threads: Pi approximation: 3.138306, time:   935,579 microseconds

Таким образом, 4 потока как минимум почти в 4 раза быстрее, чем один поток. Я почти уверен, что если вы используете rand внутри цикла у вас никогда не получится хорошо масштабироваться.

Для всех, кому не все равно, версия main, которая производит этот вывод, написана на C ++ и выглядит так:

int main(int argc, char const* argv[]) {
    using namespace std::chrono;

    std::cout.imbue(std::locale(""));

    int processors = omp_get_num_procs();

    for (int num_threads = 1; num_threads <= processors; num_threads++) {
        std::cout << "With " << std::setw(2) << num_threads << " threads: ";
        omp_set_num_threads(num_threads);

        auto start = high_resolution_clock::now();
        float pi_approx = calculate_in_circle_random_count();
        auto stop = high_resolution_clock::now();

        printf(" %.6f, time: ", pi_approx);
        std::cout << std::setw(9) << duration_cast<microseconds>(stop - start).count() << " microsecondsn";
    }
    return 0;
}

Это остается почти таким же, но с некоторым взломанным кодом времени. В частности, это нет попытка переписать код на C ++ в целом. Если бы я делал это, я бы, вероятно, сделал ряд вещей по-другому (начиная с того факта, что C ++ <random> header уже предоставляет чистый способ обработки генерации случайных чисел для каждого потока, поэтому я бы использовал его вместо того, чтобы накатывать свои собственные).

Я думаю, вы слишком усложняете это. Вы можете позволить OpenMP разделить работу между потоками за вас, используя #pragma openmp parallel for.

Вы можете использовать сокращение для in_circle вместо критического раздела — это позволяет OpenMP суммировать переменную для каждого потока и добавлять их все в конце, уменьшая конкуренцию за переменную.

Вот более простая версия (я также удалил неиспользуемые включения и сделал объявление функции прототипом):

#include <stdio.h>
#include <stdlib.h>

double calculate_in_circle_random_count(void)
{
    static const unsigned total = 1000000;
    unsigned in_circle = 0;

#pragma omp parallel for reduction(+:in_circle)
    for (unsigned i = 0;  i < total;  ++i) {
        double x = (double)rand() / RAND_MAX;
        double y = (double)rand() / RAND_MAX;

        double val = x * x + y * y;
        in_circle +=  val < 1.0;
    }

    return 4.0 * in_circle / total;
}

int main(void)
{
    const double pi_approx = calculate_in_circle_random_count();
    printf("Pi approximation: %.6fn", pi_approx);
}

Обратите внимание, что в этом масштабе преобладают накладные расходы на распараллеливание, что делает параллельную версию намного медленнее, чем последовательную. Чем больше итераций, тем ближе вы приблизитесь, но в цикле так мало, что пользы невелики.