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 的值可能非常大，无法分配l * l 的空间。

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

Answer 1

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

1. 您的内核代码在此处使用普通添加将无法正常工作：

   *count++;
   

   这是因为多个线程同时尝试执行此操作，而 CUDA 不会自动为您排序。出于解释的目的，我们将使用atomicAdd() 解决此问题，但也可以使用其他方法。

2. & 符号不属于这里：

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

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

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

3. 这种方法（有效地使用threadIdx.x 为一个维度和blockIdx.x 另一个维度创建一个线程的方形数组）只能在interactions 值导致l 值1024 或更少，因为 CUDA 线程块被限制为 1024 个线程，并且您在内核启动中使用 l 作为线程块的大小。要解决此问题，您需要了解如何创建任意维度的 CUDA 2D 网格，并适当地调整内核启动和内核内索引计算。现在，我们只需确保计算出的 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 维度的网格。您的测试实际上是在询问“该网格中的哪些点在一个圆圈内，中心位于网格的原点（角），半径为 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 的二维网格，应该至少可以工作到 10 亿个interactions。