OpenMP reduce array without thread-local copies

Question

I wish to OpenMP reduce a large array into a smaller dynamic array. For example, where

large[] = {1, 1, 1, 2, 2, 2, 3, 3, 3};

// OpenMP reduce length-3 sublists of large, to produce:

small[] = {3, 6, 9};

My requirements are similar to this question, but I have important additional constraints:

• must support OpenMP 3.1 (so I cannot use OpenMP 4.5's array reductions, as this answer does)
• cannot have a thread-private copy of small (like this answer), since small can be just as large as large, and hence thread-private copies may cause a stack overflow
• large array must be iterated in order, for good caching (since it may be very large)

Some other details:

• the summed elements of large are not necessarily adjacent
• this reduction happens within a function. The small (output) array is pre-allocated by the user, and passed by reference. Ideally, the function should not allocate a new local copy of small
• both small and large have lengths of a power of 2 (e.g. 2, 4, 8, 16 ...). Every element of small is reduced from the same number of large elements (another power of 2).

Here is example serial psuedocode, for clarity:

void myfunc(double* small, double* large, int lLen, // other args) {

    for (int i=0; i<lLen; i++) {
    
        int sInd = // determined by other args
        
        small[sInd] += large[i];
    }
}

Here is an example implementation, disqualified since it makes use of OpenMP 4.5's array reduction. It furthermore undesirably makes use of a local copy of small (as required to use reduction).

void myfunc(double* small, int sLen, double* large, int lLen, // other args) {

    double smallLocal[sLen];
    int i, sInd;

    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        smallLocal[i] = 0;
    
    #pragma omp parallel private (i, sInd) reduction (+:smallLocal) for
    for (i=0; i<largeLen; i++) {

        sInd = // determined by other args
        
        smallLocal[sInd] += large[i];
    }

    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        small[i] = smallLocal[i];
}

I'm hoping I can achieve the same thing in OpenMP 3.1, possibly avoiding even the local copy of smallLocal, by managing the array element reductions myself. How can I do this?

Answer 1

You could use OpenMPs atomic update construct, which is already present in OpenMP 3.1:

void myfunc(double* small, double* large, int lLen, // other args) {

    #pragma omp parallel for
    for (int i=0; i<lLen; i++) {
    
        int sInd = // determined by other args
        #pragma omp atomic update
        small[sInd] += large[i];
    }
}

This should be faster than using locks.

Answer 2

Here's a somewhat unsatisfying solution using locks, which requires a local locks array as big as small (which may need to be malloc'd).

void myfunc(double* small, int sLen, double* large, int lLen, // other args) {

    int i, sInd;

    // clear small
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        small[i] = 0;

    // initialise locks
#ifdef _OPENMP
    omp_lock_t locks[sLen];
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        omp_init_lock(&locks[i]);
#endif 
    
    // main loop
    #pragma omp parallel private (i, sInd) for
    for (i=0; i<largeLen; i++) {

        sInd = // determined by other args
        
        // update small with locks
#ifdef _OPENMP
        omp_set_lock(&locks[sInd]);
        small[sInd] += large[i];
        omp_unset_lock(&locks[sInd]);
#else 
        small[sInd] += large[i];
#endif
    }

    // destroy locks
#ifdef _OPENMP
    #pragma omp parallel private (i) for
    for (i=0; i<sLen; s++)
        omp_destroy_lock(&locks[i]);
#endif 
}

As expected, this was slower than using OpenMP 4.5's array reduction:

• for sLen = 8, lLen = 1048576, using locks was 26x slower
• for sLen = 32768, lLen = 1048576, using locks was 1.5x slower

It should be mentioned that OpenMP 4.5's array reduction is only compatible with local stack arrays, and not with dynamic memory. This limits the size of the reduced array to ~1 MiB, or 2^17 doubles (lLen = 131072). In contrast, this locks solution can handle any size, since locks can be malloc'd.

Answer 3

You could do this sort of thing, which is common in scientific codes:

#pragma omp parallel for
for(int i = 0; i < N/3; i++)
{
    for (int j = 0; j < 3; j++)
    {
        small[i] += large[3*i + j];
    }
}