openmp update matrix from neighboring elements (parallelise preconditioner)

Question

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...

Answer 1

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.

Answer 2

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).