Cuda - Sieve of Eratosthenes

Here is an example of how to use cuda to find all prime numbers within a range.

Its uses the Sieve of Eratosthenes which creates an array initialized to 0 ( indicates prime) then you iterate over the array starting at zero and for each multiple of zero you mark it 1. You then move to the next element and repeat. When your down iterating all the primes are still marked as 0.

As an example cuda program I wrote the following kernel

__global__ static void Sieve(int * sieve,int sieve_size)
{ 
 int idx = blockIdx.x * blockDim.x + threadIdx.x; 
 if (idx > 1) {
  for(int i=idx+idx;i < sieve_size;i+=idx) 
   sieve[i] = 1;
 }
}

The index is figured out by the special keywords blockidx, blockDim and threadIdx. I then go over the multiples for that index and mark them as not prime.

The code to call this looks like this:

/************************************************************************/
/* Init CUDA                                                            */
/************************************************************************/
bool InitCUDA(void)
{
 int count = 0;
 int i = 0;

 cudaGetDeviceCount(&count);
 if(count == 0) {
  fprintf(stderr, "There is no device.\n");
  return false;
 }

 for(i = 0; i < count; i++) {
  cudaDeviceProp prop;
  if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
   if(prop.major >= 1) {
    break;
   }
  }
 }
 if(i == count) {
  fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
  return false;
 }
 cudaSetDevice(i);
 return true;
}
int main(int argc, char** argv)
{
  
 int  *device_sieve;
 int host_sieve[1000]; 
 double sieve_size = sizeof(host_sieve)/sizeof(int);

 if(!InitCUDA()) {
  return 0;
 } 
 cudaMalloc((void**) &device_sieve, sizeof(int) * sieve_size); 

 Sieve<<<1, sqrt(sieve_size), 0>>>(device_sieve, sieve_size);
 cudaThreadSynchronize();  
 cudaMemcpy(&host_sieve, device_sieve, sizeof(int) * sieve_size, cudaMemcpyDeviceToHost);
 
 cudaFree(device_sieve);
 
 for(int i=2;i < sieve_size;++i) 
  if (host_sieve[i] == 0)
   printf("%d is prime\n",i);
 

 return 0;
}

For simplicity I simply create a new thread for each multiple. This is not very efficient as I may be using a value that is already found to not be prime.