An answer to this question on Stack Overflow.

Question

I would like to implement OpenMP to parallelize my code. I am starting from a very basic example to understand how it works, but I am missing something...

So, my example looks like this, without parallelization:

int main() {
  ...
  for (i = 0; i < n-1; i++) {
    u[i+1] = (1+h)*u[i]; // Euler
    v[i+1] = v[i]/(1-h); // implicit Euler
  }
    
  ...
  return 0;
}

Where I omitted some parts in the "..." because are not relevant. It works, and if I print the u[] and v[] arrays on a file, I get the expected results.

Now, if I try to parallelize it just by adding:

#include <omp.h>
int main() {
  ...
  omp_set_num_threads(2);
  #pragma omp parallel for
  for (i = 0; i < n-1; i++) {
    u[i+1] = (1+h)*u[i]; // Euler
    v[i+1] = v[i]/(1-h); // implicit Euler
  }
    
  ...
  return 0;
}

The code compiles and the program runs, BUT the u[] and v[] arrays are half full of zeros.

If I set omp_set_num_threads( 4 ), I get three quarters of zeros.
If I set omp_set_num_threads( 1 ), I get the expected result.

So it looks like only the first thread is being executed, while not the other ones...

What am I doing wrong?

Answer

OpenMP assumes that each iteration of a loop is independent of the others. When you write this:

for (i = 0; i < n-1; i++) {
  u[i+1] = (1+h)*u[i]; // Euler
  v[i+1] = v[i]/(1-h); // implicit Euler
}

The iteration i of the loop is modifying iteration i+1. Meanwhile, iteration i+1 might be happening at the same time.

Unless you can make the iterations independent, this isn't a good use-case for parallelism.

And, if you think about what Euler's method does, it should be obvious that it is not possible to parallelize the code you're working on in this way. Euler's method calculates the state of a system at time t+1 based on information at time t. Since you cannot knowing what's at t+1 without knowing first knowing t, there's no way to parallelize across the iterations of Euler's method.