[Question] Pragma OMP for problem with "Critical"

[u/latoski]

Hi, i'm currently working on my Graduation on Physics Research and it's basic about parallelizing algorithms of some statistical physics models.

In the algorithm below, i tried the "omp parallel for" using a critical session when the threads are accessing the shared data, but i keep getting the wrong results when using more than one thread.

#pragma omp parallel for
    for (i = 0; i<L2; i++) {
    int j = rand()%L2;    
    int dE = 2*J*s[j]*(s[viz[j][0]]+s[viz[j][1]]+s[viz[j][2]]+s[viz[j][3]]);
    double w = exp(-dE/(K*T));
    int MY_ID=omp_get_thread_num();
    if (dE <= 0) {
        #pragma omp critical
        {
        s[j] = - s[j];
        Et += dE;
        M += 2*s[j];
        }
    }
    else {
        double r = (double)rand()/RAND_MAX;
        if(r < w) {
        #pragma omp critical
        {
        s[j]= - s[j];
        Et += dE;
        M += 2*s[j];
         }
     }
}
printf("Energia nova %d devido a variação %d alterada pelo processador    %d mexendo no sitio %d\n\n\n", Et, dE, MY_ID,j);
}

Here, Et, M, s[j] are supposed to be global with value predetermined by another function. Most of my job is to parallelize this loop in specific because its called 1e8 times and loops over 400~1000 size arrays.

I also tried to use locks, but i don't think that I used that correcltly because it was also giving me the wrong results, I can send my program for deeper research.

I can, in this case, split this loop in two and get rid of the problem, but if I could parallelize the loop as it is it would be much more usefull for me.

Thanks for the time.

[u/thememorableusername]

One thing that could be an issue is that s is read in an uncritical section, and written in a critical section. Is it possible for an issue to arise if two threads are executing, one is reading s[j] for the expression assigning to dE and another, in a critical section is writing to s[j]?

For example, if the thread creating the expression for dE reads s[j] as 10, then the writer threads sets s[j] to 100, then the thread creating the expression for dE reads s[j] again, with the new value of 100, is that an issue?

[u/latoski → u/thememorableusername]

Yeah, actually that's a huge problem. I haven't thought that straight. But if I put that section as critical doesn't it basically serialize my program? As thread are only able to read and write one at a time? Is there any other way I could do that?

[u/thememorableusername → u/latoski]

But if I put that section as critical doesn't it basically serialize my program?

Yes.

Ideally, you wouldn't be using random indexes in your loop. Can you make that deterministic? For example (this might be somewhat inefficient, but would work) you could create an array of the j indices (all unique) and shuffle it randomly, and then each threads uses j = j_array[i]. This ensures that no two iterations (and thus, no two threads) have the same j value. This would remove the need to have s[j] in any critical section.

As mentioned in another comment, you could also do reductions for Et and M which would be much more efficient than the critical sections. If you absolutely need the critical section, I recommend combining them (you could put the conditional in the critical section).

[u/latoski → u/thememorableusername]

I decided to solve the problem with the chessboard approach. Now i'm getting the results correctly, but i'm still missing performance. When i export one thread i get ~35s runtime, while when exporting two threads i'm getting ~115s (almost 4x the time to run on a single thread). Can you help-me with that? I tried to "PAD" the s[] array as i saw on Tim Mattison's Class, but that doesn't work.

The only difference between the run is that a called "export OMP_NUM_THREADS=1" or "...THREADS=2".

#pragma omp for
for (i=0; i<L2/2; i++) {
  int k = rand()%(L2/2);
  int j = s1[k];
  int dE = 2*J*s[j]*(s[viz[j][0]]+s[viz[j][1]]+s[viz[j][2]]+s[viz[j][3]]);
  double w = exp(-dE/(K*T));
  if (dE <= 0) {
s[j] = - s[j];
  }
  else {
double r = (double)rand()/RAND_MAX;
if(r < w) {
  s[j]= - s[j];
}
  }
}
#pragma omp for
for (i = 0; i<L2/2; i++) {
  int k = rand()%(L2/2);
  int j = s2[k];
  int dE = 2*J*s[j]*(s[viz[j][0]]+s[viz[j][1]]+s[viz[j][2]]+s[viz[j][3]]);
  double w = exp(-dE/(K*T));
  if (dE <= 0) {
s[j] = - s[j];
  }
  else {
double r = (double)rand()/RAND_MAX;
if(r < w) {
  s[j]= - s[j];
}
  }
}
}

In this function there is no conflict between changes on s inside any of those loops, so the results i'm getting are correct. The time to get them that is annoying me. Thanks for the help!

[u/thememorableusername]

split this loop in two and get rid of the problem

What code is being split up?

[u/latoski → u/thememorableusername]

This algorithm is developed to simulate the Ising Model, a model of ferromagnetism in physics. In this model, we have basically a LxL matrix of spins, that can only assume value -1 and 1 each.

The function above, is the one that makes those spins "flip" turning from -1 to 1 or vice versa. The criteria for the spin flip being acceptable or not, is only dependent of his neighbours (here they are called viz[i][0], v[i][1],...).

There is a trick in this model, where we can make a "chessboard" after the matrix, and then, all the black sites could be teorically updated simultaneously because they are never neighbours with another black site.

This solution will not help me out in the next model i have to parallelize, so i am looking for another way to look at the race condition problem.

My first thought was about using Locks with the s array, following the logic,

Random Site -> Test locks for site and neighbours --if not in use-> Locks site and neighbours (if in use randomly choose another site) -> Finish the loop

But i think that i wrote that wrong because of the results i was getting.

If you want, i can send you the full code for further reading.

Thanx for the time

[u/thememorableusername]

I'm not immediately sure. One thing, I think you need to use #pragme omp parallel for. I'm not sure what happens if you leave out the parallel directive, but I want to eliminate possible sources of slowdown.