【发布时间】：2017-04-11 20:29:23
【问题描述】：
我想编写一个轻量级的 PIC（Particle-in-cell）程序。我所说的“轻量级”是指它不需要扩大规模：假设所有数据都可以放入单个 GPU 设备的内存和主机系统的内存中。但是我希望它尽可能快。
问题是，PIC的典型结构是两个阶段的交互：场求解器和粒子推进器。工作流程是这样的： 初始化系统 -> 推动粒子 -> 求解场 -> 推动粒子 -> 求解场... -> 输出
下一个推动粒子或求解场必须等到前一个求解场或推动粒子完成。可能需要数百万次迭代才能获得最终输出。
作为测试，省略场求解器，粒子推进器可以写成：
__device__
void push(Particle &par) {
// some routines to move a particle. same excecutiong time for every particle.
}
并像这样使用 kernel_1 来执行它：
__global__
void kernel_1(int n, Particle* parlist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
if (i < n) {
push(parlist[i]);
}
}
在主循环中
for (int i=0;i<M;i++) {
kernel_1<<<(n+255)/256, 256>>>(n, parlist);
}
M 是所需的迭代次数。但是，性能非常缓慢：在我的八核 Intel E5-2640 v3 和 Nvidia Quadro m4000 系统上，CUDA 提供与使用 openmp 的纯 CPU 版本类似的性能。对于 10,000,000 且 M=1000 的粒子数，大约需要 10 秒。
但是，如果我将循环移入内核：
void kernel_2(int n, Particle* parlist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
if (i < n) {
for (int i=0;i<M;i++) {
push(parlist[i]);
}
}
}
和
kernel_2<<<(n+255)/256, 256>>>(n, parlist);
同样的M=1000，只需要100ms，就是10000%的加速。我已经验证了两种情况的结果相同且正确。可能内核运行M次调用成本太高了。
将循环移入内核所带来的性能提升是如此令人难以置信，但却是真实的。对于第一种情况，添加字段求解器很容易：只需编写一个新内核并在主循环中顺序执行两个内核。但是性能应该是中等的。
我发现很难将场解算器例程添加到第二种情况：在没有多次调用内核的情况下，块之间似乎没有同步机制，但是场解算器必须等到所有粒子都被推送，这必须分配到不同的块（因为粒子的数量非常多）。
那么是否可以在一个内核中实现两阶段迭代？性能提升太多不容忽视。
编辑： 我发现性能差异非常令人困惑：100ms 和 10s 的差异只是一行代码甚至是循环序列。我已将 push() 修改为更复杂一点（2d Boris pusher）：
class Particle
{
public:
float x, y; //m
float vx, vy; //m/s
float m; //kg
float q; //ee
};
__device__
void run(Particle& par, float B)
{
float t, s, vpx, vpy;
t = (par.q*ee*B/par.m)*dt/2;
s = 2*t/(1+t*t);
vpx = par.vx+t*par.vy;
vpy = par.vy-t*par.vx;
par.vx += s*vpy;
par.vy -= s*vpx;
par.x += par.vx*dt;
par.y += par.vy*dt;
}
我为 Particle 创建了 1 个 n 元素数组，为 B 创建了 1 个 n 元素浮点数组。它们是在主机和 cudaMemcpy 到设备上创建的。然后我检查了以下三个内核的性能：
__global__
void kernel_A(int n, int m, Particle* parlist, float* Blist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j;
if (i<n) {
for (j=0;j<m;j++) {
run(parlist[i], Blist[i]);
}
}
}
__global__
void kernel_B(int n, int m, Particle* parlist, float* Blist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j;
float B;
if (i<n) {
B = Blist[i];
for (j=0;j<m;j++) {
run(parlist[i], B);
}
}
}
__global__
void kernel_C(int n, int m, Particle* parlist, float* Blist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j;
float B;
if (i<n) {
B = Blist[i];
for (j=0;j<m;j++) {
run(parlist[i], B);
__syncthreads();
}
}
}
__global__
void kernel_D(int n, int m, Particle* parlist, float* Blist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j;
float B;
if (i<n) {
B = Blist[i];
}
for (j=0;j<m;j++) {
if (i<n) {
run(parlist[i], B);
}
}
}
__global__
void kernel_E(int n, int m, Particle* parlist, float* Blist)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j;
float B;
if (i<n) {
for (j=0;j<m;j++) {
run(parlist[i], Blist[i]);
__syncthreads();
}
}
}
而且运行时间完全不同。对于 n=10,000,000 和 m=1000：
- 内核_A：7.6s
- 内核_B：66ms
- 内核_C：9.9s
- 内核_D：10.0s
- 内核_E：10.0s
三个内核的结果都是一样的，都是正确的（检查CPU版本）。
我从官方 CUDA 编程指南了解到，分支很昂贵，因此 kernel_C 应该比 kernel_B 慢，尽管我怀疑差异是两个数量级。我不明白为什么 kernel_B 的性能比 kernel_A 好得多。 Kernel_B 不必访问 Blist 1000 次，而 kernel_A 则需要，但是它们都需要访问 parlist 1000 次，对吗？为什么访问 Blist 这么慢？
Kernel_A、kernel_D 和 kernel_E 有相似的性能，这让我更加困惑：所以与 kernel_B 相比，额外的时间花在访问 Blist 或同步上？
我想在我的 PIC 程序中实现 kernel_B 的性能。
【问题讨论】：
标签： c++ parallel-processing cuda