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