1

I have the following code which I would like to make it parallel (pseudocode)

int na = 10000000;
int nb = na;
double A[na];
double B[2*na];
double a;

for(int j=0;j<nb;j++)
{
  i = rand() % na;
  A[i]+=5.0*i;
  B[i+10]+=6.0*i*i;
}

Of course, I cannot use #pragma omp parallel for because sometimes (which cannot be predicted) the same element will be accessed by two threads at the same time. How can this block of code be parallelized? Thank you

dimpep
  • 63
  • 1
  • 6
  • https://stackoverflow.com/questions/16789242/fill-histograms-array-reduction-in-parallel-with-openmp-without-using-a-critic – Z boson Jun 26 '18 at 13:11

1 Answers1

1

There are two ways to do that:

  • Use an atomic update on the values

    #pragma omp parallel for
for(int j=0;j<nb;j++)
{
    // make sure to declare i locally!
    int i = fun();
    #pragma omp atomic
    A[i]+=5.0*i;
}

    This is the simplest way. Each write is executed atomically and therefore more expensive. You also need to consider that accessing adjacent elements from multiple threads becomes expensive (false sharing). Use this if A is large and you do a lot of computations per uptate.

  • Use an array-reduction

    #pragma omp parallel for reduction(+:A)
for(int j=0;j<nb;j++)
{
    // make sure to declare i locally!
    int i = fun();
    A[i]+=5.0*i;
}

    This creates a local copy of A for each thread which is added together to the outside A after the parallel region. This requires more memory and some computation after, but parallel code itself can work most efficiently. Use this if A is small and the are little computations for each update.

BTW: Never use rand() in parallel applications, it is not defined as thread safe and sometimes it is implemented with a lock and becomes horribly inefficient.

EDIT: In your example with B you can safely apply either omp atomic or reduction separately to the statement since each operation only needs to be performance atomically independently.

Zulan
  • 21,896
  • 6
  • 49
  • 109
  • I don't use rand(). Rand is used here as an example to show what you have written as a fun() – dimpep Jun 25 '18 at 12:57
  • In addition, can I use `#pragma omp atomic` for more than one updates (such as `A[i]+=fun() B[i+5]+=fun()` etc? – dimpep Jun 25 '18 at 13:03
  • You cannot update two values in one atomic update - but if the updates don't depend on each other that may still be valid. That would require a specific example. – Zulan Jun 25 '18 at 13:09
  • Maybe 2 different `atomic` clauses? One for each value to be updated – dimpep Jun 25 '18 at 13:13
  • Please edit your question to include a specific example for two updates that shows the dependencies for the updates. – Zulan Jun 25 '18 at 13:33
  • @dimpep Yes, different `atomic` or `reduction` clauses are fine in your edited example. – Zulan Jun 25 '18 at 16:42
  • I don't know if I am thinking correctly today (bad cold) but do you really want to use an atomic clause in the array reduction? Shouldn't you just declare `reduction(+:A[na])` and remove the atomic cluase? – Z boson Jun 26 '18 at 13:09
  • @Zboson thanks for catching, of course you are right that was a copy-paste error. However, I think you don't need to specify `[na]` for the full array (only for array sections), do you? – Zulan Jun 26 '18 at 13:45
  • @Zulan, I don't know. I am going by this comment https://stackoverflow.com/questions/20413995/reducing-on-array-in-openmp/20421200#comment62676230_20413995. Fortran allowed array reductions for a long time but C did not because it could not determine the size so I suspect they fixed it by explicitly stating the size. – Z boson Jun 26 '18 at 13:49
  • I just checked the OpenMP 4.5 documentation. I think you are right that you don't need the array section notation I mentioned. I'm not sure why he used it [here](https://stackoverflow.com/questions/20413995/reducing-on-array-in-openmp/20421200#comment62676230_20413995). What if `A` is a pointer type and not an array? Does it have to be an array? I could not figure it out from the documentation. Maybe you need array sections when the type is a pointer? – Z boson Aug 06 '18 at 08:27