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

Question

我正在开发我的第一个Cuda应用程序，目前遇到的最大问题是内核的吞吐量低于预期。

这个内核的任务是计算一个N乘N大小的矩阵(DD)，其中包含数据矩阵上所有元素之间的平方距离。数据矩阵(Y)的大小为N乘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％
  • 内存（L2缓存）利用率约为35％
• 占用率不是问题
  • 活动线程束占理论最大值的57.59％
  • 占用率占理论最大值的90％

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

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

- casparjespersen

---

D和N的范围是什么？有哪些典型值？对于优化来说，这些事情很重要。 - Robert Crovella

通常情况下，5k < N < 30k，而 D 可以是2或3。 - casparjespersen

你可能想在“cuda”标签下搜索“euclid”。你会发现许多相关的问题和答案，包括这个。 - Robert Crovella

你可以用 float2 或 float3 替换内部循环 (for (int d = 0; d < D; ++d)) ... - Lanting

这儿难道不是有很多对称性可以利用的吗？ - talonmies

我在相关的C++问题中添加了一个讨论解决此问题。 - Parker

Answer 1

我通常会忽略这类优化问题，因为在我看来它们接近于离题。更糟糕的是，您没有提供MCVE，因此任何试图回答您的人都必须编写自己的支持代码来编译和基准测试您的内核。而这种工作确实需要基准测试和代码分析。但是由于您的问题基本上是线性代数问题（而且我喜欢线性代数），所以我还是回答了它，而不是将其关闭。

说完这些，有几件事情立即突显出代码可以改进并可能对运行时间产生实质影响。

第一件事是，内部循环的旅程计数是预先知道的。每当您遇到这种情况时，请让编译器知道。循环展开和代码重新排序是非常强大的编译器优化，而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%（在CUDA 9.1上的GTX 970上，7.1ms与4.5ms相比）。

但是，还有一个更加显而易见的优化。你正在执行的计算将产生对称输出矩阵。你只需要做(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);

这里使用了一个小帮助类来封装为上三角形矩阵寻址方案所需的三角形数。这样做可以节省大量内存和内存带宽，同时减少了计算的总FLOP数。如果你需要在之后进行其他操作，则BLAS（和CUBLAS）支持对上或下三角形矩阵进行计算。建议使用它们。当我运行此代码时，我的速度提升了约75％（与同一GTX 970上的7.1ms相比，现在只需要1.6ms）。

免责声明：所有您在此处看到的代码都是在45分钟的午餐时间内编写的，并且经过了非常轻微的测试。本回答中的任何内容都不作保证是正确的。我确认它可以编译并在运行时不会产生运行时错误以获取性能分析数据。仅此而已。

- talonmies