動的配列を持つopenaccネストループ

Question

openaccを適用して、マルチコアおよびgpuアクセラレーションされたバイナリを開発しようとしています。私はFarberの本を読んで、そこから、そしてNVIDIAが提供するいくつかのオンラインコースを通して、正常にテストプログラムを走らせました。それから、私は従来のコードを並列化しようとしました。

私が書いているサブルーチンは、3つのネストされたループ内で使用される3次元配列に依存しています。これは、nstepsの軌跡に対する典型的なN（N + 1）/ 2対の賢明距離問題である。私たちの科学の問題では、nstepsは通常1E5〜1E7、n粒子は1E4〜5E5です。

for(i=0 ; i < nsteps ; i++){ 
    pair = 0 ; 
    for(j=0 ; j < nparticles-1 ; j++){ 
     x1 = position[i][j][0] ; 
     y1 = position[i][j][1] ; 
     z1 = position[i][j][2] ; 
     for(k=j+1 ; k < nparticles ; k++){ 
      x2 = position[i][k][0] ; 
      y2 = position[i][k][1] ; 
      z2 = position[i][k][2] ; 
      dx2 = (x1 - x2) * (x1 - x2) ; 
      dy2 = (y1 - y2) * (y1 - y2) ; 
      dz2 = (z1 - z2) * (z1 - z2) ; 
      sdist = sqrt(dx2 + dy2 + dz2) ; 
      dist[pair] += sdist ; 
      pair++ ; 
     } 
    } 
} 

Iコードとしてコンパイル時に入力位置のアレイ（nsteps、nParticleの）を制御することはできませんが、C-配列データ型にPythonのnumpyの配列を変換するC++拡張にパイソンを介して実行されます。 openaccコードはソースオブジェクトライブラリとしてコンパイルされます。 C++エクステンションはopenaccコードを呼び出します。これは、この配置を必要とするレガシーコードです。 Pythonでシステムを定義し、C++拡張を呼び出してopenacc * .soファイルにアクセスする手順はうまくいきます。

この問題は本やコースの演習では容易には分かりません。そのため、問題の解決策やコメントが他の人に役立つかもしれません。

スタックオーバーフローやその他の情報源の例を使用してコードスニペットにパラレルのカーネルとループの&データ指示文を試しました。私の知る限りでは、Faberの本やその他の情報源で使用されている例は、この問題で提起されたユースケースに対処していません。私は最初は2つのループを並列化することができますが、最も内側のループは並列化されません（ループはベクトル化されません：データ依存性）。私は一般的なガイダンスを探しているので、失敗したビットを投稿していないので、より教育的なディスカッションやヒントが提供されます。

さて、今私の質問に。

1. プラグマディレクティブを使用して入力位置配列の未知の次元を処理するにはどうすればよいですか？
2. dist []配列の累積を管理するにはどうすればいいですか？（パーティクルの数によってコンパイル時に不明な長さがあります）
3. この問題にはどのようなプラグマディレクティブが推奨されていますか？
4. "j-loop"の "k-loop"の依存関係をどのように扱いますか？
5. 使用するディレクティブを定義するのに役立つように問題を平坦化する必要がありますか？

THX、

SB

はUPDATE：私はDIST配列のインデックス付けを考慮するためにコードを変更し、提案プラグマを添加@jefflarkinの提案あたりの結果を提供する

。サブルーチンは正常にコンパイルされ、並列に実行されます。ここでは、プロファイリングを開始して、リソース使用率を最大化するためにルーチンを最適化する方法を見ていきます。以下に、作業コードのコピーを示します。

#pragma acc data copyin(position[nsteps][nparticles][3]) copy(dist[npairs]) 
for(i=0 ; i < nsteps ; i++){ 
    #pragma acc parallel loop 
    for(j=0 ; j < nparticles-1 ; j++){ 
     x1 = position[i][j][0] ; 
     y1 = position[i][j][1] ; 
     z1 = position[i][j][2] ; 
     #pragma acc loop 
     for(k=j+1 ; k < nparticles ; k++){ 
      x2 = position[i][k][0] ; 
      y2 = position[i][k][1] ; 
      z2 = position[i][k][2] ; 
      dx2 = (x1 - x2) * (x1 - x2) ; 
      dy2 = (y1 - y2) * (y1 - y2) ; 
      dz2 = (z1 - z2) * (z1 - z2) ; 
      sdist = sqrt(dx2 + dy2 + dz2) ; 
      local_count = ((j*nparticles)-((j*(j+1))/2))+k-(j+1) ; 
      dist[local_count] += sdist ; 
     } 
    } 
} 

このコンパイラ結果（pgC++10）：CFLAGS = -fPIC -C -fast -acc -Minfo =アクセル-ta =マルチコア-O3 -shared）

38, Generating Multicore code 
39, #pragma acc loop gang 
45, Loop is parallelizable 

及びGPUのための（CFLAGS = -v -fPIC -C -fast -acc - MINFO =アクセル-ta =テスラ：cuda8、fastmath -03 -shared）

35, Generating copyin(coor[:nframes][:nparticles][:3]) 
    Generating copy(dist[:npairs]) 
38, Accelerator kernel generated 
    Generating Tesla code 
    39, #pragma acc loop gang /* blockIdx.x */ 
    45, #pragma acc loop vector(128) /* threadIdx.x */ 
45, Loop is parallelizable 

ライン35、I-ループ（nsteps）、線38は、jループ線45はK-ループです。

Answer 1

あなたは、配列のサイズは不明ですが、n個のパーティクルに基づいてデータ句のサイズを計算することができます。位置配列はサイズ[n個] [n個] [3]です。 kは各jについて縮小するのでdist配列は少しトリッキーですが、（n個の粒子-1）+（n個の粒子2）+ ... + 1の合計です。 ）/ 2となる。パーティクルの数は実行時まで決定されないかもしれませんが、コンパイラは配列に関する何かを知っていなければなりません。それ以外の場合、[]は動作しません。

jとkに基づいて配列にインデックスを計算できるように、dist配列へのインデックスを変更することをお勧めします。そうしないと、ペアの増分を保護するためのアトミック操作が必要になり、各ステップについて、同じ一対の点が距離配列内の同じ場所にあることを保証する。

その固定、私は、ステップループにjのループ上acc parallel loopをacc dataを考えると、k個のループ上の簡単なacc loop思いでは、（依存関係の分析でコンパイラを助けるために）あなたを得るべきだとランニング。