Calculate PI using OpenMP task directive

Question

I need to parallelize code that calculates number π using Leibniz formula for π with OpenMP task directive.

Leibniz formula

So, I got a sequential code:

double sequential_execution(long long n)
{
    long long i;
    double factor;
    double sum = 0.0;
    double startTime = omp_get_wtime();

    for (i = 0; i < n; i++) {
        factor = (i % 2 == 0) ? 1.0 : -1.0;
        sum += factor / (2 * i + 1);
    }
    double endTime = omp_get_wtime();
    printf("Sequential execution took %f seconds\n", endTime - startTime);
    sum = 4.0 * sum;
    return sum;
}

My first idea was to capture content of for loop as a single task with n = 100000000:

double parallel_execution(long long n)
{
    long long i=0;
    double factor;
    double sum = 0.0;
    long long index; 
    long squareRootN = ceil(sqrt(n));

    double startTime = omp_get_wtime();
#pragma omp parallel default(none) private(i,factor) shared(n,sum) 
{
    #pragma omp single
    {
        for ( i = 0; i < n; i++) {
            #pragma omp task
            {
                factor = (i % 2 == 0) ? 1.0 : -1.0;
                #pragma omp atomic
                sum += factor / (2 * i + 1);
            }
        }
    }
}
    double endTime = omp_get_wtime();
    printf("Parallel execution took %f seconds\n", endTime - startTime);
    sum = 4.0 * sum;
    return sum;
}

But sequential execution was way way faster.(Seq. time: 0.3 s, Par. time: 87 s)

Second idea was to increase granularity of one task and decrease number of tasks in a way where one for loop that goes from 0 do n-1 was split into two nested loops where each one goes from 0 to sqrt(n)-1. Now, each task has for loop that goes from 0 to sqrt(n)-1, and sqrt(n) tasks are generated, again for n = 100000000.

double parallel_execution(long long n)
{
    long long i=0;
    double factor;
    double sum = 0.0;
    long long index; 
    long squareRootN = ceil(sqrt(n));

    double startTime = omp_get_wtime();
#pragma omp parallel default(none) shared(sum,n,squareRootN) private(i,factor,index)
{
    #pragma omp single
    {
        for (i=0;i<squareRootN;i++)
        #pragma omp task
        {
            for (long j=0;j<squareRootN;j++)
            {
                index = i*squareRootN + j;
                if (index > n) break;
                factor = (index % 2 == 0)?1.0 : -1.0; 
                #pragma omp atomic
                sum += factor / (2*index + 1);
            }
        }
    }
}
    double endTime = omp_get_wtime();
    printf("Parallel execution took %f seconds\n", endTime - startTime);
    sum = 4.0 * sum;
    return sum;
}

Now, I got better time, but again it was way slower than sequential execution(Seq : 0.3s, Par : 11s).

At this point, I'm starting to think that it's not possible to get speed-up using task directive, but again, is there something that I did wrong or is there some way to restructure the problem to get better performances ? Thanks

Edit: Best function so far:

double parallel_execution(long long n)
{
    double factor;
    int totalThreads = 0;
    long squareRootN = ceil(sqrt(n));
    double master_sum = 0;
    double *sum;
    double startTime = omp_get_wtime();
#pragma omp parallel default(none) shared(sum,n,squareRootN,totalThreads) private(factor)
{
    #pragma omp single
    {
        totalThreads = omp_get_num_threads();
        sum = (double*)calloc(totalThreads,sizeof(double));
        for (long long i=0;i<squareRootN;i++)
        #pragma omp task
        {
            for (long long j=0;j<squareRootN;j++)
            {
                long long index = i*squareRootN + j;
                if (index > n) break;
                factor = (index % 2 == 0)?1.0 : -1.0; 
                sum[omp_get_thread_num()] += factor / (2*index + 1);
            }
        }
    }
}
    for (int i=0;i<totalThreads;i++) master_sum += sum[i];
    double endTime = omp_get_wtime();
    printf("Parallel execution took %f seconds\n", endTime - startTime);
    master_sum*=4;
    return master_sum;
}

Input size: n = 1000000000 Seq. time: 3.19 s Par. time: 4 s

Answer 1

You are paying the overheads of the atomic operation and task creation and management. You can get a better speedup with a simpler parallel for with reduction, namely:

#pragma omp parallel default(none) shared(n) reduction( + : sum ) 
for ( i = 0; i < n; i++) {
     double factor = (i % 2 == 0) ? 1.0 : -1.0;
     sum += factor / (2 * i + 1);
}

We can slightly improve the sequential code by separating beforehand the odds from the evens:

#pragma omp parallel default(none) shared(n, sum) nowait
{
     #pragma omp for reduction( + : sum ) 
     for (int i = 0; i < n; i+=2 ) {
        sum += 1.0 / (2 * i + 1);
    }
    #pragma omp for reduction( + : sum ) 
    for (int i = 1; i < n; i += 2) {
        sum += -1.0 / (2 * i + 1);
    }
}

You can improve it further by having a single loop for performing the even and odds computation for each iteration of that loop.

You do not need to make the 'i' from the loop private, it will be implicitly private in OpenMP.

If you really have to use tasks, you can try to minimize the synchronization overhead by replicating the variable sum among threads, and reduce it manually at the end of the parallel region, (I am assuming n >= 2 and n being even just for simplicity sake):

double sum[total_threads];

#pragma omp parallel default(none) shared(n, sum)
{
    int threadID = omp_get_thread_num();
    sum[threadID] = 0.0;
    #pragma omp single
    {
        for ( i = 0; i < n; i+=2) {
            #pragma omp task
            {
                sum[threadID] += 1.0 / (2 * i + 1);
                sum[threadID] += -1.0 / (2 * (i + 1) + 1);
            }
        }
    }
  }

double master_sum = 0.0;
for(int i = 0; i < total_threads; i++)
    master_sum += sum[i];

If you are using a C compiler that supports OpenMP 4.5 you can use a more sophisticated constructor, namely taskloop Construct, and combined it with the reduction of the variable sum.