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

Question

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

我正在编写的子程序依赖于在三个嵌套循环中使用的三维数组。对于 nsteps 的轨迹，这是一个典型的 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 运行到 C++ 扩展，该扩展将 python numpy 数组转换为 C 数组数据类型。 openacc 代码将被编译为源对象库。 C++ 扩展调用 openacc 代码。这是需要这种安排的遗留代码。在 Python 中定义系统、调用 C++ 扩展和访问 openacc *.so 文件的步骤工作正常。

问题在于，这个问题在本书或课程练习中并不明显，因此该问题的解决方案和/或 cmets 可能对其他人有所帮助。

我已经使用 Stack Overflow 上的一些示例和我找到的其他来源尝试了代码 sn-p 上的并行、内核和循环和数据指令。据我所知，Faber 书和其他来源中使用的示例并未解决此问题中提出的用例。充其量我可以让前两个循环并行化，但最里面的循环不并行化（循环未矢量化：数据依赖性）。当我正在寻找一般指导时，我不会发布失败的位，以便提供更多的教学讨论/提示。

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

1. 如何使用 pragma 指令处理输入位置数组的未知维度？
2. 如何管理 dist[] 数组的累积（由于粒子对的数量，它本身的长度在编译时是未知的）？
3. 一般建议使用哪些编译指示指令来解决此问题？
4. 如何处理“k-loop”对“j-loop”的依赖性？
5. 是否应该将问题展平以帮助定义要使用的指令？

谢谢，

SB

更新：

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

#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 计算该数组的索引，否则您将需要一个原子操作来保护对的增量，而这不会保证对于每一步，同一对点将位于距离数组中的同一位置。

有了这个问题，我认为steps循环上的acc data，j循环上的acc parallel loop和k循环上的简单acc loop（帮助编译器进行依赖分析）应该让您启动并运行。