如何在CUDA中有效地计算矩阵的所有列和行的总和？答案

【问题标题】：How to calculate the sum of all columns and rows of a matrix efficiently in CUDA?如何在CUDA中有效地计算矩阵的所有列和行的总和？
【发布时间】：2015-01-17 17:54:48
【问题描述】：

我想计算 CUDA 中矩阵的所有列的总和和所有行的总和。一种方法是使用 BLAS 中的SGEMV 子例程，将矩阵乘以 1 的向量。

但是，这会导致对矩阵进行两次扫描，假设它比 L1 缓存大得多：一次用于行，另一次用于列。此外，我计划进一步修改其他运算符的代码，所以这就是我编写自己的内核的原因。

到目前为止，我的方法是将矩阵分解为大小为32 x 32 的子矩阵。每个线程块将这样的子矩阵加载到共享内存中，计算子矩阵的行和列的总和，并将它们原子地添加到适当的输出（下面的row 和col）。这样，矩阵数据只需从 VRAM 中读取一次。

为简单起见，代码假设矩阵为n x n、n % 32 == 0，线程块为32 x 32

__global__ void sum_cols_and_rows(size_t n, const float* matrix, float* col, float* row)
{   
    __shared__ float sh[32][32];

    size_t x = blockDim.x * blockIdx.x + threadIdx.x;
    size_t y = blockDim.y * blockIdx.y + threadIdx.y;

    float sum = matrix[x + n * y];
    sh[threadIdx.x][threadIdx.y] = sum;

    for(unsigned w = 16; w >= 1; w /= 2)
        sum += __shfl_down(sum, w);
    const size_t laneID = threadIdx.x & 0x1f; // 32-1
    if(laneID == 0)
        atomicAdd(row + y, sum);
    __syncthreads();

    sum = sh[threadIdx.y][threadIdx.x]; // swapped indexes
    for(unsigned w = 16; w >= 1; w /= 2)
        sum += __shfl_down(sum, w);
    if(laneID == 0)
        atomicAdd(col + blockDim.x * blockIdx.x + threadIdx.y, sum);
}

// launch :
sum_cols_and_rows<<<dim3(n/32, n/32), dim3(32, 32), 32*32*sizeof(float)>>>(n, matrix, col, row);

然而，表现相当令人失望。我在 GTX 980 上看到大约 20% 的理论 224GB/s 内存带宽，即使在大型矩阵上，例如 16384x16384。

有没有办法让这种方法达到理论带宽限制？

【问题讨论】：

您可以尝试sh[32][33]; - 这可能有助于解决共享内存库冲突。除此之外，我不确定你是否受益于每个 NxN 块拥有 N^2 个线程，你可以尝试使用 N（也许使用更大的 N），根本不需要共享内存。
@zch sh[32][33] 给了我 50% 的加速，尽管我仍然处于理论极限的 30%。谢谢！我认为我需要共享内存来将数据从线程 (x,y) 发送到线程 (y, x) 在一个块中，并避免从 VRAM 重新读取该值。
并非如此。我的建议是每个块有 N 个线程，并且每个线程计算一个垂直总和。类似于：for(i 0..N-1) { float v = matrix[i][threadIdx]; vertical += v; horizontalShuffle(v); if(threadIdx==0) AtomicAdd(v); } AtomicAdd(vertical).
或者你可以保留共享内存并在没有任何洗牌的情况下执行for(i 0..N-1) { float v = matrix[i][threadIdx]; vertical += v; sh[i][threadIdx] = v; } AtomicAdd(vertical); sync(); for(i 0..N-1) { horizontal += sh[threadIdx][i]; }; AtomicAdd(horizontal);。对于足够大的矩阵，它可能是最有效的。
@zch 您的“无共享内存”建议让我获得了 71% 的带宽，这与我从 CUDA 样本中看到的“一维缩减”接近，所以这可能和它得到的一样好. 谢谢！（顺便说一句：起初我以为你在谈论列主要（FORTRAN）布局，它表现不佳——我猜是非合并访问）

标签： cuda gpgpu nvidia

【解决方案1】：

在您的解决方案中，矩阵的每个 NxN 块都由单独的 NxN 线程块处理。实际上，每个单独的线程所做的工作很少，因此开销支配了实际计算。您可以通过让线程块处理多个矩阵块来改进它。

但是有一个更简单的解决方案，每个矩阵块只使用 N 个线程，其中一个线程对整列求和。

实现将与此类似：

__global__ void sum_cols_and_rows(size_t n, const float* matrix, float* col, float* row)
{   
    size_t laneID = threadIdx.x & 31;

    size_t x = blockDim.x * blockIdx.x + threadIdx.x;
    size_t y = N_ITERATIONS * blockIdx.y;

    size_t idx = y * n + x;

    float vertical = 0;

    for(int i = 0; i < N_ITERATIONS; i++) {
        float v = matrix[idx];
        vertical += v;
        for(unsigned w = 16; w >= 1; w /= 2)
            v += __shfl_down(v, w);
        if(laneID == 0)
            atomicAdd(&row[y], v);
        y++;
        idx += n;
    }

    atomicAdd(&col[x], vertical);
}

这里的可调参数是每个线程组的扭曲数和每个矩阵块中的行数 (N_ITERATIONS)。较大的值可能会减少开销，但会以并行性为代价。

另一个可以试验的想法是vectorized loading - 其中之一：

float2 v2 = reinterpret_cast<float2*>(matrix)[idx];
float v = v2.x + v2.y;

float4 v4 = reinterpret_cast<float4*>(matrix)[idx];
float v = (v4.x + v4.y) + (v4.z + v4.w);

【讨论】：