CUDA - Eratosthenes 筛分法

Question

我正在 GPU 上编写 Sieve of Eratosthenes (https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes) 的实现。但是没有这样的 - http://developer-resource.blogspot.com/2008/07/cuda-sieve-of-eratosthenes.html

方法：

1. 创建具有默认值 0/1（0 - 素数，1 - 否）的 n 元素数组并将其传递到 GPU（我知道它可以直接在内核中完成，但目前没有问题）。李>
2. 块中的每个线程检查单个数字的倍数。每个块检查总 sqrt(n) 可能性。每个块 == 不同的间隔。
3. 将倍数标记为 1 并将数据传回主机。

代码：

#include <stdio.h>
#include <stdlib.h>
#define THREADS 1024

__global__ void kernel(int *global, int threads) {
    extern __shared__ int cache[];

    int tid = threadIdx.x + 1;
    int offset = blockIdx.x * blockDim.x;
    int number = offset + tid;

    cache[tid - 1] = global[number];
    __syncthreads();

    int start = offset + 1;
    int end = offset + threads;

    for (int i = start; i <= end; i++) {
        if ((i != tid) && (tid != 1) && (i % tid == 0)) {
            cache[i - offset - 1] = 1;
        }
    }
    __syncthreads();
    global[number] = cache[tid - 1];
}


int main(int argc, char *argv[]) {
    int *array, *dev_array;
    int n = atol(argv[1]);
    int n_sqrt = floor(sqrt((double)n));

    size_t array_size = n * sizeof(int);
    array = (int*) malloc(n * sizeof(int));
    array[0] = 1;
    array[1] = 1;
    for (int i = 2; i < n; i++) {
        array[i] = 0;
    }

    cudaMalloc((void**)&dev_array, array_size);
    cudaMemcpy(dev_array, array, array_size, cudaMemcpyHostToDevice);

    int threads = min(n_sqrt, THREADS);
    int blocks = n / threads;
    int shared = threads * sizeof(int);
    kernel<<<blocks, threads, shared>>>(dev_array, threads);
    cudaMemcpy(array, dev_array, array_size, cudaMemcpyDeviceToHost);

    int count = 0;
    for (int i = 0; i < n; i++) {
        if (array[i] == 0) {
            count++;
        }
    }
    printf("Count: %d\n", count);
    return 0;
}

跑步： ./筛 10240000

当 n = 16, 64, 1024, 102400 时它可以正常工作......但是对于 n = 10240000 我得到不正确的结果。问题出在哪里？

Answer 1

在我看来，这段代码有很多问题。

1. 您基本上是在访问超出范围的项目。在你的内核中考虑这个序列：

   int tid = threadIdx.x + 1;
   int offset = blockIdx.x * blockDim.x;
   int number = offset + tid;
   
   cache[tid - 1] = global[number];
   

   您（在某些情况下——见下文）已经启动了一个与您的global 数组大小完全相同的线程数组。那么当最高编号的线程运行上述代码时会发生什么？ number = threadIdx.x+1+blockIdx.x*blockDim.x。这个number 索引将超出数组末尾。这对于n 的许多可能值都是正确的。如果您使用proper cuda error checking 或使用cuda-memcheck 运行您的代码，这个问题对您来说很明显。当您遇到 CUDA 代码问题时以及在寻求他人帮助之前，您应该始终做这些事情。

2. 只有当输入 n 是一个完美的正方形时，代码才有可能正常工作。其原因包含在这些代码行中（以及内核中的依赖项）：

   int n = atol(argv[1]);
   int n_sqrt = floor(sqrt((double)n));
   ...
   int threads = min(n_sqrt, THREADS);
   int blocks = n / threads;
   

   （注意这里正确的函数是atoi 而不是atol，但我离题了...）除非n 是一个完美的正方形，否则生成的n_sqrt 会比 n 的实际平方根。这将导致您计算出比所需大小小 的总线程数组。 （此时如果你不相信我也没关系。运行我将在下面发布的代码并输入一个像 1025 这样的大小，然后查看线程数 * 块的大小是否足以覆盖 1025 的数组。）

3. 正如你所说：

   每个块检查总 sqrt(n) 个可能性。

   希望这也指出了非完美平方n 的危险，但我们现在必须问“如果n 大于最大线程块大小的平方（1024）怎么办？答案是代码在许多情况下将无法正常工作 - 您选择的输入 10240000 虽然是一个完美的正方形，但超过了 1024^2 (1048576)，因此它不起作用。您的算法（我声称它不是Eratosthenes) 要求每个块都能够检查sqrt(n) 可能性，正如您在问题中所述。当由于每个块的线程限制而不再可能时，您的算法就会开始中断。

这是一个尝试修复上述问题 #1 的代码，并且至少解释了与 #2 和 #3 相关的故障：

#include <stdio.h>
#include <stdlib.h>
#define THREADS 1024
#define MAX 10240000

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)


__global__ void kernel(int *global, int threads) {
    extern __shared__ int cache[];

    int tid = threadIdx.x + 1;
    int offset = blockIdx.x * blockDim.x;
    int number = offset + tid;

    if ((blockIdx.x != (gridDim.x-1)) || (threadIdx.x != (blockDim.x-1))){
      cache[tid - 1] = global[number];
      __syncthreads();

      int start = offset + 1;
      int end = offset + threads;

      for (int i = start; i <= end; i++) {
        if ((i != tid) && (tid != 1) && (i % tid == 0)) {
            cache[i - offset - 1] = 1;
        }
      }
      __syncthreads();
      global[number] = cache[tid - 1];}
}


int cpu_sieve(int n){
    int limit = floor(sqrt(n));
    int *test_arr = (int *)malloc(n*sizeof(int));
    if (test_arr == NULL) return -1;
    memset(test_arr, 0, n*sizeof(int));
    for (int i = 2; i < limit; i++)
      if (!test_arr[i]){
        int j = i*i;
        while (j <= n){
          test_arr[j] = 1;
          j += i;}}
    int count = 0;
    for (int i = 2; i < n; i++)
      if (!test_arr[i]) count++;
    return count;
}

int main(int argc, char *argv[]) {
    int *array, *dev_array;
    if (argc != 2) {printf("must supply n as command line parameter\n"); return 1;}
    int n = atoi(argv[1]);
    if ((n < 1) || (n > MAX)) {printf("n out of range %d\n", n); return 1;}
    int n_sqrt = floor(sqrt((double)n));

    size_t array_size = n * sizeof(int);
    array = (int*) malloc(n * sizeof(int));
    array[0] = 1;
    array[1] = 1;
    for (int i = 2; i < n; i++) {
        array[i] = 0;
    }

    cudaMalloc((void**)&dev_array, array_size);
    cudaMemcpy(dev_array, array, array_size, cudaMemcpyHostToDevice);

    int threads = min(n_sqrt, THREADS);
    int blocks = n / threads;
    int shared = threads * sizeof(int);
    printf("threads = %d, blocks = %d\n", threads, blocks);
    kernel<<<blocks, threads, shared>>>(dev_array, threads);
    cudaMemcpy(array, dev_array, array_size, cudaMemcpyDeviceToHost);
    cudaCheckErrors("some error");
    int count = 0;
    for (int i = 0; i < n; i++) {
        if (array[i] == 0) {
            count++;
        }
    }
    printf("Count: %d\n", count);
    printf("CPU Sieve: %d\n", cpu_sieve(n));
    return 0;
}

Answer 2

我认为有几个问题，但这里有一个指向实际问题的指针：Eratosthenes 的筛子迭代地删除已经遇到的素数的倍数，并且您希望将工作负载分成线程块，其中每个线程块都在一块共享内存（在您的示例中为缓存）上运行。然而，线程块通常独立于所有其他线程块，并且不能轻松地相互通信。一个例子来说明这个问题：索引为0的线程块中索引为0的线程删除了2的倍数。索引> 0的线程块无法知道这一点。