在CUDA中使用小M的两个MxN矩阵执行逐矢量点积的最快方法是什么？

Question

我有两个矩阵,每个矩阵都MxN在哪里,M = 16并且N更大(n = 262144比如说).我的目标是产生一个长度向量N,其中每个元素对应nth于每个矩阵中向量的点积.

我尝试了以下方法,其中cIdx对应于每个矩阵中列向量的列索引.毫不奇怪,NVIDIA Visual Profiler告诉我这种方法主要是内存带宽限制.

    public static void MatrixDotProduct(
        float* matrix1,
        float* matrix2,
        float* dotProduct,
        int2 matrixDimensions)
    {
        int i = blockIdx.x * blockDim.x + threadIdx.x;
        int stride = gridDim.x * blockDim.x;
        float sum;

        for (int cIdx = i; cIdx < matrixDimensions.y; cIdx += stride)
        {
            int ci = cIdx * matrixDimensions.x;
            sum = 0f;

            for (int j = 0; j < matrixDimensions.x; j++)
            {
                sum += matrix1[ci + j] * matrix2[ci + j];
            }

            dotProduct[cIdx] = sum;
        }
    }

我还发现了一个矢量矢量点积的版本,我试图并行化.不幸的是,这比上面的实现慢了20%(也许是因为我的M = 16？).有没有更好的方法来解决我在这里缺少的问题？

Answer 1

由于缺乏数据重用和较小的矢量大小(小于扭曲),这不是一个容易优化的情况.代码应该是内存绑定的,缓存性能至关重要.

要尝试的一件事是通过使用载荷的矢量类型来减少体积并提高内存事务的效率.我希望有类似的东西:

__device__ float vdot(float2 v1, float2 v2) {
    return (v1.x * v2.x) + (v1.y * v2.y);
}

__device__ float vdot(float4 v1, float4 v2) {
    return (v1.x * v2.x) + (v1.y * v2.y) + (v1.z * v2.z) + (v1.w * v2.w);
}

template<typename VT, int NT>
__device__ float vector_dotprod(const VT* v1, const VT* v2) {
    float sum = 0.f;
#pragma unroll
    for (int j = 0; j < NT; j++) {
        sum += vdot(v1[j], v2[j]);
    }
    return sum;
}

template<typename VT, int Nrows>
__global__
void MatrixDotProductPlus(float* matrix1, float* matrix2, float* dotProduct, int2 matrixDimensions)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    int stride2 = stride * Nrows;

    VT* m1 = reinterpret_cast<VT*>(matrix1) + i * Nrows;
    VT* m2 = reinterpret_cast<VT*>(matrix2) + i * Nrows;
    for (; i < matrixDimensions.y; i += stride, m1 += stride2, m2 += stride2) {
        dotProduct[i] = vector_dotprod<VT,Nrows>(m1, m2);
    }
}

[警告:只有非常轻微的测试 - 使用风险自负]

大约是float2案例代码的两倍,float4对于长度为16的向量,Maxwell或Pascal架构的大约快四倍.这假设您知道编译时向量的长度,它们是圆的倍数2或4,虽然我怀疑循环展开与矢量类型本身一样重要.