OPENMP - 前提条件付きシュワルツアルゴリズムを並列化する

Question

Schwarzアルゴリズムを並列化する必要がありますが、前提条件を処理する方法と入れ子になったループがあることはわかりません。 OpenMPまたはMPIを使用する必要があります。

void ssor_forward_sweep(int n, int i1, int i2, int j1, int j2, int k1, int k2, double* restrict Ax, double w) 
{ 
#define AX(i,j,k) (Ax[((k)*n+(j))*n+(i)]) 

int i, j, k; 
double xx, xn, xe, xu; 

for (k = k1; k < k2; ++k) { 
    for (j = j1; j < j2; ++j) { 
     for (i = i1; i < i2; ++i) { 
     xx = AX(i,j,k); 
     xn = (i > 0) ? AX(i-1,j,k) : 0; 
     xe = (j > 0) ? AX(i,j-1,k) : 0; 
     xu = (k > 0) ? AX(i,j,k-1) : 0; 
     AX(i,j,k) = (xx+xn+xe+xu)/6*w; 
     } 
    } 
} 
#undef AX 
} 

各ループが前回のループからの値を使用することを考慮して、この機能を並列化して最適な時間を得る方法。

私はすでに（ステンシルヤコビ3Dのような）ブロックで2または分割によりループ2の並列化を試みたが、成功せず...

はどうもありがとうございました！

Answer 1

ループ間データ依存関係は、ネストされたループで得られる並列度を制限します。

依存性のあるタスクを使用することができます。これは最も簡単な方法です。 OpenMPランタイム・ライブラリはスケジューリングを行い、アルゴリズムにのみ焦点を当てます。もう一つの良い点は、ループの最後には同期がなく、コードの依存する部分間だけに同期があることです。

#pragma omp parallel 
#pragma omp single 
for (int k = 0; k < k2; k += BLOCK_SIZE) { 
    for (int j = 0; j < j2; j += BLOCK_SIZE) { 
    for (int i = 0; i < i2; i += BLOCK_SIZE) { 
     #pragma omp task depend (in: AX(i-1,j,k),AX(i,j-1,k),AX(i,j,k-1)) \ 
         depend (out: AX(i,j,k)) 
     // your code here 
    } 
    } 
} 

タスクは、時には並列ループ（負荷と同期の粒度に応じて）よりも少し高価であるので、別の代替は、基本的には内側ループの要素が独立しているように、反復空間を変換wavefront parallelizationパターン、ですお互いに（そのためにあなたのために並列を使用することができます）間。

あなたはどちらの方法にしたい場合は、私は強くあなたがブロック1にあなたのアルゴリズムを有効にすることをお勧め：2段階に分けて計算を行うためにあなたの3 - ネストされたループを展開：

1. 反復固定サイズのブロック/キューブ（新しい誘導変数をii、jj、kkと呼ぶ）の中から選択します。
2. ブロックごとに、ループの元のシリアルバージョンを呼び出します。

ブロックする目的は、パラレル化のオーバーヘッドが目立たないように、パラレルな部分の粒度を増やすことです。ここ

が遮断部のためのいくつかの擬似コードである：波面並列の場合

#define min(a,b) ((a)<(b)?(a):(b)) 
// Inter block iterations 
for (int kk = 0; kk < k2; kk += BLOCK_SIZE) { 
    for (int jj = 0; jj < j2; jj += BLOCK_SIZE) { 
    for (int ii = 0; ii < i2; ii += BLOCK_SIZE) { 

     // Intra block iterations 
     for (int k = kk; k < min(k2,k+BLOCK_SIZE); k++) { 
     for (int j = jj; j < min(j2,j+BLOCK_SIZE); j++) { 
      for (int i = ii; i < min(i2,i+BLOCK_SIZE); i++) { 
       // Your code goes here 
      } 
     } 
     } 

    } 
    } 
} 

、最後のステップでは、そのように、波面にアウターループ（インターブロック反復）を旋回しています相互に依存しない要素を反復処理します。 3D反復空間では、基本的に（0,0,0）から（i2、j2、k2）に進む対角面です。下のイメージでは、赤で強調表示されているようなものです。

理解しやすいので、私は、2D波面の一例を置くつもりです。

#define min(a,b) ((a)<(b)?(a):(b)) 
#pragma omp parallel 
for (int d = 1; d < i2+j2; d++) { 
    int i = min(d,i2) - 1; 
    int j = 0; 
    // Iterations in the inner loop are independent 
    // Implicit thread barrier (synchronization) at the end of the loop 
    #pragma omp for 
    for (; i >= 0 && j < min(d,j2); i--, j++) { 
     // your code here 
    } 
}