带有动态数组的 openacc 嵌套循环

Question

我正在尝试应用 openacc 来开发多核和 gpu 加速二进制文件。我已经阅读了 Farber 的书并通过 NVIDIA 提供的一些在线课程成功地运行测试了程序。然后，我尝试并行化我们的遗留代码。

我正在编写的子程序依赖于在三个嵌套循环中使用的三维数组。这是 n 步轨迹的典型 N(N+1)/2 成对距离问题。对于我们的科学问题，nsteps 通常是 1E5 到 1E7，nparticles 是 1E4 到 5E5。

for(i=0 ; i < nsteps ; i++){
    pair = 0 ;
    for(j=0 ; j < nparticles-1 ; j++){
        x1 = position[i][j][0] ;
        y1 = position[i][j][1] ;
        z1 = position[i][j][2] ;
        for(k=j+1 ; k < nparticles ; k++){
            x2 = position[i][k][0] ;
            y2 = position[i][k][1] ;
            z2 = position[i][k][2] ;
            dx2 = (x1 - x2) * (x1 - x2) ;
            dy2 = (y1 - y2) * (y1 - y2) ;
            dz2 = (z1 - z2) * (z1 - z2) ;
            sdist = sqrt(dx2 + dy2 + dz2) ;
            dist[pair] += sdist ;
            pair++ ;
        }
    }
}

我在编译时无法控制输入位置数组（nsteps、nparticles），因为代码是运行到 python 到将 python numpy 转换为 C++ 扩展数组到 C 数组数据类型。 openacc 代码将被编译为源对象库。 C++ 扩展调用 openacc 代码。这是一个遗留代码，需要这个 ar运行gement。 Python 中定义系统的步骤，调用 C++ 扩展并访问 openacc *.so 文件工作正常。

问题是这个问题在书中或课程中的练习中并不明显，因此解决方案 and/or 对该问题的评论可能会对其他人有所帮助。

我已经使用 Stack Overflow 和我找到的其他来源的一些示例在显示的代码片段上尝试了并行、内核和循环和数据指令。据我所知，Faber 书中使用的示例和其他来源并未解决该问题中提出的用例。充其量我可以让前两个循环并行化，但最内层的循环不会并行化（循环未矢量化：数据依赖性）。由于我正在寻找一般指导，因此我不会发布失败的内容，以便提供更具教学意义的讨论/技巧。

好的，现在回答我的问题。

如何使用 pragma 指令处理输入位置数组的未知维度？
如何管理 dist[] 数组的累积（由于粒子对的数量，它本身的长度在编译时是未知的）？
对于这个问题，通常推荐哪些 pragma 指令？
如何处理 "k-loop" 对 "j-loop" 的依赖？
是否应该将问题展平以帮助定义要使用的指令？

谢谢，

SB

更新：

为了根据@jefflarkin 的建议提供结果，我更改了代码以说明 dist 数组的索引并添加了建议的编译指示。该子例程可以并行编译运行s。我现在将开始分析以了解如何优化例程以最大限度地利用资源。这是工作代码的副本：

#pragma acc data copyin(position[nsteps][nparticles][3]) copy(dist[npairs])
for(i=0 ; i < nsteps ; i++){
    #pragma acc parallel loop
    for(j=0 ; j < nparticles-1 ; j++){
        x1 = position[i][j][0] ;
        y1 = position[i][j][1] ;
        z1 = position[i][j][2] ;
        #pragma acc loop
        for(k=j+1 ; k < nparticles ; k++){
            x2 = position[i][k][0] ;
            y2 = position[i][k][1] ;
            z2 = position[i][k][2] ;
            dx2 = (x1 - x2) * (x1 - x2) ;
            dy2 = (y1 - y2) * (y1 - y2) ;
            dz2 = (z1 - z2) * (z1 - z2) ;
            sdist = sqrt(dx2 + dy2 + dz2) ;
            local_count = ((j*nparticles)-((j*(j+1))/2))+k-(j+1) ;
            dist[local_count] += sdist ;
        }
    }
}

我得到了这个编译器结果（pgc++ (16.10): CFLAGS= -fPIC -c -fast -acc -Minfo=accel -ta=multicore -O3 -shared）

 38, Generating Multicore code
 39, #pragma acc loop gang
 45, Loop is parallelizable

对于 GPU (CFLAGS= -v -fPIC -c -fast -acc -Minfo=accel -ta=tesla:cuda8,fastmath -O3 -shared)

 35, Generating copyin(coor[:nframes][:nparticles][:3])
     Generating copy(dist[:npairs])
 38, Accelerator kernel generated
     Generating Tesla code
     39, #pragma acc loop gang /* blockIdx.x */
     45, #pragma acc loop vector(128) /* threadIdx.x */
 45, Loop is parallelizable

其中第 35 行 i-loop (nsteps)，第 38 行是 j-loop，第 45 行是 k-loop。

Answer 1

你说数组的大小是未知的，但是你可以根据nparticles为你的数据子句计算它们。位置数组的大小为 [nparticles][nparticles][3]。 dist 数组有点棘手，因为 k 对于每个 j 都会缩小，但它是 (nparticles - 1) + (nparticles - 2) + ... + 1 的总和，我认为是 (nparticles * (nparticles - 1) )/2。尽管直到运行时才能确定粒子数，但编译器必须了解有关数组的一些信息，否则 [] 将无法工作。

我建议将索引更改为 dist 数组，以便您可以根据 j 和 k 计算该数组的索引，否则您将需要一个原子操作来保护 pair 的增量，这将'保证对于每一步，同一对点将位于距离数组中的相同位置。

固定后，我想在步骤循环上使用 acc data，在 j 循环上使用 acc parallel loop，在 k 循环上使用简单的 acc loop（以帮助编译器与依赖分析）应该让你起来运行.

带有动态数组的 openacc 嵌套循环

openacc nested loops with dynamic arrays

c++

multidimensional-array

openacc