CUDA

Question

我正在编写埃拉托色尼筛法的实现（https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes) on GPU. But no sth like this - http://developer-resource.blogspot.com/2008/07/cuda-sieve-of-eratosthenes.html

方法：

创建默认值为 0/1（0 - 质数，1 - 否）的 n 元素数组并将其传递给 GPU（我知道它可以直接在内核中完成，但目前这不是问题）。
块中的每个线程检查单个数字的倍数。每个块检查总的 sqrt(n) 可能性。每个块==不同的间隔。
将倍数标记为 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

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

Answer 2

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

您基本上是在访问超出范围的项目。在你的内核中考虑这个序列：
```
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 代码时遇到问题时以及在寻求他人帮助之前，您应该始终做这些事情。
只有当输入 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 的实际平方根。这将导致您计算出比必要大小小的总线程数组。（如果你现在不相信我也没关系。运行我将在下面 post 的代码并输入 1025 之类的大小，然后查看线程数 * 块数是否足够大覆盖一个 1025 的数组。)
正如您所说：

Each block checks in total sqrt(n) possibilities.

希望这也指出非完美平方 n 的危险，但我们现在必须问“如果 n 大于最大线程块大小 (1024) 的平方怎么办？答案是该代码在许多情况下无法正常工作——您选择的 10240000 输入虽然是一个完美的正方形，但超过了 1024^2 (1048576)，因此无法正常工作。您的算法（我声称是 not a Sieve of 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;
}

CUDA - 埃拉托色尼筛分法

CUDA - Sieve of Eratosthenes division into parts

c

parallel-processing

primes

sieve-of-eratosthenes