3
$\begingroup$

Issue of data dependency with stencil code...

How to parallelise this using openmp? I looked at the openmp manual, I figured out how to use DO ORDERED to get the same result as the serial version, but when run with openmp, it is x200 or so times slower.

This is example code from the internet doing the update from neighboring cells just to demonstrate the kind of problem I am after. I am just looking for a conceptual description of the correct approach, with examples if you like, but you don't have to spend time on implementing for this specific case.

  // Computation
  for ( i=1 ; i < N-1 ; i++ ) {
     for ( j=1 ; j < N-1 ; j++ ) {
        prev(i,j) = ((prev(i-1,j-1)+foo*prev(i-1,j)+bar*prev(i-1,j+1)+prev(i,j-1)+prev(i,j)+prev(i,j+1)+prev(i+1,j-1)+prev(i+1,j)+prev(i+1,j+1)));
     }
  }

You can also see an example of another similar problem here, on page 6.

And another example from here (in C):

inline void NavierCalc::PSOR(int i,int j,int k,NS_REAL omg1,NS_REAL lomg) {

// Relaxiertes Gaus-Seidel-Verfahren
//
IFFLUID(flag[i][j][k]) {
SetPLocalBorder(P,i,j,k);
P[i][j][k]=omg1*P[i][j][k] - lomg/(S.ddPstar[0][3][i]+
     S.ddPstar[1][4][j]+S.ddPstar[2][5][k]) *  
  (S.ddPstar[0][6][i]*P[i+1][j][k] + S.ddPstar[0][0][i]*P[i-1][j][k] + 
   S.ddPstar[1][7][j]*P[i][j+1][k] + S.ddPstar[1][0][j]*P[i][j-1][k] +
   S.ddPstar[2][8][k]*P[i][j][k+1] + S.ddPstar[2][0][k]*P[i][j][k-1] 
   - RHS[i][j][k]);
}
}

A relevant example

From this page, except in my case it is not PX, it is X at the left-hand side, so it does depend on itself, and the loop is split into two (one walks from left to right and uses SOMETHING-1 indices at the right-hand side, while the other walks from right to left and uses SOMETHING+1 indices at the right-hand side): 

 do j=2,size(X,dim=2)-1
 do i=2,size(X,dim=1)-1
    PX(i,j,bid) = precondNE    (i,j,bid)*X(i+1,j+1,bid) + &
                  precondNW    (i,j,bid)*X(i-1,j+1,bid) + &
                  precondSE    (i,j,bid)*X(i+1,j-1,bid) + &
                  precondSW    (i,j,bid)*X(i-1,j-1,bid) + &
                  precondNorth (i,j,bid)*X(i  ,j+1,bid) + &
                  precondSouth (i,j,bid)*X(i  ,j-1,bid) + &
                  precondEast  (i,j,bid)*X(i+1,j  ,bid) + &
                  precondWest  (i,j,bid)*X(i-1,j  ,bid) + &
                  precondCenter(i,j,bid)*X(i  ,j  ,bid)
 end do
 end do

One more relevant example

this page - SSOR preconditioner:

!   Execute SSOR preconditioner.
    do j= 1,n-1
      do i = 1,n-1
        rhat(i,j) = w*(r(i,j)-aw(i,j)*rhat(i-1,j)   &
                     -as(i,j)*rhat(i,j-1))/ac(i,j)
      end do
    end do
    do j= 1,n-1
      do i = 1,n-1
        rhat(i,j) =  ((2.-w)/w)*ac(i,j)*rhat(i,j)
      end do
    end do
    do j= n-1,1,-1
      do i = n-1,1,-1
        rhat(i,j) = w*(rhat(i,j)-ae(i,j)*rhat(i+1,j)  &
                      -an(i,j)*rhat(i,j+1))/ac(i,j)
      end do
    end do

And this presentation outlines the concept or reordering by coloring. More examples may be needed to understand the concept fully...

Improve this question
$\endgroup$
3
$\begingroup$

Since you're traversing your mesh left to right, top to bottom, most steps will depend on results computed in the directly previous step. If you want to maintain the order, this leaves very little room for parallelism. The best you could do is to get each thread to do a row of your matrix, i.e. parallelize only over i, and somehow try to prevent the threads overtaking each other.

In general, synchronization points, i.e. parts of the parallel code where one thread waits for another, will always have a negative impact on parallel performance, and should be avoided.

You may therefore want to ask yourself how crucial it is that your parallelized code replicates the serial results exactly? Such stencil updates are usually used as part of a Gauss-Seidel iteration, or any other iterative method, which should converge to the correct result and is usually not dependent on the exact order of the updates.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The block in question appears to involve preconditioner, so the order matters (otherwise the solver takes a very long time to converge). What examples of not-overtaking threads could you find? $\endgroup$ – user5273 Oct 27 '13 at 22:51
  • $\begingroup$ @Svetlana That would depend on how your preconditioner works, i.e. Is it sensitive to the order in which the variables are solved for? In any case, I would strongly advise against synchronizing the threads such that they don't overtake each other as this will probably cause massive overheads and bad scaling. $\endgroup$ – Pedro Oct 28 '13 at 9:34
  • $\begingroup$ It is sensitive to the order. Agreed "not overtaking" is a poor idea. What else to try? $\endgroup$ – user5273 Nov 6 '13 at 6:16
  • $\begingroup$ (Note: Question is updated with relevant snippets.) $\endgroup$ – user5273 Nov 6 '13 at 6:16
3
$\begingroup$

I think you're overlooking a major point in parallelizing your code, and that is the issue of data dependencies. All of the examples you've given are of the form

for i=1..N
   for j=1..N
      A(i,j) = foo*A(i-1,j-1) + foo*A(i-1,j+1) + foo*A(i+1,j+1) + ...

This won't parallelize since updating A depends on other loop iterations. One approach to dealing with the dependency would be to use temporary arrays. For example:

#pragma omp parallel
   #pragma omp for
      for i=1..N
         for j=1..N
            temp(i,j) = foo*A(i-1,j-1) + foo*A(i-1,j+1) + foo*A(i+1,j+1) + ...

   #pragma omp for
      for i=1..N
         for j=1..N
            A(i,j) = temp(i,j)

Besides restructuring your algorithm, I would recommend using some kind of reporting tool with your compiler to find what is being vectorized. Intel's compilers use -vec-report, and I believe GNU has -opt-report (double-check that).

Improve this answer
$\endgroup$
2
  • $\begingroup$ I'm not 'overlooking' the point. I put it into the question title. $\endgroup$ – user5273 Oct 29 '13 at 22:09
  • 1
    $\begingroup$ When calculating temp(i,j), it NEEDS to know fresh A(i-1,j-1). Not the old value. $\endgroup$ – user5273 Nov 3 '13 at 22:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy