Parallel Reduction in CUDA for calculating primes

Question

I have a code to calculate primes which I have parallelized using 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++;

    }

I am trying to port it to GPU, and I would want to reduce the count as I did for the OpenMP example above. Following is my code, which apart from giving incorrect results is also slower:

__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];

}

First of all, how do I make this kernel fast, I think the bottleneck is the for loop, and I am not sure how to replace it. And next, my counts are not correct. I did change the '%' operator and noticed some benefit.

In the flags array, I have marked the primes from 2 to sqroot(N), in this kernel I am calculating primes from sqroot(N) to N, but I would need to check whether each number in {sqroot(N),N} is divisible by primes in {2,sqroot(N)}. The o_flags array stores the partial sums for each block.

---

EDIT: Following the suggestion, I modified my code (I understand about the comment on syncthreads now better); I realized that I do not need the flags array and just the global indexes work in my case. What concerns me at this point is the slowness of the code (more than correctness) that could be attributed to the for loop. Also, after a certain data size (100000), the kernel was producing incorrect results for subsequent data sizes. Even for data sizes less than 100000, the GPU reduction results are incorrect (a member in the NVidia forum pointed out that that may be because my data size is not of a power of 2). So there are still three (may be related) questions -

1. How could I make this kernel faster? Is it a good idea to use shared memory in my case where I have to loop over each tid?

2. Why does it produce correct results only for certain data sizes?

3. How could I modify the reduction?

   __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

While I profess to know nothing about sieving for primes, there are a host of correctness problems in your GPU version which will stop it from working correctly irrespective of whether the algorithm you are implementing is correct or not:

1. __syncthreads() calls must be unconditional. It is incorrect to write code where branch divergence could leave some threads within the same warp unable to execute a __syncthreads() call. The underlying PTX is bar.sync and the PTX guide says this:

   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.

2. Your code unconditionally sets s_flags to one after conditionally loading some values from global memory. Surely that cannot be the intent of the code?
3. The code lacks a synchronization barrier between the sieving code and the reduction, this can lead to a shared memory race and incorrect results from the reduction.
4. If you are planning on running this code on a Fermi class card, the shared memory array should be declared volatile to prevent compiler optimization from potentially breaking the shared memory reduction.

If you fix those things, the code might work. Performance is a completely different issue. Certainly on older hardware, the integer modulo operation was very, very slow and not recommended. I can recall reading some material suggesting that Sieve of Atkin was a useful approach to fast prime generation on GPUs.