计算素数的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数组存储每个块的部分和。

---

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

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

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

3. 如何修改减量?

   __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版本中存在许多正确性问题，无论您正在实现的算法是否正确，都将阻止它正常工作:

1. __syncthreads()呼叫必须为无条件呼叫。编写分支分歧可能导致同一分支内的一些线程无法执行__syncthreads()调用的代码是不正确的。底层PTX是bar.sync, PTX指南是这样说的:

   barrier在每个warp的基础上执行，就好像a中的所有线程翘曲已激活。因此，如果经纱中的任何线程执行条指令，就好像warp中的所有线程都执行了条指令。经纱中的所有线都停止转动，直到到达屏障完成，并且屏障的到达计数增加经纱大小(不是经纱中活动线程的数量)。在有条件执行的代码，条形指令只能在以下情况下使用众所周知，所有线程对条件的求值都是相同的曲速不会发散)。因为barrier是在per-warp上执行的基础上，可选的支数必须是经纱尺寸的倍数。

你的代码在有条件地从全局内存加载一些值后无条件地将s_flags设置为1。这肯定不是法规的意图吧?
2. 代码在筛选代码和缩减代码之间缺乏同步屏障，这可能导致共享内存竞争和缩减的错误结果。
3. 如果你计划在Fermi类卡上运行这段代码，共享内存数组应该声明为volatile，以防止编译器优化可能破坏共享内存减少。

如果你修复了这些东西，代码可能会工作。性能是一个完全不同的问题。当然，在旧的硬件上，整数模操作非常非常慢，不推荐使用。我记得我读过一些资料，表明阿特金筛法是在gpu上快速生成质数的有效方法。