バイトニックソートの並列処理

前回、並列処理の実験として奇偶転置ソートを実装しました。今回は、平均アルゴリズム計算量がO(N・log^2(N))で高速なバイトニックソート(Bitonic sort)を実装しました。

前回の続きの記事なので、重複する解説は省きます。

1. アルゴリズム

実はよくわかっていません(:P

バイトニックソートのベースは(たぶん)マージソートで、分割と統治の繰り返しを行っています。


// 基本的なバイトニックソート
void bitonicSort( int num, int *ary1 )
{
  // 2の冪乗でなければ計算しない(5)
  if( ( ( num - 1 ) & num ) != 0 )
     return;

  for( int block = 2; block <= num; block *= 2 ) {
    for( int step = block / 2; step >= 1; step /= 2 ) {
      for( int idx = 0; idx < num; idx++ ) {
        int e = idx ^ step; // (1)
        if( e > idx ) { // (2)
          int v1 = ary1[ idx ];
          int v2 = ary1[ e ];
          if( ( idx & block ) != 0 ) { // (3)
            if( v1 < v2 ) { // (4)
              ary1[ e ] = v1;
              ary1[ idx ] = v2;
            }
          } else {
            if( v1 > v2 ) {
              ary1[ e ] = v1;
              ary1[ idx ] = v2;
            }
          }
        }
      }
    }
  }
}

具体的なアルゴリズムはググってください。

一部解説を入れておきます。(1)で、比較する要素を求めています(なぜかxorすると比較する要素番号と一致します)。(2)で、不要な比較を飛ばしています。実際の比較はnum(全要素数)の半分で(なぜ半分なのかは英語版wikipediaの図を見ればわかります)、3重目のfor文で飛ばしきれない分を飛ばしています。詳しくは後で説明します。(3)は、昇順か降順かを決めています。マージソートでは順序がどちらか一方向ですが、バイトニックソートは昇順と降順が定期的に入れ替わります。(4)は、交換するか否かの比較です。

先頭にある(5)は、2の冪乗か否かを調べています。マージソート系は、(工夫しなければ)2の冪乗でないとソートできません。2の冪乗はbitが1個しか立たず、1を引くとそのbit以下の全bitが立つ性質があるので、それを利用しています。まじめにfor文を回して2、4、8、16…と調べていくよりはるかに高速です。

2のべき乗でないソートはこちらを参考にしてください。ただ、並列処理はかなり難しそうなので、ここでは使っていません。代わりに、ソートする前に余分な要素すべてにINT_MAXを入れておくことで、2のべき乗でないソートができます。

2. GPGPU最適化

バイトニックソートはそのままGPUに移植できます。しかし、そのまま移植すると効率が悪すぎるので、GPUに適した形に変形します。

ここでは、分岐の除去を行いました。GPUは分岐に弱いことで有名ですが、これは複数のスレッドに対して同時に1命令しか実行できないことが原因です。分岐があると、真になるパスと偽になるパスの両方を実行することになり、結果として実行効率が落ちます(CUDAではWarp divergenceといって、避けろ避けろとうるさく言われます)。元のバイトニックソートはif文が4つもありますが、頑張って1つだけにしました。

// 条件分岐をできるだけ排したバイトニックソート
void bitonicSort2( int num, int *ary1 )
{
  // 2の冪乗でなければ計算しない
  if( ( ( num - 1 ) & num ) != 0 )
    return;

  for( int block = 2; block <= num; block *= 2 ) {
    for( int step = block / 2; step >= 1; step /= 2 ) {
      unsigned int maskLow = (unsigned int)step - 1;
      unsigned int maskHigh = 0xFFFFFFFF ^ maskLow;
      for( int i = 0; i < num / 2; i++ ) { // (A)
        int idx = ( ( i & maskHigh ) << 1 ) | ( i & maskLow );
        int v1 = ary1[ idx ];
        int v2 = ary1[ idx + step ];
        int isAsc = ( ( idx & block ) != 0 ) ? 1 : 0; // (B)
        int isBigger = ( v1 > v2 ) ? 1 : 0; // (C)
        bool needSwap = ( ( isAsc + isBigger ) == 1 ) ? true : false; // (D)
        if( needSwap ) {
          ary1[ idx ] = v2;
          ary1[ idx + step ] = v1;
        }
      }
    }
  }
}

(2)の部分に大きく手を加えました。改良前では、block==2、step==1のとき、(idx, e)はfor文によって(0, 1)、(1, 0)、(2, 3)、(3, 2)、(4, 5)、(5, 4)…となります。このとき、(1, 0)、(3, 2)、(5, 4)など、idxよりeが小さい場合、for文によってスキップされます(もう比較&並び替えは済んでいますから)。これはもったいないので、初めからidx = 0, 2, 4…となるよう工夫します。

step=1のときは1つ置きで飛ばされますが、step==2の場合、idx = 0, 1, 4, 5, 8, 9…というように、2回処理して2回飛ばすの繰り返しになります。このidxの値を条件分岐なして作り出すため、マスクビットとシフトを組み合わせています。要するに、step==2なら2bit目を常に0にし、本来2bit目以降にくる数字を3bit目以降にずらします。step==3なら3bit目以降をずらします。これでidxが作り出せます。その結果不要なループが減るので、(A)のループが半分で済みます。

(B)(C)(D)に3項演算子が出てきます。これも比較といえば比較なのですが、実は条件付き選択(移動)命令をに置き換えることで、条件分岐が発生しないようにすることができます。GeForce GPUのPTXならsetp.cc.s32命令、HLSLならmovc命令、x86 CPUならCMOV命令とかSETcc命令とかANDNPD命令とかに相当します(ちなみに今回のコードをVisualStdio 2012でコンパイルすると、x86ではテーブル参照(たぶん)とかSBB命令とかを使ったややこしいコードに、x64ではSETne命令やSETg命令を使ったコードになりました。CMOV/SETcc/ANDNPDは対応CPUが限られているので、代替命令になるのはしょうがないですね)。もちろん、3項演算子の2項目と3項目で何らかの計算を行っていると、そのパスを実行しなければならないので効率が低下しますが、今回は定数なので問題ありません。

ほかにも共有メモリを使えだとかコアレスアクセスにしろだとか色々ありますが(特にBitonic sortはNVIDIAのReductionサンプルと同じアイデアである程度コアレスアクセス化とカーネル関数呼び出しの削減できるはずです)、面倒なのでやりません。CUDA使いなら、GPU Computing SDKにBitonic sortのサンプルがあるので(Sorting Network)参考にしてください。ちなみに奇偶マージソートや奇数ソート版もあります。DirectXのサンプルにもBitonic sortのサンプル(ComputingShaderSort11)があります。参考までに、OpenCL版もAMDのサイトにあります。

3. ソースコード

長いのでpastebinに貼っておきます。(Githubにも用意する予定)

OpenMPはVisualStudio 2005、PPLはVisualStudio 2010、C++ AMPはVisualStudio 2012が必要です。OpenMPの使用にはコンパイラオプションの変更が必要です。

4. 実行時間

Sortが基本バイトニックソート、Sort2が無駄なループを削減したソート、AmpがSortのGPU移植、Amp2がSort2のGPU移植です。GPUは転送コストを含みます(PCI Express Gen2 x16)。

N=2097152

Core i7 2640M(SandyBridge, 2core/4thread) + Radeon HD 6630M

bench x86 x64
Sort 0.986 0.898
Sort2 1.267 0.836
OpenMP 0.865 0.899
PPL 1.843 1.064
Amp 0.283(FAILED) 0.375(FAILED)
Amp2 0.162 0.195

FAILED : 計算は完了したが結果が不正

Core i7 2600(SandyBridge, 2core/4thread) + GeForce GTX560Ti

bench x86 x64
Sort 0.823 0.868
Sort2 1.099 0.776
OpenMP 0.970 0.799
PPL 0.863 0.517
Amp 0.079 0.098
Amp2 0.047 0.036

.

5. だらだら考察

GPUはほんとに速い。基本バイトニックソートの6倍速いです。なぜかRadeonでのAmpは計算に失敗してしまいましたが、エミュレーション実行(後述)やGeForceでは問題なかったので、GPUかドライバの問題だと思います。

(表にはありませんが)Radeonについては実は速いのはNが十分大きいときだけで、N=65536だとAmp2が約7倍Sortより遅いです。N=131072と262144の境目くらいでCPUとGPUの速度が入れ替わります。これはおそらく関数の呼び出しコストでしょう。関数呼び出しを減らす、Nが小さいときはCPUで処理するなどの工夫が必要です。

なお、GeForceでは小さい数でもCPUより遅くなることはありませんでした(clock()関数での計測なので、マイクロ秒単位で測ったら差が出る可能性はありますが)。GeForceよりRadeonの方が関数呼び出しが遅いと聞いていましたが、嘘ではなさそうです。

前の記事の奇偶転置ソートもそうでしたが、OpenMPの効果が薄く、PPLはx86だけ大幅に遅くなってしまいます。使い方が間違っているような気がしますが、よくわかっていません。

なお、配列の並べ替えではなく、構造体配列の並べ替えをしたいこともあると思いますが、2.のCUDAのソースなどを参考にしてください(作るの面倒…)。

6. C++ AMPのデバッグと、コンパイルされたHLSLコードを見る方法

C++AMPはDirectComputeのラッパーで、DirectComputeはHLSLを使うことを前の記事で説明しました。HLSLコードはGPU上で走るので、普通の方法ではデバッグできません。

VisualStudio 2012でC++AMPコードをデバッグするには、プロジェクトのプロパティを次のように変更します。

構成プロパティ->デバッグ->「デバッガーの種類」を「GPU のみ」に変更

こうするとC++ AMPコードがCPU用にコンパイルされ、エミュレーション実行でデバッグすることができます。ただし、残念ながらCPUコードのデバッグと同時にはできません。

コンパイル後のHLSLコードを確認するには、次のように変更します。

1. 「構成」を「Release」に変更

2. 構成-> C/C++ ->コード生成->「ランタイム ライブラリ」を「マルチスレッド DLL デバッグ(/MDd)」に変更

3. C++AMPコードの先頭にブレークポイントを置いてReleaseモードで実行

4. メニューバーのデバッグ->ウィンドウ->逆アセンブル

このように表示されます

ただし、ReleaseビルドなのにDebugランタイムを使う状態なので、配布用にビルドするときは設定を戻してください(Debugランタイムは配布不可)。

情報元 : C++ AMP Disassembly?

7. PPLで高速ソート

実はMicrosoft公式のPPLのサンプルにBitonic sortがあります。こちらはparallel_for_eachではなくparallel_invokeを使った再帰のコードになっています。実行速度は、N=2097152のとき、i7 2640Mのシングルスレッドで1.052[秒]、マルチスレッドで0.512[秒]、i7 2600のシングルで0.624[秒]、マルチで0.187[秒]でした。…あれ?すごく速い…

…よくわかりませんが私のコードよりこっちを使えばいいんじゃないでしょうか(逃)

広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中