Какой самый быстрый способ выполнения векторно-точечных произведений для двух матриц MxN с маленьким M в CUDA?

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. Предполагается, что вы знаете длину векторов во время компиляции, и они кратны 2 или 4, хотя я сомневаюсь, что развертывание цикла столь же критично, как и сам тип вектора.