使用OpenACC collapse子句折叠过多循环层级时出现错误结果

问题背景 #

在基于MPI与OpenACC的大规模分析工作中，遇到由OpenACC #pragma指令中collapse子句引发的异常：使用collapse(2)可得到预期结果，使用collapse(3)时结果完全错误，尽管两者循环执行次数一致，且错误具有确定性。以下针对三个核心问题逐一解答：

核心问题 #

1. 为何collapse(3)会产生错误结果，即便循环求和次数相同？
2. 为何collapse(2)可避免该问题？
3. 使用collapse(2)时，是否存在可加速最内层j循环且保证结果正确的OpenACC编译指令？

问题解答 #

1. collapse(3)导致错误的原因 #

当使用collapse(3)时，编译器会将x、y、j三层循环合并为一个并行循环空间，分配给多个GPU线程同时执行。此时，同一个sum_XY[x][y]会被多个线程并发写入：不同j值对应的线程会同时执行sum_XY[x][y] += a[xx][j]操作，而sum_XY并未被标记为归约（reduction）变量，也没有同步机制。这种无保护的并发写入会引发数据竞争，导致sum_XY的累加结果被随机覆盖，最终得到错误值。

虽然nbrSums的结果正确，是因为编译器自动识别到它的累加操作并生成了隐式归约（从编译输出的Generating implicit reduction(+:nbrSums)可验证），但sum_XY的依赖关系更复杂，编译器无法自动为其添加归约，因此出现错误。

2. collapse(2)避免问题的原因 #

使用collapse(2)时，仅合并x和y两层循环，最内层j循环被编译器标记为串行执行（从编译输出的46, #pragma acc loop seq可验证）。这意味着每个(x,y)对应的j循环由单个线程完整执行，sum_XY[x][y]的累加是单线程内的顺序操作，不存在多线程并发写入的数据竞争，因此结果完全正确。

编译器通过依赖分析检测到：sum_XY[x][y]的更新依赖于自身的前一次值（循环携带依赖），如果并行化j循环会导致竞争，因此自动将j设为串行。

3. 加速最内层j循环的正确指令 #

可以通过显式添加归约（reduction）子句来并行化j循环，同时保证累加结果正确。有两种实现方式：

方式一：在并行指令中添加全局归约

修改#pragma指令，将sum_XY声明为归约变量：

#pragma acc parallel loop gang device_type(acc_device_nvidia) vector_length(256) collapse(2) private(x,y,j,xx,yy) reduction(+:sum_XY[:blockSize][:maxSize])
for ( x=0; x<blockSize; x++ ) {
    for ( y=0; y<maxSize; y++ ) {
        for ( j=0; j<ncols; j++ ) {
            xx = x + blockStart;
            yy = xx - y;
            if ( yy >= 0 ) {
                sum_XY[x][y] += a[xx][j];
                nbrSums += 1;
            }
        }
    }
}

方式二：在j循环前添加局部归约

在j循环上方单独添加归约指令：

#pragma acc parallel loop gang device_type(acc_device_nvidia) vector_length(256) collapse(2) private(x,y,j,xx,yy)
for ( x=0; x<blockSize; x++ ) {
    for ( y=0; y<maxSize; y++ ) {
        #pragma acc loop reduction(+:sum_XY[x][y])
        for ( j=0; j<ncols; j++ ) {
            xx = x + blockStart;
            yy = xx - y;
            if ( yy >= 0 ) {
                sum_XY[x][y] += a[xx][j];
                nbrSums += 1;
            }
        }
    }
}

两种方式都会让编译器为每个sum_XY[x][y]创建局部副本，并行完成j循环的累加后再合并结果，既避免了数据竞争，又实现了j循环的并行加速。相比原子操作（atomic），归约的性能更优，适合大规模循环场景。

复现代码与编译结果 #

复现代码 #

#include <mpi.h>

void* allocMatrix (int nRow, int nCol) {
    void* restrict m = malloc (sizeof(int[nRow][nCol]));
    return(m);
}

int main(int argc, char *argv[]) {
    int nrows = 100;
    int ncols = 15;
    int blockSize = 10;
    int blockOffset = blockSize/2;
    int blockStart;
    int maxSize = 5;

    int (*a)[ncols] = (int (*)[ncols])allocMatrix(nrows, ncols);
    int (*sum_XY)[maxSize] = (int (*)[maxSize])allocMatrix(blockSize, maxSize);

    int nbrBlocks = (nrows-blockSize)/(blockSize/2) + 1;
    int nbrSums = 0;

    int blockNbr, x, y, j, xx, yy;
    int sumAll = 0;

    memset( a, 1, sizeof(int)*nrows*ncols );

    for ( blockNbr=0; blockNbr<nbrBlocks; blockNbr++ ) {
        blockStart = (blockOffset) * blockNbr;

        for ( x=0; x<blockSize; x++ ) {
            for ( y=0; y<maxSize; y++ ) {
                sum_XY[x][y] = 0;
            }
        }
     
#pragma acc parallel loop gang device_type(acc_device_nvidia) vector_length(256) collapse(2) private(x,y,j,xx,yy)
//#pragma acc parallel loop gang device_type(acc_device_nvidia) vector_length(256) collapse(3) private(x,y,j,xx,yy)
        for ( x=0; x<blockSize; x++ ) {
            for ( y=0; y<maxSize; y++ ) {
                for ( j=0; j<ncols; j++ ) {
                    xx = x + blockStart;
                    yy = xx - y;
                    if ( yy >= 0 ) {
                        sum_XY[x][y] += a[xx][j];
                        nbrSums += 1;
                    }
                }
            }
        }

        for ( x=0; x<blockSize; x++ ) {
            for ( y=0; y<maxSize; y++ ) {
                sumAll += sum_XY[x][y];
            }
        }

    }

    printf( "sumAll : %d\n", sumAll );
    printf( "nbrSums: %d\n", nbrSums );
}

使用collapse(2)的编译与运行结果 #

$ mpic++ -mcmodel=medium -acc -ta=tesla:managed -Minfo=accel MRE_collapse.c -o MRE_collapse
nvc++-Warning-CUDA_HOME has been deprecated. Please, use NVHPC_CUDA_HOME instead.
main:
     44, Generating NVIDIA GPU code
         44, #pragma acc loop gang, vector(256) collapse(2) /* blockIdx.x threadIdx.x */
             Generating implicit reduction(+:nbrSums)
         45,   /* blockIdx.x threadIdx.x collapsed */
         46, #pragma acc loop seq
     44, Generating implicit copy(sum_XY[:10][:5],nbrSums) [if not already present]
         Generating implicit copyin(a[blockStart:10][:15]) [if not already present]
     46, Complex loop carried dependence of sum_XY->,a-> prevents parallelization
         Loop carried dependence of sum_XY-> prevents parallelization
         Loop carried backward dependence of sum_XY-> prevents vectorization
$ mpirun  -np 1 MRE_collapse
sumAll : 1263225620
nbrSums: 14100
$

使用collapse(3)的编译与运行结果 #

$ mpic++ -mcmodel=medium -acc -ta=tesla:managed -Minfo=accel MRE_collapse.c -o MRE_collapse
nvc++-Warning-CUDA_HOME has been deprecated. Please, use NVHPC_CUDA_HOME instead.
main:
     44, Generating NVIDIA GPU code
         44, #pragma acc loop gang, vector(256) collapse(3) /* blockIdx.x threadIdx.x */
             Generating implicit reduction(+:nbrSums)
         45,   /* blockIdx.x threadIdx.x collapsed */
         46,   /* blockIdx.x threadIdx.x collapsed */
     44, Generating implicit copy(sum_XY[:10][:5],nbrSums) [if not already present]
         Generating implicit copyin(a[blockStart:10][:15]) [if not already present]
$ mpirun  -np 1 MRE_collapse
sumAll : -1347440724
nbrSums: 14100
$

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