Fill histograms (array reduction) in parallel with OpenMP without using a critical section

Question

I would like to fill histograms in parallel using OpenMP. I have come up with two different methods of doing this with OpenMP in C/C++.

The first method proccess_data_v1 makes a private histogram variable hist_private for each thread, fills them in prallel, and then sums the private histograms into the shared histogram hist in a critical section.

The second method proccess_data_v2 makes a shared array of histograms with array size equal to the number of threads, fills this array in parallel, and then sums the shared histogram hist in parallel.

The second method seems superior to me since it avoids a critical section and sums the histograms in parallel. However, it requires knowing the number of threads and calling omp_get_thread_num(). I generally try to avoid this. Is there better way to do the second method without referencing the thread numbers and using a shared array with size equal to the number of threads?

void proccess_data_v1(float *data, int *hist, const int n, const int nbins, float max) {
    #pragma omp parallel 
    {
        int *hist_private = new int[nbins];
        for(int i=0; i<nbins; i++) hist_private[i] = 0;
        #pragma omp for nowait
        for(int i=0; i<n; i++) {
            float x = reconstruct_data(data[i]);
            fill_hist(hist_private, nbins, max, x);
        }
        #pragma omp critical 
        {
            for(int i=0; i<nbins; i++) {
                hist[i] += hist_private[i];
            }
        }
        delete[] hist_private;
    }
}

void proccess_data_v2(float *data, int *hist, const int n, const int nbins, float max) {
    const int nthreads = 8;
    omp_set_num_threads(nthreads);
    int *hista = new int[nbins*nthreads];

    #pragma omp parallel 
    {
        const int ithread = omp_get_thread_num();
        for(int i=0; i<nbins; i++) hista[nbins*ithread+i] = 0;
        #pragma omp for
        for(int i=0; i<n; i++) {
            float x = reconstruct_data(data[i]);
            fill_hist(&hista[nbins*ithread], nbins, max, x);
        }

        #pragma omp for
        for(int i=0; i<nbins; i++) {
            for(int t=0; t<nthreads; t++) {
                hist[i] += hista[nbins*t + i];
            }
        }

    }
    delete[] hista;
}

Edit: Based on a suggestion by @HristoIliev I have created an improved method called process_data_v3

#define ROUND_DOWN(x, s) ((x) & ~((s)-1))
void proccess_data_v2(float *data, int *hist, const int n, const int nbins, float max) {
    int* hista;
    #pragma omp parallel 
    {
        const int nthreads = omp_get_num_threads();
        const int ithread = omp_get_thread_num();

        int lda = ROUND_DOWN(nbins+1023, 1024);  //1024 ints = 4096 bytes -> round to a multiple of page size
        #pragma omp single
        hista = (int*)_mm_malloc(lda*sizeof(int)*nthreads, 4096);  //align memory to page size

        for(int i=0; i<nbins; i++) hista[lda*ithread+i] = 0;
        #pragma omp for
        for(int i=0; i<n; i++) {
            float x = reconstruct_data(data[i]);
            fill_hist(&hista[lda*ithread], nbins, max, x);
        }

        #pragma omp for
        for(int i=0; i<nbins; i++) {
            for(int t=0; t<nthreads; t++) {
                hist[i] += hista[lda*t + i];
            }
        }

    }
    _mm_free(hista);
}

Could you please explain why you are using nested parallel regions? (I am referring to your process_data_v1 approach). Maybe I am not understanding something, but according to your code, it seems to me that you are asking for Nthreads**2. It is to say, you are asking for more resources than the available ones. Is that correct? In other words, could you explain the behaviour of parallel regions inside the outer one? Thanks... — Alejandro
Hi @user2088790 , isn't proccess_data_v1 the fastest one? Because we don't need shared memory. I try version2 and 3, they are slower than v1. Any suggestion? — Ardian

Hristo Iliev Hristo Iliev · Accepted Answer · 2013-05-28T12:05:37

You could allocate the big array inside the parallel region, where you can query about the actual number of threads being used:

int *hista;
#pragma omp parallel 
{
    const int nthreads = omp_get_num_threads();
    const int ithread = omp_get_thread_num();

    #pragma omp single
    hista = new int[nbins*nthreads];

    ...
}
delete[] hista;

For better performance I would advise that you round the size of each thread's chunk in hista to a multiple of the system's memory page size, even if this could potentially leave holes between the different partial histograms. This way you will prevent both false sharing and remote memory access on NUMA systems (but not in the final reduction phase).

Fill histograms (array reduction) in parallel with OpenMP without using a critical section

1 Answers