在 cuda C++ 中加速成对力计算答案

【问题标题】：accelerate pairwise force calculation in cuda C++在 cuda C++ 中加速成对力计算
【发布时间】：2022-08-14 20:16:41
【问题描述】：

最近我写了一个分子动力学代码，使用 CUDA 并行计算计算离子-电子力。内核列表如下：

__global__ void (*x,*y,*z,N){
    int i = (blockIdx.x * blockDim.x) + threadIdx.x;
    while(i<N) {
        double dx;
        double dy;
        double dz;
        double dr;
        double Fx;
        double Fy;
        double Fz;

        for (int j = 0; j < N; j++){
            dx=x[i]-x[j];
            dy=y[i]-y[j];
            dz=z[i]-z[j];
            dr=sqrt(dx*dx+dy*dy+dz*dz) 
            dr=dr*dr*dr
            Fx+=k*q*q*dx/dr
            Fy+=k*q*q*dy/dr
            Fz+=k*q*q*dz/dr        //force=kq^2r/r^3 written in Cartesian coordinate
            }
        //rest of the code manipulate force is irrelevant to my question and I want to keep my code short
        i += blockDim.x * gridDim.x;
    }
}

x,y,z 是粒子的位置，dx,dy,dz 是 xyz 距离，for 循环中的 Fx,Fy,Fz 是施加在第 i 个粒子上的力的总和，更具体地说，您需要计算 x[i ]-x[j] 并遍历所有 js 以找到合力，让内核并行执行所有 i。

我发现这很慢，因为我知道 GPU 正在从全局内存中读取数组。当我将 x[i] 更改为一个数字时，它会变得快 10 倍，因为它正在从寄存器（L1 缓存）中读取。我的数组太大（超过 20000 个元素，双浮点数）无法放入寄存器。但是使用其他内存还能快一点吗？我知道有常量内存和共享内存，但我不知道如何实现。我认为 x[i] 位于全球内存中，导致它很慢，并且所有线程都试图同时读取 x[i]。有什么办法可以提高速度？

我会担心正确性而不是性能。您的代码无法计算可重复的正确结果。它甚至不将任何内容存储到全局内存中，这意味着如果您在优化的情况下编译代码，它应该编译为执行时间为零的空内核
我在发布此代码时确实更改了我的代码，原始代码很长，需要对这些变量和算法进行更仔细的处理，抱歉，只显示我的部分代码，它过于简单，Fx Fy 和 Fz 绝对需要存储在某个地方，我想念那部分。我的问题是由于循环，每个线程都在读取相同的 x[i] N 次并读取 x[j] N^2 次。有什么方法可以减少读取相同变量的次数或加快读取变量的速度
分块平铺方法可以正常工作。基本上将 i 和 j 视为矩阵中的行和列。使用与优化矩阵-矩阵乘法相同的分块评估方案。如果我有时间，我可能稍后会写一个正确的答案
发布无法编译的、损坏的代码并询问优化策略有什么意义？细节很重要
请注意，在这样的 N 体代码中，由于1 / dr 因子对于“长”距离非常小（它以O(1 / (n^3)) 的速率减小），因此假设某些力可以忽略是很常见的。因此，您通常可以丢弃大部分计算而不存储它。四叉树和 AMR 方法有助于做到这一点（尽管它并不简单）。此外，存储结果通常不是一个好主意：您需要动态计算它以获得快速代码。现在内存比计算单元慢得多，而且它不会很快变得更好（恰恰相反）。

标签： performance memory cuda pairwise

【解决方案1】：

这是一个使用共享内存来优化访问模式的基本版本。

#define KERNEL_BLOCKSIZE 256

__global__ void __launch_bounds__(KERNEL_BLOCKSIZE)
kernel(const double* x, const double* y, const double* z, int N,
       double k, double q, double* fake_out)
{
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    /*
     * threads beyond the bound still participate in value fetching, so we cannot
     * return early
     */
    const bool active = i < N;
    double xi, yi, zi;
    if(active)
        xi = x[i], yi = y[i], zi = z[i];
    const double kqq = k * q * q;
    double Fx = 0., Fy = 0., Fz = 0.;
    __shared__ double xt[KERNEL_BLOCKSIZE];
    __shared__ double yt[KERNEL_BLOCKSIZE];
    __shared__ double zt[KERNEL_BLOCKSIZE];
    for(int j = 0; j < N; j += blockDim.x) {
        __syncthreads();
        const int thread_j = j + threadIdx.x;
        if(thread_j < N) {
            xt[threadIdx.x] = x[thread_j];
            yt[threadIdx.x] = y[thread_j];
            zt[threadIdx.x] = z[thread_j];
        }
        __syncthreads();
        for(int l = 0, M = min(KERNEL_BLOCKSIZE, N - j); l < M; ++l) {
            const double dx = xi - xt[l], dy = yi - yt[l], dz = zi - zt[l];
            // 1 / sqrt(dx*dx + dy+dy + dz*dz)
            const double rnorm = rnorm3d(dx, dy, dz);
            const double dr = rnorm * rnorm * rnorm;
            const double scale = kqq * dr;
            Fx += scale * dx;
            Fy += scale * dy;
            Fz += scale * dz;
        }
    }
    if(active)
        fake_out[i] = norm3d(Fx, Fy, Fz);
}

这没什么花哨的，也没有解决 O(N²) 运行时的固有问题。我做了以下更改

摆脱 while 循环。循环计数器被声明为int i。所有 CUDA 设备中的最大网格尺寸为 2^31-1。这意味着我们总是可以启动整个网格，每个线程只有一个循环。

考虑到二次运行时，无论如何，我们没有机会运行如此巨大的网格。但是如果我们确实有一个更大的，只需启动多个在子集上运行的内核

使用共享内存来缓冲块。我选择了 256 作为固定块大小。这往往运作良好。 512可能是另一个值得尝试的尺寸
整个dr 计算可以折叠成一个预定义的数学函数
为了得到至少可以编译成合理代码的东西，我添加了一个输出

双缓冲

我们可以通过使用双缓冲来减少所需的__syncthreads() 的数量。但是，这会使共享内存使用量翻倍。只有 64 kiB 共享内存的平台将受到有限的占用。它需要进行基准测试以查看哪个版本效果更好。

__global__ void __launch_bounds__(KERNEL_BLOCKSIZE)
kernel_dbuf(const double* x, const double* y, const double* z, int N,
            double k, double q, double* fake_out)
{
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    const bool active = i < N;
    double xi, yi, zi;
    if(active)
        xi = x[i], yi = y[i], zi = z[i];
    const double kqq = k * q * q;
    double Fx = 0., Fy = 0., Fz = 0.;
    __shared__ double xt[2][KERNEL_BLOCKSIZE];
    __shared__ double yt[2][KERNEL_BLOCKSIZE];
    __shared__ double zt[2][KERNEL_BLOCKSIZE];
    int dbuf = 0;
    for(int j = 0; j < N; dbuf ^= 1, j += blockDim.x) {
        const int thread_j = j + threadIdx.x;
        if(thread_j < N) {
            xt[dbuf][threadIdx.x] = x[thread_j];
            yt[dbuf][threadIdx.x] = y[thread_j];
            zt[dbuf][threadIdx.x] = z[thread_j];
        }
        __syncthreads();
        for(int l = 0, M = min(KERNEL_BLOCKSIZE, N - j); l < M; ++l) {
            const double dx = xi - xt[dbuf][l];
            const double dy = yi - yt[dbuf][l];
            const double dz = zi - zt[dbuf][l];
            // 1 / sqrt(dx*dx + dy+dy + dz*dz)
            const double rnorm = rnorm3d(dx, dy, dz);
            const double dr = rnorm * rnorm * rnorm;
            const double scale = kqq * dr;
            Fx += scale * dx;
            Fy += scale * dy;
            Fz += scale * dz;
        }
    }
    if(active)
        fake_out[i] = norm3d(Fx, Fy, Fz);
}

像这样启动内核：

__host__ void
launch(const double* x, const double* y, const double* z, int N,
      double k, double q, double* fake_out, cudaStream_t stream)
{
    const int numBlocks = (N + KERNEL_BLOCKSIZE - 1) / KERNEL_BLOCKSIZE;
    kernel<<<numBlocks, KERNEL_BLOCKSIZE, 0, stream>>>(x, y, z, N, k, q, fake_out);
}

其他想法

人们已经评论了算法固有的低效率
我想k 和q 是单独的变量是有充分理由的，而且您不只是将预先计算的k * q * q 传递给内核
在我看来，在 GPU 上进行计算时，使用双打应该始终是最后的手段。降低精度的可能途径，至少对于部分算法：
- 将dr 计算替换为不易发生溢出的计算。像这样：
```
float scale = 1.f / max(max(abs(dx), abs(dy)), abs(dz));
float rnorm = rnorm3df(dx * scale, dy * scale, dz * scale) * scale;
float dr = rnorm * rnorm * rnorm;
```
- 对Fx、Fy、Fz 使用 Kahan 求和
- 仅对 Fx、Fy、Fz 使用双精度，但不适用于 x、y、z 位置或其他计算

【讨论】：