11

行列 A を更新するループがあり、それを openmp にしたいのですが、どの変数を共有して非公開にするべきかわかりません。ii と jj だけでうまくいくと思っていましたが、うまくいきません。どこかで !$OMP ATOMIC UPDATE も必要だと思います...

ループは、N 個と N-1 個の粒子間の距離を計算し、行列 A を更新するだけです。

            !$OMP PARALLEL DO PRIVATE(ii,jj)
            do ii=1,N-1
                    do jj=ii+1,N
                            distance_vector=X(ii,:)-X(jj,:)
                            distance2=sum(distance_vector*distance_vector)
                            distance=DSQRT(distance2)
                            coff=distance*distance*distance
                            PE=PE-M(II)*M(JJ)/distance
                            A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
                            A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
                    end do
            end do
            !$OMP END PARALLEL DO
4

2 に答える 2

25

OpenMP の黄金律は、外側のスコープで定義されたすべての変数 (いくつかの除外を含む) が、デフォルトで並列領域で共有されることです。2008 年より前の Fortran にはローカル スコープがないため (つまりBLOCK ... END BLOCK、以前のバージョンにはありません)、すべての変数 (変数を除くthreadprivate) が共有されます。これは私にとって非常に自然なことです (Ian Bush とは異なりdefault(none)、さまざまな複雑な科学コードで 100 以上のローカル変数すべての可視性を再宣言します)。

各変数の共有クラスを決定する方法は次のとおりです。

  • N- すべてのスレッドで同じにする必要があり、その値を読み取るだけであるため、共有されます。
  • ii- ワークシェアリング ディレクティブの対象となるループのカウンターであるため、その共有クラスはprivate. PRIVATE句で明示的に宣言しても問題はありませんが、実際には必要ありません。
  • jjjj- ワークシェアリング ディレクティブの対象とならないループのループ カウンターは、したがってprivate.
  • X- すべてのスレッドが参照し、そこからのみ読み取るため、共有されます。
  • distance_vectorprivate-各スレッドが異なるペアのパーティクルで動作するため、明らかにそうあるべきです。
  • distance、distance2、およびcoff- 同上。
  • M- と同じ理由で共有する必要がありますX。
  • PE-アキュムレータ変数として機能し（これはシステムのポテンシャルエネルギーだと思います）、リダクション操作の対象になる必要があります。つまり、REDUCTION(+:....)節に入れる必要があります。
  • A-これはトリッキーです。これは共有さA(jj,:)れ、同期構造で保護されて更新されるか、リダクションを使用できます (OpenMP では、C/C++ とは異なり、Fortran では配列変数のリダクションが可能です)。A(ii,:)複数のスレッドによって変更されることはないため、特別な処理は必要ありません。

リダクションを適切に行うとA、各スレッドは のプライベート コピーを取得しA、これはメモリを大量に消費する可能性がありますが、この直接的な O(N 2 ) シミュレーション コードを使用して、非常に多数の粒子を含むシステムを計算するとは思えません。リダクションの実装に関連する特定のオーバーヘッドもあります。この場合、句Aのリストに追加するだけです。REDUCTION(+:...)

同期構造には 2 つのオプションがあります。ATOMICコンストラクトまたはコンストラクトのいずれかを使用できますCRITICAL。スカラー コンテキストにのみ適用されるためATOMIC、代入ループを「ベクトル化解除」し、ATOMIC各ステートメントに個別に適用する必要があります。次に例を示します。

!$OMP ATOMIC UPDATE
A(jj,1)=A(jj,1)+(M(ii)/coff)*(distance_vector(1))
!$OMP ATOMIC UPDATE
A(jj,2)=A(jj,2)+(M(ii)/coff)*(distance_vector(2))
!$OMP ATOMIC UPDATE
A(jj,3)=A(jj,3)+(M(ii)/coff)*(distance_vector(3))

これをループとして書き直すこともできます。ループ カウンターを宣言することを忘れないでくださいprivate。

CRITICALループをベクトル化解除する必要はありません。

!$OMP CRITICAL (forceloop)
A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
!$OMP END CRITICAL (forceloop)

重要な領域に名前を付けるのはオプションであり、この特定のケースでは少し不必要ですが、一般的には関連のない重要な領域を分離することができます。

どちらが速いですか？ATOMICまたはCRITICAL?で展開 それは多くのことに依存します。CRITICALアトミックインクリメントは、少なくとも x86 では、ロックされた加算命令で実装されますが、OpenMP ランタイムへの関数呼び出しが頻繁に含まれるため、通常はかなり遅くなります。彼らがよく言うように、YMMV。

要約すると、ループの作業バージョンは次のようになります。

!$OMP PARALLEL DO PRIVATE(jj,kk,distance_vector,distance2,distance,coff) &
!$OMP& REDUCTION(+:PE)
do ii=1,N-1
   do jj=ii+1,N
      distance_vector=X(ii,:)-X(jj,:)
      distance2=sum(distance_vector*distance_vector)
      distance=DSQRT(distance2)
      coff=distance*distance*distance
      PE=PE-M(II)*M(JJ)/distance
      do kk=1,3
         !$OMP ATOMIC UPDATE
         A(jj,kk)=A(jj,kk)+(M(ii)/coff)*(distance_vector(kk))
      end do
      A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
   end do
end do
!$OMP END PARALLEL DO

あなたのシステムは 3 次元であると仮定しました。

以上のことを踏まえて、私はイアン・ブッシュに、位置行列と加速度行列がメモリにどのように配置されるかを再考する必要があることを支持します。キャッシュを適切に使用すると、コードが向上し、X(:,ii)-X(:,jj)ベクトル SIMD 命令を使用して実装されるベクトル化などの特定の操作が可能になります。

于 2013-01-17T14:30:13.130 に答える
3

書かれているように、競合状態を避けるために同期が必要になります。2 スレッドの場合を考えてみましょう。スレッド 0 が ii=1 で始まるとすると、jj=2,3,4, .... と見なされ、スレッド 1 は ii=2 で始まると見なされ、jj=3,4,5,6 と見なされます。したがって、スレッド 0 が ii=1,jj=3 を検討し、スレッド 1 が ii=2,jj=3 を同時に検討している可能性があります。これは明らかにラインで問題を引き起こす可能性があります

                        A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)

両方のスレッドの jj の値が同じであるためです。そうです、競合を避けるために更新を A に同期する必要がありますが、良い方法がすぐにはわからないことを認めなければなりません。考えて、何かあったら編集します。

ただし、他に 3 つのコメントがあります。

1）あなたのメモリアクセスパターンはひどいです。これを修正すると、少なくともopenmpと同じくらい高速になり、面倒が少なくなると思います。Fortran では、最初のインデックスを最速で下る必要があります。これにより、メモリ アクセスが空間的にローカルになり、メモリ階層を適切に使用できるようになります。これが最新のマシンで良好なパフォーマンスを得るために最も重要なことであることを考えると、これを正しく行うように努める必要があります。したがって、上記が次のように記述できるように配列を配置できれば、上記のほうがよいでしょう。

        do ii=1,N-1
                do jj=ii+1,N
                        distance_vector=X(:,ii)-X(:jj)
                        distance2=sum(distance_vector*distance_vector)
                        distance=DSQRT(distance2)
                        coff=distance*distance*distance
                        PE=PE-M(II)*M(JJ)/distance
                        A(:,jj)=A(:,jj)+(M(ii)/coff)*(distance_vector)
                        A(:,ii)=A(:,ii)-(M(jj)/coff)*(distance_vector)
                end do
        end do

これが、2番目のインデックスではなく、最初のインデックスにどのように適用されるかに注意してください。

2) openmp を使用する場合は、default(None) を使用することを強くお勧めします。これは厄介なバグを回避するのに役立ちます。もしあなたが私の生徒の一人だったら、これをしないとたくさんの点数を失うでしょう!

3) Dsqrt は古風です - 現代の Fortran (つまり 1967 年以降のもの) では、いくつかのあいまいなケースを除いて、sqrt で十分であり、より柔軟です。

于 2013-01-16T15:15:50.183 に答える