优化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%
  • 内存(二级缓存)利用率约为35%
• 占用不是问题
  • 主动变形为理论值64的57.59
  • 占据率为理论值的90%

对于我的应用,典型值是:

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

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.1毫秒,而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);

这使用一个小辅助类来封装上三角输出矩阵的寻址方案所需的三角形数.这样做可以节省大量的内存和内存带宽,并减少计算的总FLOP数量.如果之后需要做其他事情,BLAS(和CUBLAS)支持上三角矩阵或下三角矩阵的计算.使用它们.当我运行这个时,我获得了大约75%的加速(在相同的GTX 970上为7.1毫秒与1.6毫秒).

巨大的免责声明:您在这里看到的所有代码都是在45分钟的午休时间内编写的,并且经过了非常轻微的测试.我绝对没有声称这个答案中的任何内容实际上是正确的.我已经确认它在运行它以获取分析数据时编译并且不会产生运行时错误.这就对了.Cavaet Emptor等等.