为三角矩阵计算优化 CUDA 内核的执行

Question

我正在开发我的第一个 Cuda 应用程序，我有一个 "below-expected throughput" 的内核，这似乎是目前最大的瓶颈。

内核的任务是计算一个 N × N 大小的矩阵 (DD)，其中包含数据矩阵上所有元素之间的平方距离。数据矩阵 (Y) 的大小为 N x D（以支持多维数据）并存储为行优先。

来源：

__global__ void computeSquaredEuclideanDistance(const float * __restrict__ Y, float * __restrict__ DD, const int N, const int D) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;

    for (int i = index; i < N * N; i += stride) {
        const int m = i / N;
        const int n = i % N;
        float tmp = 0;
        for (int d = 0; d < D; ++d) {
            const float Ynd = Y[d + D * n];
            const float Ymd = Y[d + D * m];
            const float Ydiff = Ynd - Ymd;
            tmp += Ydiff * Ydiff;
        }
        DD[n + N * m] = tmp;
    }
}

这是用 size_t blockSize = 256 和 size_t numBlocks = (N*N + blockSize - 1)/blockSize 调用的。

如何优化这个内核？我最初的想法是，最耗时的部分是在不利用某种共享内存的情况下读取数据，但是谁能指导我如何处理这个问题？

nvvc 分析工具的备注：

延迟分析：
- 计算利用率约为 40%
- 内存（二级缓存）利用率约为 35%
入住不是问题
- 理论 64 的 57.59 的活动扭曲
- 理论 100 的 90% 的入住率

对于我的应用程序，典型值为：

5k < N < 30k
D 是 2 或 3

Answer 1

我通常会忽略这些类型的优化问题，因为在我看来，它们处于离题边缘。更糟糕的是，您没有提供 MCVE 因此任何试图回答的人都必须编写他们自己的所有支持代码来编译和基准测试您的内核。而这类工作确实需要基准测试和代码分析。但是因为你的问题基本上是一个线性代数问题（而且我喜欢线性代数），所以我回答了它而不是关闭投票它太宽泛......

随心所欲。代码中有几件事可以立即跳出，可以改进，并且可能会对运行时间产生 material 影响。

首先是内循环的行程计数是先验已知的。任何时候遇到这种情况，请让编译器知道。循环展开和代码重新排序是一种非常强大的编译器优化，NVIDIA 编译器非常擅长。如果将 D 移动到模板参数中，您可以这样做：

template<int D>
__device__ float esum(const float *x, const float *y)
{
    float val = 0.f;
#pragma unroll
    for(int i=0; i<D; i++) { 
        float diff = x[i] - y[i];
        val += diff * diff;
    }
    return val;
}

template<int D>
__global__ 
void vdistance0(const float * __restrict__ Y, float * __restrict__ DD, const int N)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    for (int i = index; i < N * N; i += stride) {
        const int m = i / N;
        const int n = i % N;
        DD[n + N * m] = esum<D>(Y + D * n, Y + D * m);
    }
}

template __global__ void vdistance0<2>(const float *, float *, const int);
template __global__ void vdistance0<3>(const float *, float *, const int);

编译器将内联 esum 并展开内部循环，然后它可以使用其重新排序试探法更好地交错加载和触发器以提高吞吐量。生成的代码也具有较低的寄存器占用空间。当我运行 N=10000 和 D=2 时，我得到了大约 35% 的加速（7.1 毫秒对带有 CUDA 9.1 的 GTX 970 上的 4.5 毫秒）。

但是还有比这更明显的优化。您正在执行的计算将产生一个对称的输出矩阵。您只需要执行 (N*N)/2 操作来计算完整矩阵，而不是您在代码中执行的 N*N [技术上 N(N/2 -1) 因为对角线条目为零，但让我们忘记对角线为了本次讨论的目的]。

所以采用不同的方法并使用一个块来计算上三角输出矩阵的每一行，然后你可以这样做：

struct udiag
{
    float *p;
    int m;

    __device__ __host__ udiag(float *_p, int _m) : p(_p), m(_m) {};
    __device__ __host__ float* get_row(int i) { return p + (i * (i + 1)) / 2; };
};


template<int D>
__global__ 
void vdistance2(const float * __restrict__ Y, float * __restrict__ DD, const int N)
{
     int rowid = blockIdx.x;
     int colid = threadIdx.x;
     udiag m(DD, N);

     for(; rowid < N; rowid += gridDim.x) {
         float* p = m.get_row(rowid);
         const float* y = Y + D * rowid;
         for(int i=colid; i < (N-rowid); i += blockDim.x) {
             p[i] = esum<D>(y, y + D * i);
         }
    }
}
template __global__ void vdistance2<2>(const float *, float *, const int);
template __global__ void vdistance2<3>(const float *, float *, const int);

这使用了一个小帮手class来封装上三角输出矩阵的寻址方案所需的三角数。这样做可以节省大量内存和内存带宽，并减少计算的总 FLOP 计数。如果您之后需要做其他事情，BLAS（和 CUBLAS）支持对上三角矩阵或下三角矩阵进行计算。使用它们。当我运行这样做时，我获得了大约 75% 的加速（7.1 毫秒对同一 GTX 970 上的 1.6 毫秒）。

重要免责声明：您在此处看到的所有代码都是在 45 分钟的午休期间编写的，并且非常进行了轻微测试。我绝对不声称此答案中的任何内容实际上是正确的。当我运行它获取分析数据时，我已经确认它编译并且不会产生运行时间错误。这就对了。 Cavaet Emptor 等等。

为三角矩阵计算优化 CUDA 内核的执行

Optimizing execution of a CUDA kernel for Triangular Matrix calculation

c++

cuda

distance-matrix