Uniform and Poisson distributions: Generating random numbers in parallel mode, where "bad" is using just one common generator, and "good" is using undependent generator for each thread.

poisson_2d.cpp

#include <assert.h>
#include <chrono>
#include <iostream>
#include <omp.h>
#include <random>
#include <vector>

using namespace std;

#define Nx 50
#define Ny 25
#define ITERATIONS (Nx * Ny * Nx * Ny)

// must be divisible by number of processor cores or less
#define THREAD_NUMS 3

int index(int x, int y)
{
    return y * Nx + x;
}

void show(int *values)
{
    const int cn = 1024;
    char symbols[cn] =
        " .-~+=ox%X#8O0@:::::::::^^^^^^^^^^*****************";

    for (int i = 50; i < 100; ++i)
        symbols[i] = '!';

    for (int i = 100; i < 350; ++i)
        symbols[i] = '?';

    for (int i = 350; i < cn; ++i)
        symbols[i] = '$';


    for (int y = 0; y < Ny; ++y)
    {
        cout << "// ";
        for (int x = 0; x < Nx; ++x)
        {
            int value = (double)values[index(x, y)] * Nx * Ny /
                ((double)ITERATIONS / THREAD_NUMS);
            assert(value < cn);
            cout << symbols[value];
        }
        cout << "\n";
    }
    cout << endl;
}

int main()
{
    omp_set_num_threads(THREAD_NUMS);

    cout.precision(2);

    int values[Nx * Ny];
    int rs[Nx * Ny];
    int badV[Nx * Ny];
    int badR[Nx * Ny];

    for (int i = 0; i < Nx * Ny; ++i)
    {
        values[i] = 0;
        rs[i] = 0;

        badV[i] = 0;
        badR[i] = 0;
    }

    auto d = std::chrono::high_resolution_clock::now().time_since_epoch();
    std::uniform_int_distribution<int> firstDist(0, d.count());

    int commonRx = 0, commonRy = 0;
    std::mt19937_64 firstGen;
    firstGen.seed(d.count());

    int badRx = 0, badRy = 0;
    std::mt19937_64 badGen;
    badGen.seed(d.count());

    std::mt19937_64 _randomGenerator[THREAD_NUMS];

#pragma omp parallel
    {
#pragma omp critical
        _randomGenerator[omp_get_thread_num()].seed(firstDist(firstGen));

        int localValues[Nx * Ny];
        int localBad[Nx * Ny];
        for (int i = 0; i < Nx * Ny; ++i)
        {
            localValues[i] = 0;
            localBad[i] = 0;
        }

        for (int i = 0; i < ITERATIONS; ++i)
        {
            std::uniform_int_distribution<int> distX(0, Nx - 1);
            std::uniform_int_distribution<int> distY(0, Ny - 1);

            // good
            int x = distX(_randomGenerator[omp_get_thread_num()]);
            int y = distY(_randomGenerator[omp_get_thread_num()]);

#pragma omp atomic
            commonRx += x;
#pragma omp atomic
            commonRy += y;

            int k = index(x, y);
            assert(k < Nx * Ny);
#pragma omp atomic
            ++localValues[k];

            // bad
            int bx = distX(badGen);
            int by = distY(badGen);

#pragma omp atomic
            badRx += bx;
#pragma omp atomic
            badRy += by;

            int bk = index(bx, by);
            assert(bk < Nx * Ny);
#pragma omp atomic
            ++localBad[bk];

#pragma omp master
            {
                // good
                int rx = commonRx / omp_get_num_threads();
                int ry = commonRy / omp_get_num_threads();

                if (rx >= Nx) rx = Nx - 1;
                if (ry >= Ny) ry = Ny - 1;

                int j = index(rx, ry);
                assert(j < Nx * Ny);

                ++rs[j];

                commonRx = 0;
                commonRy = 0;

                // bad
                int rbx = badRx / omp_get_num_threads();
                int rby = badRy / omp_get_num_threads();

                if (rbx >= Nx) rbx = Nx - 1;
                if (rby >= Ny) rby = Ny - 1;

                j = index(rbx, rby);
                assert(j < Nx * Ny);

                ++badR[j];

                badRx = 0;
                badRy = 0;
            }

#pragma omp barrier
        }

        for (int i = 0; i < Nx * Ny; ++i)
        {
#pragma omp atomic
            values[i] += localValues[i];

#pragma omp atomic
            badV[i] += localBad[i];
        }
    }

    cout << "//   BAD:" << endl;
    show(badV);
    show(badR);

    cout << "\n\n\n//   GOOD:" << endl;
    show(values);
    show(rs);

    return 0;
}


//   BAD:
// ======================o=================
// ========================================
// ========================================
// ========================================
// ===o=============================o==o===
// ========================================
// ==o=====================================
// ========================================
// ========================================
// ================o==o=====o==============
// ======================================o=
// ========================================
// ========================================
// ========================================
// ============================o====o======
// ========================================
// ========================================
// o=======================================
// ========================================
// ========================================

// ........................ ..
// ..........-.................
// ......------------...........
// .....------------------.......
// .....----------------------..... .  .. .
// ....-------~~~-~-~~~~-------........   .
// .....-----~~~~~~~~~~~~~--~-----.........
//  .....-------~~~~~~~~~~~~~----..........
// ......------~~~~~~~~~~~~~~------...... .
//  ......------~~~~~~~~~~~-~------........
// .......------~~~~~~~~~~~~-~-----........
//  .......----~-~~~~~~~~~~~--------.......
// ........------~~~~~~~~~~~~--------......
//   .......----------~~-~~~---------......
//       ......---------~------------......
//         .........----------------.......
//           ...........----------.........
// .        . .............................
//             ............................
//        .    .. ....... .. ......... ...-




 //   GOOD:
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ========================================
// ====================o===================
// ========================================

//
//
//                 .........
//            ......-----......
//         ....---~~~~~~~~----...
//        ...---~~+++++++++~~~--...
//       ..--~~++==ooooooo==++~--....
//      ...--~+==oxxxx%%%xxoo=+~~--..
//      ..--~~+=oxx%%XXX%%xxo==+~~-..
//      ..--~++=ox%%XXXXX%%xoo=+~--..
//      ..--~~==oxx%%XX%X%%xxo=+~~-..
//      ...-~~+=ooxx%%%%%%xoo=++~--..
//       ..--~~+==ooxxxxxooo=+~~--...
//        ..---~~+++======++~~---..
//         ....--~~~~~~~~~~----...
//           .....----------....
//              .............
//
//
//