OpenMP並列プログラミング

提供: floatingexception
移動: 案内、 検索

並列化できるループと並列化できないループ

原理的に並列化できないループ

ループの前回の結果を利用して、次の結果を求めるループは、当然のことながら、並列化できません。これをループの依存性(dependency)といいます。例えば、次のフィボナッチ数列を求めるプログラムは、並列化できません。 

     f(1) = 1
     f(2) = 1
     do i = 3, 100
        f(i) = f(i-1) + f(i-2)
     enddo

プログラムが複雑過ぎたり、副プログラム呼び出しがあったりすると、依存性の有無をコンパイラが判断できなくて、自動並列化できないことがあります。 vastfomp(後述)を使う場合には、CVD$ NODEPCHKというディレクティブ(指示行)をdo文の直前に書くと、ループの中に依存性の疑いがあっても、プログラマの責任で自動並列化されます。

次のプログラム断片は、潜在的な誤りのために、自動並列化できない例です。 

     do i = 1, n
        if (x(i) < 0.0d+0) then
           y = -1.0d+0
        else if (x(i) > 0.0d+0) then
           y = 1.0d+0
        endif
        z(i) = y
     enddo

IFブロックのどれも成立しない場合に、ループの中で変数yへの代入が行なわれず、ループの前回の値が使われ、依存性が発生します。 

     do i = 1, n
        if (x(i) < 0.0d+0) then
           y = -1.0d+0
        else if (x(i) >= 0.0d+0) then
           y = 1.0d+0
        endif
        z(i) = y
     enddo

このプログラムは、IFブロックのどちからかが必ず実行されるはずですが、そこまではvastfompが判定できずに、自動並列化できません。 

     do i = 1, n
        if (x(i) < 0.0d+0) then
           y = -1.0d+0
        else if (x(i) >= 0.0d+0) then
           y = 1.0d+0
        else
           y = 0.0d+0 ! To optimize
        endif
        z(i) = y
     enddo

こうすれば、自動並列化できます。

vastfompのリスティングファイル等を読み、ループを並列化できなかった理由を調べると、潜在的な誤りに気づくことがあります。

並列呼び出しできる副プログラムとできない副プログラム

この文書で、Fortran言語のサブルーチンと関数、 C/C++言語の関数を、まとめて副プログラムと呼びます。

副プログラムには、並列プログラムの中で同時に呼び出されると、正常に動くものと、異状を起こすものがあります。副プログラムが正常に並列に動くことを、 スレッドセーフ(thread-safe)といいます。

副プログラムが正常に並列に動くためには、副プログラムの局所変数が並列数だけ別々の インスタンス(変数のメモリー上の実体)を持っている必要があります。また、孫副プログラムも並列に動くようになっている必要があります。入出力を伴う副プログラムを並列に動かすのは困難です。

C/C++言語では、全ての局所変数をautoにして、広域変数の書き換えを行なわなければ、比較的簡単にスレッドセーフな副プログラムを作れます。

Fortran言語では、コンパイラオプションによって局所変数を動的に確保し、共通ブロックの内容の書き換えを避ければ、スレッドセーフな副プログラムを作れますが、 C/C++言語の場合ほど簡単ではありません。

なお、OpenMPのthreadprivateという機能を使うと、共通ブロックのコピーをスレッドの数だけ作って、共通ブロックを書き換えるプログラムを並列化できることがあります。

副プログラム呼び出しを含むループの並列化

前述のように、副プログラムには、スレッドセーフなものと、そうでないものがあるので、副プログラム呼び出しを行なうループは自動並列化できません。

ただし、SQRT等のFortranの標準組込み関数は、スレッドセーフであるとコンパイラが知っているので、ループの中にあってもそのループを自動並列化できます。

vastfompを使う場合には、CVD$ CNCALLというディレクティブ(指示行)をdo文の直前に書くと、ループの中に副プログラム呼び出しがあっても、プログラマの責任で自動並列化されます。例を示します。 

CVD$ CNCALL
     do i = 1, N
        y(i) = my_func(x(i))
     enddo

但し、副プログラムの引数になっている変数を、 privateとすべきかsharedとすべきか、 vastfompが判断を誤ることがありますので、ご注意ください。次の例の変数zのように、副プログラムによって値が変えられるので、 privateとすべきなのに、誤ってsharedにされることがありました。 

CVD$ CNCALL
     do i = 1, N
        call my_sub(x(i), z) ! サブルーチンがzを書き換え
        y(i) = z
     enddo

並列領域でのエラー検出

ループの中で実行時エラーを検出して、エラーメッセージを表示して終了するプログラムがしばしばあります。このようなプログラムは、そのままでは自動並列化できません。 

     do j = 1, N
        do i = 1, N
           if (a(i,j) <= 0.0d+0) then
              write(*,*) \"a(i,j) must be positive.\"
              stop
           endif
           ...
        enddo
     enddo

エラーが発生したかどうかのフラグか、エラー回数のカウンターを使うと、並列化できることがあります。 

     numerror = 0
!$omp parallel do reduction(+:numerror)
     do j = 1, N
        do i = 1, N
           if (a(i,j) <= 0.0d+0) then
              numerror = numerror + 1
           endif
           ...
        enddo
     enddo
     if (0 < numerror) then
        write(*,*) \"a(i,j) must be positive.\"
        stop
     endif

三角行列と対称行列の添え字

三角行列と対称行列の処理では、しばしば次のようなコードが現れます。 

     k = 0
     do i = 1, n
       do j = 1, i
         k = k + 1
         ...
       enddo
     enddo

このループは、変数kの依存性のために並列化できません。次のように書き換えると、並列化できることがあります。 

!$omp parallel do private(k)
     do i = 1, n
       do j = 1, i
         k = i * (i - 1) / 2 + j
         ...
       enddo
     enddo

Fortranプログラム自動並列化ツールvastfompは、このように自動的に変換してくれることもあります。

青山さんの経験則

理化学研究所情報基盤センターの青山さんから教えていただいた経験則です。正攻法ではありませんが、次の経験則があります。

逆むきに回しても結果が変わらないループは、たいてい並列化できる。

例えば、 

     do i = 1, n

というループがあった場合に、 

     do i = n, 1, -1

と逆向きに回してみて、計算結果が変わらなければ、そのループを並列化できる可能性が高いです。実際に、青山さんは、他人が書いたプログラムを並列化する際に、この方法で目星をつけているそうです。

OpenMPプログラムの高速化

多次元配列における偽共有現象

現代のCPU、だいたいPentium3以降では、 CPUの内部クロックが極めて高速なのに対して、メモリーの速度がそれほど速くないので、キャッシュメモリーによってメモリーの遅さを補っています。

メモリー共有型マルチCPUコンピューターでは、ハードウェアとOSが、キャッシュの整合性を保障するのが、重要な働きです。

キャッシュは、1バイトごとに管理されているのではなく、 キャッシュラインというブロックごとに管理されています。キャッシュラインの大きさはCPUによって異なり、例えばItanium2では128バイトです。

Fortranで2次元配列、例えばa(2,3)を宣言すると、次のようにメモリーに配置されます。 

a(1,1)
a(2,1)
a(1,2)
a(2,2)
a(1,3)
a(2,3)

左側の添字を1づつ増やすと、メモリーの連続する場所にアクセスできます。

次のプログラムは悪い例です。 

!$omp parallel do
     do i = 1, m
        do j = 1, n
           a(i, j) = f(i, j)
        enddo
     enddo

1番目のCPUがa(1,1)に書き込むのと同時に、 2番目のCPUがa(2,1)に書き込もうとします。これら2つの配列要素は、メモリー内で隣接し、必ずではありませんが同じキャッシュラインに含まれるので、本当に同時に書き込みが起きると、キャッシュの不整合が発生します。そこで、処理系が、CPUの書き込みを待たせて、キャッシュの整合性を保つための処理が発生します。

このようなオーバーヘッドを、 偽共有(false shareing)といいます。自動並列化またはOpenMPによって直列実行よりも遅くなった場合には、まずこの現象を疑うべきです。

ファイル:Falsesharering.gif

次のプログラムは、良い例です。 

!$omp parallel do
     do j = 1, n
        do i = 1, m
           a(i, j) = f(i, j)
        enddo
     enddo

一般的に、Fortranのプログラムでは、最も左側の添え字を最も内側のループで回すと、キャッシュの利用効率が上がって、処理が高速になります。

なお、C/C++言語では、多次元配列の配置がFortranと逆なので、最も右側の添え字を最も内側のループで回すと、処理が高速になります。

転置行列を計算する場合には、次の2種類のコードが考えられます。 

!$omp parallel do
     do j = 1, n
        do i = 1, n
           b(j,i) = a(i,j)
        enddo
    enddo



!$omp parallel do
     do j = 1, n
        do i = 1, n
           b(i,j) = a(j,i)
        enddo
    enddo

偽共有現象は、メモリーに書き込む際に問題になるので、後者のプログラム例のほうがより高速に動きます。

ループの合併

次のコード断片は、自乗共役勾配(Conjugate Gradient Square)法の一部分です。 

!$omp parallel do
     do i = 1, n
        e(i) = r(i) + beta * h(i)
     enddo
!$omp parallel do
     do i = 1, n
        p(i) = e(i) + beta * (h(i) + beta * p(i))
     enddo

ベクトル計算機では、このままでよいのですが、並列計算機では、次のようにループを合併すると、並列化領域の開始と終了のオーバーヘッドが減り、キャッシュミスも減り、より高速になります。 

!$omp parallel do
     do i = 1, n
        e(i) = r(i) + beta * h(i)
        p(i) = e(i) + beta * (h(i) + beta * p(i))
     enddo

この例ではこれ以上改良できませんが、もし、このループのあとでe(i)の値が不要ならば、 eを配列変数からスカラー変数に変えて、より高速にできます。

多重ループとリストベクトル

ベクトルコンピューターにとっては、多重ループよりも1重ループのほうが都合がよいので、次のプログラム断片のようなリストベクトルがしばしば使われます。 

     real(8) :: x(maxi, maxj, maxk) ! 3次元のデータ
     integer :: i, j, k, l
     integer :: lv(maxi * maxj * maxk, 3) ! これがリストベクトル

     ! リストベクトルを作る
     l = 0
     do k = 1, maxk
       do j = 1, maxj
         do i = 1, maxi
           l = l + 1
           lv(l, 1) = i
           lv(l, 2) = j
           lv(l, 3) = k
         enddo
       enddo
     enddo
     ! リストベクトルを使って1重ループで処理する
!$omp parallel do private(i, j, k)
     do l = 1, maxi * maxj * maxk
       i = lv(l, 1)
       j = lv(l, 2)
       k = lv(l, 3)
       ! x(i, j, k) を処理する
       ...
     enddo

しかし、並列計算機では、次のプログラム断片のように、多重ループを使って最も外側のループを並列化すると、より高速になります。 

!$omp parallel do
     do k = 1, maxk
       do j = 1, maxj
         do i = 1, maxi
           ! x(i, j, k) を処理する
           ...
         enddo
       enddo
     enddo

ここでも、Fortran言語の原則で、左側の添字を内側のループで処理すると効率がよくなります。


OpenMP 2.5の新機能 workshare

OpenMP 2.5にworkshareという新機能があります。並列化ブロック内の配列演算を自動並列化します。例を示します。 

!$omp parallel workshare
     a(:,:) = b(:,:)
!$omp end parallel workshare


single構文の活用

並列化できる部分とできない部分が混在しているコードでは、 しばしば次のように書きます。 

!$omp parallel
      並列化できる部分
!$omp end parallel
      並列化できない部分
!$omp parallel
      並列化できる部分
!$omp end parallel

しかし、次のように書くほうが、スレッドの生成消滅のオーバーヘッドが減って、 速くなると、参考書に書いてあります。 アプリケーションプログラムのコードでは、あまり変わらないこともあります。 

!$omp parallel
      並列化できる部分
!$omp single
      並列化できない部分
!$omp end single
      並列化できる部分
!$omp end parallel

reductionのタイミング

前述の「single」構文を実験している際に、 自乗共役勾配法が収束しなくなるという問題が発生し、 いつリダクションが起きるかが問題になりました。

次の例は、 数列の和を求めるプログラムで、最初の並列ループでは1から10の和、 次の並列ループでは1から100の和、次の並列ループは1から1000の和を 計算しています。 

      program reduction1
      implicit none
      integer :: a(1000), i, s1, s2, s3

      do i = 1, 1000
         a(i) = i
      enddo
      s1 = 0
      s2 = 0
      s3 = 0
!$omp parallel reduction(+:s1,s2,s3)
!$omp do
      do i = 1, 10
         s1 = s1 + a(i)
      enddo
!$omp single
      write(*,*) s1, s2, s3
!$omp end single
!$omp do
      do i = 1, 100
         s2 = s2 + a(i)
      enddo
!$omp single
      write(*,*) s1, s2, s3
!$omp end single
!$omp do
      do i = 1, 1000
         s3 = s3 + a(i)
      enddo
!$omp end parallel
      write(*,*) s1, s2, s3
      stop
      end program reduction1

Intel Fortran Compiler 9.1 ia64での実行結果 

OMP_NUM_THREADS=2
          15           0           0
          15        1275           0
          55        5050      500500

「end parallel」のタイミングで、正しい合計が求まっています。

「reduction」と書く場所を変えてみます。 

      program reduction2
      implicit none
      integer :: a(1000), i, s1, s2, s3

      do i = 1, 1000
         a(i) = i
      enddo
      s1 = 0
      s2 = 0
      s3 = 0
!$omp parallel
!$omp do reduction(+:s1)
      do i = 1, 10
         s1 = s1 + a(i)
      enddo
!$omp single
      write(*,*) s1, s2, s3
!$omp end single
!$omp do reduction(+:s2)
      do i = 1, 100
         s2 = s2 + a(i)
      enddo
!$omp single
     write(*,*) s1, s2, s3
!$omp end single
!$omp do reduction(+:s3)
      do i = 1, 1000
         s3 = s3 + a(i)
      enddo
!$omp end parallel
      write(*,*) s1, s2, s3
      stop
      end program reduction1

OMP_NUM_THREADS=2
          55           0           0
          55        5050           0
          55        5050      500500

「do」が終わったタイミングで、正しい合計が求まっています。 PGIでも同様の結果が出ました。

上記の挙動が、仕様なのか、実験した環境ででたまたまこうなるのか、 調査中です。 Open-MPの仕様書には、次のように書かれています。

The reduction clause specifies an operator and one or more list items. For each list item, a private copy is created on each thread, and is initialized appropriately for the operator. After the end of the region, the original list item is updated with the values of the private copies using the specified operator.

「region」の終わり、すなわち「end parallel」のタイミングで、 リダクションされた値が保障されるようです。

関連項目

  • VAST-F/toOpenMP --- FortranのソースコードにOpenMPディレクティブを自動的に挿入するツール

参考書

  • 牛島省、
    「OpenMPによる並列プログラミングと数値計算法」、
    丸善、2006年、ISBN4-621-07717-1 アマゾンで詳細を見る
    (プログラム例はFortranで書かれています。入門から始まって、LU分解の並列化や差分方程式の並列化のような例が豊富です。)

外部リンク

「http://granada.sakura.ne.jp/floatingexception/index.php?title=OpenMP並列プログラミング&oldid=20」から取得