OpenMP Monte Carlo \(\pi\)

Let’s parallelize our Monte Carlo example of computing pi.

Important

We want each thread to have independent random numbers. We do this by creating a random number generator separately on each thread with different seeds.

See, e.g., https://stackoverflow.com/questions/37305104/how-to-generate-random-numbers-in-a-thread-safe-way-with-openmp

There is a worry still—if you are sampling a large fraction of the period of the pseudo-random number generator, than it is possible that the sequences will overlap. This is discussed in the SC08 presentation by Tim Mattson & Larry Meadows (look around slide 130).

Tip

The Merseene twister has a period of \(2^{19937}-1\), so it is unlikely for overlap to be an issue.

Note

True parallel random number generators exist, where a single thread controls the generator and gives the numbers to each thread as needed.

See also: Improved Long-Period Generators Based on Linear Recurrences Modulo 2 by Panneton and L’Ecuyer.

Then this calculation can proceed in a trivially parallel fashion. The only thing we need to do is mark the sum for a reduction at the end of the parallel region.

Listing 131 pi.cpp

#include <iostream>
#include <iomanip>

#include <random>
#include <cmath>

#ifdef _OPENMP
#include <omp.h>
#endif

const int N = 1000000;

int main() {

#ifdef _OPENMP
    auto start = omp_get_wtime();
#endif

    int N_circle{0};

    #pragma omp parallel reduction(+:N_circle)
    {

        // each thread will have its own random number
        // engine

        std::random_device rd;
#ifdef _OPENMP
        std::seed_seq seed{static_cast<int>(rd()), omp_get_thread_num()};
#else
        std::seed_seq seed{static_cast<int>(rd())};
#endif
        std::mt19937 generator(seed);
        std::uniform_real_distribution<double> uniform(0.0, 1.0);

        #pragma omp for
        for (int n = 0; n < N; ++n) {

            double x = uniform(generator);
            double y = uniform(generator);

            double r2 = std::pow(x - 0.5, 2) + std::pow(y - 0.5, 2);

            if (r2 <= 0.25) {
                N_circle++;
            }

        }

    }

    std::cout << "pi = " << std::setprecision(10)
              << 4 * static_cast<double>(N_circle) / N << std::endl;

#ifdef _OPENMP
    std::cout << "time = " << omp_get_wtime() - start << std::endl;
#endif

}

Tip

There is one further trick here—all of the OpenMP stuff is wrapped in #ifdef _OPENMP—the _OPENMP preprocessor variable is defined by the standard if we are building with OpenMP. By wrapping those features, this allows the code to be compiled without OpenMP if we desire.

Tip

Other types of random number generators exist, beyond those in the C++ SL. Here’s an example of a header-only library that works on CPUs and GPUs: random123.