使用OpenMP并行优化循环时pcounter计数异常问题求助

问题描述 #

我尝试使用OpenMP对循环进行并行优化，以下是相关代码：

int temppcounter=0;
int pcounter=-1;
int nextX=0;
int nextY=0;
#pragma omp parallel for shared(matArr,visArr,plots) private(nextX,nextY) firstprivate(pcounter) lastprivate(pcounter)
for(int i = 0; i < rows; i++) {
   
   for(int j = 0; j < cols; j++) {
        if(matArr[i][j]!= ZERO) {
            continue;
        }
       
        for(int h = 0; h < NSIZE; h++) {
            nextX = i + dx[h];
            nextY = j + dy[h];

            if( nextX < 0 || nextY < 0 || nextX >= rows || nextY >= cols ) {
                continue;
            }
            
            if( !visArr[nextX][nextY] && matArr[nextX][nextY]== ONE){ //(int) image.at<uchar>(nextX, nextY) == ONE ) {
                pcounter++;
                visArr[nextX][nextY] = true;
                matArr[nextX][nextY]=pixelThreshold;
                plots[nextX][nextY]=1;
                plotx[pcounter]=nextX;
                ploty[pcounter]=nextY;
            }
        }
    }
}

#pragma omp parallel for   reduction(+:temppcounter) 
for(int i=0;i<rows;i++){
    #pragma omp parallel for
    for(int j=0;j<cols;j++){
        if(plots[i][j]==1){

            temppcounter++;
        }
    }

实际标记为1的plots[i][j]数量为382771，但pcounter的计数仅为147434；若将pcounter设为共享变量，虽能得到38万+的计数，但结果不正确。猜测问题出在pcounter上，但尚未找到根本原因，恳请提供相关提示或建议。

问题根本原因分析 #

• firstprivate(pcounter)的独立计数缺陷：firstprivate会给每个线程复制一份pcounter的初始值（-1），各线程独立递增自己的副本，lastprivate仅将最后完成循环的线程的副本值写回全局变量，导致最终pcounter只是单个线程的计数结果，远小于实际总数。
• 共享pcounter的竞态条件：设为共享变量后，多个线程同时执行pcounter++会引发竞态——线程间的读、递增、写回操作重叠，导致部分递增被覆盖；同时plotx[pcounter]和ploty[pcounter]会因pcounter值混乱出现数据覆盖，结果错误。
• 嵌套并行的冗余问题：第二个并行循环内嵌套#pragma omp parallel for会创建过量线程，增加调度开销，且内层并行与外层reduction存在逻辑冲突。

解决方案建议 #

1. 原子操作保证计数与索引安全
   将pcounter设为共享变量，用OpenMP原子操作确保递增和数组赋值的原子性：

int temppcounter=0;
int pcounter=-1;
int nextX=0;
int nextY=0;
#pragma omp parallel for shared(matArr,visArr,plots,pcounter,plotx,ploty) private(nextX,nextY)
for(int i = 0; i < rows; i++) {
    for(int j = 0; j < cols; j++) {
        if(matArr[i][j]!= ZERO) {
            continue;
        }
        for(int h = 0; h < NSIZE; h++) {
            nextX = i + dx[h];
            nextY = j + dy[h];
            if( nextX < 0 || nextY < 0 || nextX >= rows || nextY >= cols ) {
                continue;
            }
            if( !visArr[nextX][nextY] && matArr[nextX][nextY]== ONE){
                // 原子递增并获取当前索引，避免竞态
                #pragma omp atomic capture
                int idx = ++pcounter;
                visArr[nextX][nextY] = true;
                matArr[nextX][nextY]=pixelThreshold;
                plots[nextX][nextY]=1;
                plotx[idx]=nextX;
                ploty[idx]=nextY;
            }
        }
    }
}

2. 线程局部计数合并（高性能方案）
   减少原子操作开销，每个线程先维护局部计数器，最后合并到全局：

int temppcounter=0;
int pcounter=-1;
int nextX=0;
int nextY=0;
#pragma omp parallel shared(matArr,visArr,plots,pcounter,plotx,ploty) private(nextX,nextY)
{
    int local_pcounter = 0;
    // 可选：线程内临时存储坐标，用于后续拷贝到全局plotx/ploty
    std::vector<std::pair<int, int>> local_coords;

    #pragma omp for
    for(int i = 0; i < rows; i++) {
        for(int j = 0; j < cols; j++) {
            if(matArr[i][j]!= ZERO) {
                continue;
            }
            for(int h = 0; h < NSIZE; h++) {
                nextX = i + dx[h];
                nextY = j + dy[h];
                if( nextX < 0 || nextY < 0 || nextX >= rows || nextY >= cols ) {
                    continue;
                }
                if( !visArr[nextX][nextY] && matArr[nextX][nextY]== ONE){
                    // 原子标记已访问，避免重复处理
                    #pragma omp atomic write
                    visArr[nextX][nextY] = true;
                    local_pcounter++;
                    matArr[nextX][nextY]=pixelThreshold;
                    plots[nextX][nextY]=1;
                    local_coords.emplace_back(nextX, nextY);
                }
            }
        }
    }

    // 合并局部计数到全局，获取线程起始索引
    #pragma omp atomic capture
    int start_idx = pcounter += local_pcounter;
    // 将线程内坐标拷贝到全局数组
    for(int k = 0; k < local_coords.size(); k++) {
        plotx[start_idx + k] = local_coords[k].first;
        ploty[start_idx + k] = local_coords[k].second;
    }
}

此方案大幅减少原子操作次数，性能更优，适合大规模数据处理。

3. 修复第二个并行循环的嵌套问题
   移除内层并行，保留单层并行+reduction即可：

#pragma omp parallel for reduction(+:temppcounter) 
for(int i=0;i<rows;i++){
    for(int j=0;j<cols;j++){
        if(plots[i][j]==1){
            temppcounter++;
        }
    }
}

嵌套并行会引发不必要的线程调度开销，单层并行足以完成二维循环的并行化。

内容的提问来源于stack exchange，提问作者Lance