Why OpenMP is slower than a sequential program for a simple reduction?

Question

I am trying to find sum of elements in an array as I show below. However, the OpenMP implementation surprisingly is slower than the sequential implementation. I tried with both heap allocated and stack allocated arrays and got similar results. Any help is greatly appreciated.

#include <iostream>
#include <omp.h>
int main() {
  int N = 10000;
  int * ary = new int[N];
  for (int i = 0; i < N; i++) { input_file >> ary[i]; }
  int sum = 0;
  clock_t begin = clock();
  for (int i = 0; i < N; i++) { sum += ary[i]; }
  clock_t end = clock();
  cout << sum;
  double elapsed_time = double(end - begin) / CLOCKS_PER_SEC;
  sum = 0;
  begin = clock();
  #pragma omp parallel
  {
    int thread_id = omp_get_thread_num();
    int total_threads = omp_get_num_threads();
    int elem_per_thread = N / total_threads;
    int base = thread_id * elem_per_thread;
    int internal_sum = 0;
    for (int i = base; i < (base + elem_per_thread); i++) {
      internal_sum += ary[i];
    }
    #pragma omp critical
    {
      sum += internal_sum;
    }
  }
  end = clock();
  cout << sum;
  elapsed_time = double(end - begin) / CLOCKS_PER_SEC;    
}

The sequential program takes 5e-06 (s) to finish and the parallel one takes 0.001733 (s). I am compiling on Ubuntu 16.04 using g++ -std=c++11 main.cpp -fopenmp -O3 && ./a.out

Answer 1

The sequential program optimizes down to doing nothing. This is because the only side effect is the value of sum, and the value of sum is not observable in your program.

With the OpenMP, the complexity of threading things off prevents the compiler from realizing you aren't doing anything.

A simple way that could avoid this is add return sum; now it shows up as an exit code, which is observable, and hence the calculation cannot be optimized away.

Now, the compiler is still free to never allocate ary, because it can prove that ary[i]==i for all i, and replace reading ary[i] with just i, then calculate at compile time that the sum of i from 1 to 10000 is 50005000, eliminate your entire loop and make it a sum=50005000 and still take zero time.

Answer 2

Remark beforehand:
Handling the way the loop is divided "by hand" I believe is counterproductive (unless you want to understand how OpenMP works). That's why I first propose you use the more standard approach with reduction operation. You can always check that it gives the same result in terms of performance.
Another remark is that using throughout your code omp_ functions will not able you to compile it without -openmp option.

Benching

So I benched with the following code:

headers

#include <iostream>
#include <fstream>
#include <omp.h>
#include <cmath>
#include <chrono>
#include <iomanip>

. test function with a very simple add operation

void test_simple(long long int N, int * ary, double & sum, long long int & elapsed_milli)
{
  std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
  start = std::chrono::system_clock::now();
  double local_sum = 0.0;
  #pragma omp parallel
  {
    #pragma omp for reduction(+:local_sum)
    for (long long int i = 0; i < N; i++) {
      local_sum += ary[i];
    }
  }
  sum = local_sum;
  end = std::chrono::system_clock::now();
  elapsed_milli = std::chrono::duration_cast<std::chrono::microseconds>
                             (end-start).count();
}

. test function with a complex, CPU intensive operation sign(x) atan(sqrt(cos(x)^2 + sin(0.5x)^2)

void test_intensive(long long int N, int * ary, double & sum, long long int & elapsed_milli)
{
  std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
  start = std::chrono::system_clock::now();
  double local_sum = 0.0;
  #pragma omp parallel
  {
    double c, s;
    #pragma omp for reduction(+:local_sum)
    for (long long int i = 0; i < N; i++) {
      c = cos(double(ary[i]));
      s = sin(double(ary[i])*0.5);
      local_sum += atan(sqrt(c*c+s*s));
    }
  }
  sum = local_sum;
  end = std::chrono::system_clock::now();
  elapsed_milli = std::chrono::duration_cast<std::chrono::microseconds>
                             (end-start).count();  
}

. Main function

using namespace std;
int main() {
  long long int N = 1073741825,i;
  int * ary = new int[N];
  srand (0);
  for (i = 0; i < N; i++) { ary[i] = rand()-RAND_MAX/2; }
  double sum = 0.0;
  sum = 0.0;
  long long int  elapsed_milli;
  cout <<"#"<<setw(19)<<"N"<<setw(20)<<"µs"<< endl;
  for(i=128; i<N; i=i*2)
  {
      test_intensive(i, ary, sum, elapsed_milli);
      //test_simple(i, ary, sum, elapsed_milli);
      cout << setw(20)<<i<<setw(20)<<elapsed_milli << setw(20)<<sum<<endl;
  }
}

Compile (using icpc)
Sequential (No OpenMP) version is compiled with :

icpc test_omp.cpp -O3 --std=c++0x  

OpenMP (OpenMP) version is compiled with :

icpc test_omp.cpp -O3 --std=c++0x -openmp

Measurement
Time measurements are done with chrono using high_precision_clock and the limit precision on my machine is microseconds hence the use of std::chrono::microseconds (no point looking for higher precision)

Graph for the simple operation (axis are in log scale !)

Graph for the complex operation (axis are in log scale !)

Conclusions drawn

• There is an offset the first time using OpenMP (the first #pragma omp crossed) because the pool thread must be set in place.
  If we take a closer look at the 'intensive case' the first time we enter the test_ function (with i=128) the time cost is way higher in the OpenMP case than the in No OpenMP case. At the second call (with i=256) we dont see the benefit of using OpenMP but the timings are coherent.

• We can see that we do not observe scalability with a small number of samples. It's clearer in the simple test case. In other words the amount of operations inside a parallel section must be high enough to render the time needed for thread pool management negligable. Otherwise there is no point dividing the operation into threads.

• In this case (with the processor I used) the minimum number of samples is around 100000. But if I were to use 256 threads it would surely be around 6000000.

• However for more CPU intensive operations using OpenMP can induce speed up even with 1000 samples (with the processor I used)

Summary

• If you bench an OpenMP code try to set up the pool thread beforehand with a simple operation with #pragma omp parallel. In your test case the setting up takes most of the time.
• Using OpenMP is a catch only if you parallelize sufficently CPU intensive functions (which is not really the case of a simple array sum...). For example this is the reason why in nested loops the #pragma omp for should always be in the outermost "possible" loop.

Answer 3

As Max Langhof and user463035818 suggested, the program is memory bounded. I changed the program to do something more than accumulation. That is, I changed sum += ary[i] to sum += (pow(ary[i], 1.1) + pow(ary[i], 1.2)) / 100000000.0 and performed the same change in the parallel program and measured the time. The parallel program became 2X faster. If the program is IO bounded, I guess there is not much I can do about it to make it faster with OpenMP. Please let me know if you think otherwise.