Find prime numbers from a given range of numbers parallel program

Question

I tried to implement the Sieve of Eratosthenes method to find prime numbers.

This is my code block

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<omp.h>


char isPrime[1000];     

// odd-only sieve
int eratosthenesOdd(int N, int useOpenMP)
{
        /* enable/disable OpenMP */
        omp_set_num_threads(useOpenMP ? omp_get_num_procs() : 1);

        /* instead of i*i <= N we write i <= sqr(N) to help OpenMP */
        const int NSqrt = (int)sqrt((double)N);
        int memorySize = (N-1)/2;


        int i, j;
        /* Let all numbers be prime initially */
        #pragma omp parallel for
        for (i = 0; i <= memorySize; i++){
                isPrime[i] = 1;
        }


        /* find all odd non-primes */
        #pragma omp parallel for schedule(dynamic)
        for (i = 3; i <= NSqrt; i += 2){
                if (isPrime[i/2]){
                        for (j = i*i; j <= N; j += 2*i){
                                isPrime[j/2] = 0;
                        }
                }
        }

        printf("2\n")   
        for(int k=3;k<=N;k+=2){
                if(isPrime[k/2]==1){
                        printf("%d ", k);
                }
        }

        /* sieve is complete, count primes */
        int total = N >= 2 ? 1 : 0;
        #pragma omp parallel for reduction(+:total)
        for (i = 1; i <= memorySize; i++){
                total += isPrime[i];
        }

        return total;
}

int main()
{
        double start, finish;
        start =  omp_get_wtime();
        int total = eratosthenesOdd(100, 1);
        finish = omp_get_wtime();
        printf("\n %d", total);
        printf("\n Elapsed time=%e seconds", finish-start);
        return 0;
}

I got reference from here. The code is running well and I can find how many prime numbers are in the given range.

I am skeptical about the elapsed time that I got after several trial with same number of term.

Lets suppose I want to see the prime numbers between 1 to 100. I also want to find out the elapsed time for various threads.

1st trial
N=100
Number of Thread  1, elapsed time = 5.008094e-04
Number of Thread  8, elapsed time = 4.649349e-04
Number of Thread 16, elapsed time = 4.652534e-04

2nd trial
N=100
Number of Thread  1, elapsed time = 4.668552e-04sec
Number of Thread  8, elapsed time = 5.837623e-04sec
Number of Thread 16, elapsed time = 5.835127e-04sec

 3rd trial
 N=100
Number of Thread  1, elapsed time = 4.530195e-04 sec
Number of Thread  8, elapsed time = 4.66317e-04sec
Number of Thread 16, elapsed time = 6.141420e-04 sec

I wonder is this program really implement parallel program? If so, how could I get this variation in times? Also, when task divide among threads the time should have decrease with compare to serial time. But here this was not happened.

Answer 1

I don't have the OMP software installed, so I had to remove the parallelization. However, the code works once the 'print out the prime numbers' loop is adjusted to deal with the fact that the algorithm knows that 2 is the only even prime number and that it stores the primality of an odd number X in isPrime[X / 2].

This leads to this code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <omp.h>

static char isPrime[1000];

// odd-only sieve
//static int eratosthenesOdd(int N, int useOpenMP)
static int eratosthenesOdd(int N)
{
    /* enable/disable OpenMP */
    //omp_set_num_threads(useOpenMP ? omp_get_num_procs() : 1);

    /* instead of i*i <= N we write i <= sqr(N) to help OpenMP */
    const int NSqrt = (int)sqrt((double)N);
    int memorySize = (N - 1) / 2;

    int i, j;
    /* Let all numbers be prime initially */
    //#pragma omp parallel for
    for (i = 0; i <= memorySize; i++)
    {
        isPrime[i] = 1;
    }


    /* find all odd non-primes */
    //#pragma omp parallel for schedule(dynamic)
    for (i = 3; i <= NSqrt; i += 2)
    {
        if (isPrime[i / 2])
        {
            for (j = i * i; j <= N; j += 2 * i)
            {
                isPrime[j / 2] = 0;
            }
        }
    }


    printf("2\n");
    for (int k = 3; k <= N; k += 2)
    {
        if (isPrime[k / 2] == 1)
        {
            printf("%d\n", k);
        }
    }

    /* sieve is complete, count primes */
    int total = N >= 2 ? 1 : 0;
    //#pragma omp parallel for reduction(+:total)
    for (i = 1; i <= memorySize; i++)
    {
        total += isPrime[i];
    }

    return total;
}

int main(void)
{
    //int total = eratosthenesOdd(100, 1);
    int total = eratosthenesOdd(100);
    printf("Number of primes: %d\n", total);
    return 0;
}

And this output:

2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97
Number of primes: 25

By inspection, this is correct for the primes up to 100.