OpenMP with parallel reduction in for loop

Question

I have a for-loop to iterate over a rather large amount of points (ca. 20000), for every point it is checked whether or not the point is inside some cylinder (that cylinder is the same for every point). Furthermore, I would like to have the highest Y coordinate from the set of points. Since I have to do this calculation a lot, and it's quite slow, I want to use OpenMP to parallelize the loop. Currently I have (somewhat reduced):

#pragma omp parallel for default(shared) private(reducedCloudSize, reducedCloud, cylinderBottom, cylinderTop) reduction(+:ptsInside, ptsInsideLarger)
for (int i = 0; i < reducedCloudSize; i++){
    highestYCoord = highestYCoord > testPt.y ? highestYCoord : testPt.y;

    if (CylTest_CapsFirst(cylinderBottom,cylinderTop,cylinderHeight*cylinderHeight,cylinderRadius*cylinderRadius,testPt) != -1){
        ptsInside++;
    }

}

Where the CylTest_CapsFirst will check whether the point is inside of the cylinder. However, this code does not work. If I leave out the reduction(+:ptsInside, ptsInsideLarger) part it actually works, but is much slower than the non-parallelized version. If I include the reduction clause, the program never even seems to enter the for-loop!

What am I doing wrong?

Thanks!

Answer 1

Assuming your function CylTest_CapsFirst does not write to anything (only reads) the only variables that need to be shared are highestYCoord and ptsInside. The only variable that needs to be private is i. You don't need to explictly declare these. But you do need to make sure that no thread writes to the shared variables at the same time. To do this efficiently you should make private versions of highestYCoord and ptsInside which you write in the parallel loop. Then you can merge the private versions with the shared version in a critical section. This is efficient as long as reducedCloudSize >> number_of_threads.

#pragma omp parallel
{
    double highestYCoord_private = highestYCoord;
    int ptsInside_private = 0;
    #pragma omp for
    for (int i = 0; i < reducedCloudSize; i++){
        highestYCoord_private = highestYCoord_private > testPt.y ? highestYCoord_private : testPt.y;
        if (CylTest_CapsFirst(cylinderBottom,cylinderTop,cylinderHeight*cylinderHeight,cylinderRadius*cylinderRadius,testPt) != -1) {
                ptsInside_private++;
        }
    }
    #pragma omp critical 
    {
        highestYCoord = highestYCoord_private > highestYCoord : highestYcoord_private ? highestYCoord
        ptsInside += ptsInside_private;
    }
}