CUDA-Kernel 应该是动态崩溃取决于块大小

Question

我想做一个稀疏矩阵，密集向量乘法。让我们假设压缩矩阵中条目的唯一存储格式是压缩行存储 CRS。

我的内核如下所示：

__global__ void
krnlSpMVmul1(
        float *data_mat,
        int num_nonzeroes,
        unsigned int *row_ptr,
        float *data_vec,
        float *data_result)
{
    extern __shared__ float local_result[];
    local_result[threadIdx.x] = 0;

    float vector_elem = data_vec[blockIdx.x];

    unsigned int start_index = row_ptr[blockIdx.x];
    unsigned int end_index = row_ptr[blockIdx.x + 1];

    for (int index = (start_index + threadIdx.x); (index < end_index) && (index < num_nonzeroes); index += blockDim.x)
        local_result[threadIdx.x] += (data_mat[index] * vector_elem);

    __syncthreads();

   // Reduction

   // Writing accumulated sum into result vector
}

正如你所看到的，内核应该尽可能地幼稚，它甚至会做一些错误的事情（例如，vector_elem 并不总是正确的值）。我知道这些事情。

现在我的问题： 假设我正在使用 32 或 64 个线程的块大小。只要我的矩阵中的一行有超过 16 个非零值（例如 17 个），就只有前 16 个乘法完成并保存到共享内存中。我知道local_result[16] 的值是第 17 次乘法的结果只是零。使用 16 或 128 个线程的块大小可以解决已解释的问题。

由于我是 CUDA 的新手，我可能忽略了最简单的事情，但我无法编造更多情况来查看。

非常感谢您的帮助！

---

编辑对 talonmies 评论：

我在计算后直接打印了local_result[16] 中的值。它是0。然而，这里是缺少的代码：

还原部分：

int k = blockDim.x / 2;
while (k != 0)
{
    if (threadIdx.x < k)
        local_result[threadIdx.x] += local_result[threadIdx.x + k];
    else
        return;

    __syncthreads();

    k /= 2;
}

以及我如何将结果写回全局内存：

data_result[blockIdx.x] = local_result[0];

这就是我的全部。

现在我正在测试一个矩阵，该矩阵由单行和 17 个元素组成，这些元素都是非零的。缓冲区在伪代码中如下所示：

float data_mat[17] = { val0, .., val16 }
unsigned int row_ptr[2] = { 0, 17 }
float data_vec[17] = { val0 } // all values are the same
float data_result[1] = { 0 }

这是我的包装函数的摘录：

float *dev_data_mat;
unsigned int *dev_row_ptr;
float *dev_data_vec;
float *dev_data_result;

// Allocate memory on the device
HANDLE_ERROR(cudaMalloc((void**) &dev_data_mat, num_nonzeroes * sizeof(float)));
HANDLE_ERROR(cudaMalloc((void**) &dev_row_ptr, num_row_ptr * sizeof(unsigned int)));
HANDLE_ERROR(cudaMalloc((void**) &dev_data_vec, dim_x * sizeof(float)));
HANDLE_ERROR(cudaMalloc((void**) &dev_data_result, dim_y * sizeof(float)));

// Copy each buffer into the allocated memory
HANDLE_ERROR(cudaMemcpy(
        dev_data_mat,
        data_mat,
        num_nonzeroes * sizeof(float),
        cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(
        dev_row_ptr,
        row_ptr,
        num_row_ptr * sizeof(unsigned int),
        cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(
        dev_data_vec,
        data_vec,
        dim_x * sizeof(float),
        cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(
        dev_data_result,
        data_result,
        dim_y * sizeof(float),
        cudaMemcpyHostToDevice));

// Calc grid dimension and block dimension
dim3 grid_dim(dim_y);
dim3 block_dim(BLOCK_SIZE);

// Start kernel
krnlSpMVmul1<<<grid_dim, block_dim, BLOCK_SIZE>>>(
        dev_data_mat,
        num_nonzeroes,
        dev_row_ptr,
        dev_data_vec,
        dev_data_result);

我希望这很简单，但如果有任何兴趣，我会解释一下。

还有一件事：我刚刚意识到使用 128 的 BLOCK_SIZE 和 33 个非零值也会导致内核失败。同样只是最后一个值没有被计算。

Answer 1

您动态分配的共享内存大小不正确。现在你正在这样做：

krnlSpMVmul1<<<grid_dim, block_dim, BLOCK_SIZE>>>(.....)

共享内存大小应以字节为单位。每个块使用 64 个线程，这意味着您将为 16 个浮点大小的字分配足够的共享内存，并解释了为什么每行的魔术 17 个条目会导致失败 - 您有共享缓冲区溢出，这将触发保护错误GPU 并中止内核。

你应该这样做：

krnlSpMVmul1<<<grid_dim, block_dim, BLOCK_SIZE * sizeof(float)>>>(.....)

这将为您提供正确的动态共享内存大小，并且应该可以消除问题。