cuda - 用于计算素数的 CUDA 中的并行缩减

Question

我有一个代码来计算我使用 OpenMP 并行化的素数:

    #pragma omp parallel for private(i,j) reduction(+:pcount) schedule(dynamic)
    for (i = sqrt_limit+1; i < limit; i++)
    {
            check = 1;
            for (j = 2; j <= sqrt_limit; j++)
            {

                    if ( !(j&1) && (i&(j-1)) == 0 )
                    {
                            check = 0;
                            break;
                    }

                    if ( j&1 && i%j == 0 )
                    {
                            check = 0;
                            break;
                    }

            }
            if (check)
                pcount++;

    }

我正在尝试将它移植到 GPU，并且我想像上面的 OpenMP 示例一样减少计数。以下是我的代码，除了给出不正确的结果外，速度也较慢:

__global__ void sieve ( int *flags, int *o_flags, long int sqrootN, long int N) 
{
    long int gid = blockIdx.x*blockDim.x+threadIdx.x, tid = threadIdx.x, j;
    __shared__ int s_flags[NTHREADS];

    if (gid > sqrootN && gid < N)
            s_flags[tid] = flags[gid];
    else
            return;
    __syncthreads();

    s_flags[tid] = 1;

    for (j = 2; j <= sqrootN; j++)
    {
            if ( gid%j == 0 )
            {
                    s_flags[tid] = 0;
                    break;
            }
    }
    //reduce        
    for(unsigned int s=1; s < blockDim.x; s*=2)
    {
            if( tid % (2*s) == 0 )
            {
                    s_flags[tid] += s_flags[tid + s];
            }
            __syncthreads();
    }
    //write results of this block to the global memory
    if (tid == 0)
            o_flags[blockIdx.x] = s_flags[0];

}

首先，我如何让这个内核快速，我认为瓶颈是for循环，我不知道如何替换它。接下来，我的计数不正确。我确实更改了 '%' 运算符并注意到了一些好处。

在 flags数组，我已经标记了从 2 到 sqroot(N) 的素数，在这个内核中我正在计算从 sqroot(N) 到 N 的素数，但是我需要检查 {sqroot(N),N} 中的每个数字是否可整除通过 {2,sqroot(N)} 中的素数。 o_flags数组存储每个块的部分和。

编辑:按照建议，我修改了我的代码(我现在更好地理解了关于同步线程的评论)；我意识到我不需要 flags 数组，只需要全局索引就可以在我的情况下工作。在这一点上我担心的是代码的缓慢(不仅仅是正确性)，这可能归因于 for 循环。此外，在特定数据大小 (100000) 之后，内核会为后续数据大小生成不正确的结果。即使对于小于100000的数据大小，GPU缩减结果也是不正确的(NVidia论坛的一位成员指出，这可能是因为我的数据大小不是2的幂)。
所以还有三个(可能是相关的)问题——

我怎样才能使这个内核更快？在我必须遍历每个 tid 的情况下，使用共享内存是个好主意吗？

为什么它只对某些数据大小产生正确的结果？

我怎样才能修改减少？

__global__ void sieve ( int *o_flags, long int sqrootN, long int N )
{
unsigned int gid = blockIdx.x*blockDim.x+threadIdx.x, tid = threadIdx.x;
volatile __shared__ int s_flags[NTHREADS];

s_flags[tid] = 1;
for (unsigned int j=2; j<=sqrootN; j++)
{
       if ( gid % j == 0 )
            s_flags[tid] = 0;
}

__syncthreads();
//reduce
reduce(s_flags, tid, o_flags);
}

Answer 1

虽然我声称对筛选素数一无所知，但您的 GPU 版本中存在许多正确性问题，无论您实现的算法是否正确，这些问题都会阻止它正常工作:

__syncthreads()调用必须是无条件的。编写代码是不正确的，因为分支发散会导致同一经线中的某些线程无法执行 __syncthreads()称呼。底层PTX是bar.sync和 PTX 指南是这样说的:

Barriers are executed on a per-warp basis as if all the threads in a warp are active. Thus, if any thread in a warp executes a bar instruction, it is as if all the threads in the warp have executed the bar instruction. All threads in the warp are stalled until the barrier completes, and the arrival count for the barrier is incremented by the warp size (not the number of active threads in the warp). In conditionally executed code, a bar instruction should only be used if it is known that all threads evaluate the condition identically (the warp does not diverge). Since barriers are executed on a per-warp basis, the optional thread count must be a multiple of the warp size.

您的代码无条件设置 s_flags在有条件地从全局内存中加载一些值后变为一个。这肯定不是代码的意图吗？

代码在筛选代码和归约之间缺少同步屏障，这可能导致共享内存竞争和归约导致的错误结果。

如果您计划在费米类卡上运行此代码，则共享内存数组应声明为 volatile以防止编译器优化可能破坏共享内存减少。

如果您修复了这些问题，代码可能会起作用。性能是一个完全不同的问题。当然，在较旧的硬件上，整数模运算非常非常慢，不推荐使用。我记得读过一些 Material 暗示 Sieve of Atkin是一种在 GPU 上快速生成素数的有用方法。