CUDA 中的二元矩阵缩减

Question

我必须遍历一个假想矩阵的所有单元格 m * n 和 add + 1 满足特定条件的所有单元格。

我天真的解决方案如下：

#include <stdio.h>

__global__ void calculate_pi(int center, int *count) {
    int x = threadIdx.x;
    int y = blockIdx.x;

    if (x*x + y*y <= center*center) {
        *count++;
    }
}

int main() {
    int interactions;
    printf("Enter the number of interactions: ");
    scanf("%d", &interactions);

    int l = sqrt(interactions);

    int h_count = 0;
    int *d_count;

    cudaMalloc(&d_count, sizeof(int));
    cudaMemcpy(&d_count, &h_count, sizeof(int), cudaMemcpyHostToDevice);

    calculate_pi<<<l,l>>>(l/2, d_count);

    cudaMemcpy(&h_count, d_count, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(d_count);

    printf("Sum: %d\n", h_count);

    return 0;
}

在我的用例中，interactions 的值可能非常大，因此无法分配 space 的 l * l。

有人可以帮助我吗？欢迎提出任何建议。

Answer 1

您的代码至少有 2 个问题：

您的内核代码将无法在此处进行普通添加：
```
*count++;
```
这是因为多个线程试图同时执行此操作，而 CUDA 不会自动为您解决这个问题。出于解释的目的，我们将使用 atomicAdd() 解决此问题，尽管其他方法也是可能的。

& 号不属于这里：

cudaMemcpy(&d_count, &h_count, sizeof(int), cudaMemcpyHostToDevice);
           ^

我认为这只是一个错字，因为您在随后的 cudaMemcpy 操作中正确地做到了：

cudaMemcpy(&h_count, d_count, sizeof(int), cudaMemcpyDeviceToHost);

这种方法（使用 threadIdx.x 一维和 blockIdx.x 另一维有效地创建方形线程阵列）最多只能使用 interactions导致 l 值为 1024 或更小的值，因为 CUDA 线程块限制为 1024 个线程，并且您在内核启动中使用 l 作为线程块的大小。要解决此问题，您需要学习如何创建任意维度的 CUDA 二维网格，并适当调整内核启动和内核内索引计算。现在我们将确保计算出的 l 值在您的代码设计范围内。

下面是解决上述问题的示例：

$ cat t1590.cu
#include <stdio.h>

__global__ void calculate_pi(int center, int *count) {
    int x = threadIdx.x;
    int y = blockIdx.x;

    if (x*x + y*y <= center*center) {
        atomicAdd(count, 1);
    }
}

int main() {
    int interactions;
    printf("Enter the number of interactions: ");
    scanf("%d", &interactions);

    int l = sqrt(interactions);
    if ((l > 1024) || (l < 1)) {printf("Error: interactions out of range\n"); return 0;}
    int h_count = 0;
    int *d_count;

    cudaMalloc(&d_count, sizeof(int));
    cudaMemcpy(d_count, &h_count, sizeof(int), cudaMemcpyHostToDevice);

    calculate_pi<<<l,l>>>(l/2, d_count);

    cudaMemcpy(&h_count, d_count, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(d_count);
    cudaError_t err = cudaGetLastError();
    if (err == cudaSuccess){
      printf("Sum: %d\n", h_count);
      printf("fraction satisfying test:  %f\n", h_count/(float)interactions);
      }
    else
      printf("CUDA error: %s\n", cudaGetErrorString(err));
    return 0;
}
$ nvcc -o t1590 t1590.cu
$ ./t1590
Enter the number of interactions: 1048576
Sum: 206381
fraction satisfying test:  0.196820
$

我们看到代码表明计算出的分数约为 0.2。这看起来是正确的吗？我声称根据您的测试，它看起来确实是正确的。您正在有效地创建一个表示 lxl 维度的网格。您的测试实际上是在询问 "which points in that grid are within a circle, with the center at the origin (corner) of the grid, and radius l/2 ?"

从图片上看，它看起来像这样：

假设红色阴影区域略小于总面积的 0.25 是合理的，因此 0.2 是对该区域的合理估计。

作为奖励，这里是减少上面第 3 项中列出的限制的代码版本：

#include <stdio.h>

__global__ void calculate_pi(int center, int *count) {
    int x = threadIdx.x+blockDim.x*blockIdx.x;
    int y = threadIdx.y+blockDim.y*blockIdx.y;

    if (x*x + y*y <= center*center) {
        atomicAdd(count, 1);
    }
}

int main() {
    int interactions;
    printf("Enter the number of interactions: ");
    scanf("%d", &interactions);

    int l = sqrt(interactions);
    int h_count = 0;
    int *d_count;
    const int bs = 32;
    dim3 threads(bs, bs);
    dim3 blocks((l+threads.x-1)/threads.x, (l+threads.y-1)/threads.y);

    cudaMalloc(&d_count, sizeof(int));
    cudaMemcpy(d_count, &h_count, sizeof(int), cudaMemcpyHostToDevice);

    calculate_pi<<<blocks,threads>>>(l/2, d_count);

    cudaMemcpy(&h_count, d_count, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(d_count);
    cudaError_t err = cudaGetLastError();
    if (err == cudaSuccess){
      printf("Sum: %d\n", h_count);
      printf("fraction satisfying test:  %f\n", h_count/(float)interactions);
      }
    else
      printf("CUDA error: %s\n", cudaGetErrorString(err));
    return 0;
}

这是基于 l 推出的 2D 网格，应该至少达到 10 亿 interactions 。

CUDA 中的二元矩阵缩减

Binary Matrix Reduction in CUDA

parallel-processing

cuda

nvidia

pi