PEZY-SC/SC2を使った話 【KMCアドベントカレンダー 2日目】
この記事はKMC Advent Calendar 2018 - Adventarの2日目の記事です。
概要
理研に設置されているスーパーコンピューター、菖蒲(Shoubu)・菖蒲SystemBでPEZY Computing社のPEZY-SC/SC2を使わせてもらいました。
性能とか使いやすさについて感想とか知見を述べます。
PEZY-SC/SC2について
PEZY Computing社の開発したプロセッサで、メニーコアのMIMD(Multiple Instruction, Multiple Data)プロセッサであることが特徴です。
MIMDなので各コアで全く異なる命令を独立に実行できる点でGPUなどのSIMDプロセッサと異なります。
使えるようになるまで
理研に設置されているスーパーコンピューター、菖蒲(Shoubu)は一般の人でも申請して認められれば無料で使うことができます。 (CUDAまたはOpenCLでの開発経験があったほうが良さそうですが)
http://accc.riken.jp/shoubu_info/application/
私は一昨年12月頃(ちょうどPEZY Computing社の社長(当時)が逮捕された頃です)に、並列オセロ終盤ソルバーの開発の名目で申請しました。
申請して通ったら書類にサインしたのをスキャンして送ったりsshの公開鍵を送ったりといった手続きをすればアカウントが作成されて使えるようになります。
また、申請して使えるようになるのはPEZY-SCを搭載した菖蒲ですが、PEZY-SC2でしか使えない機能(アトミックなど)を使いたい場合はPEZY-SC2を搭載した菖蒲SystemBを使わせてもらえるようです。
使ってみる
基本的な使い方はGPUと同じで、ホストとなるCPU側のコードとデバイス(アクセラレータ)側のコードの両方を書きます。
デバイス側のプログラムはclangに食わせるので普通のC++14が書けます(標準ライブラリは使えませんが)。 この点は古いOpenCLで書かざるを得ない環境と比べるとだいぶ楽だと思います。
メモリ空間はホストとデバイスで別なので明示的に転送しないといけません(CUDAのUnified Memoryのようなことはいまのところできない)。
今回はサンプルとしてMIMDの恩恵が受けられるN-Queen問題を実装してみます。
問題
N*Nの盤面上にチェスのクイーンを互いに効きがないように置く方法は何通りあるか? チェスのクイーンは上下左右と斜めに動ける(飛車と角行を合わせたような動き)。
詳しくは エイト・クイーン - Wikipedia などを見てください。
この問題は各行についてどこに置けるかを試すバックトラックで解けることが知られています。 また、駒の利きをビット列で管理することで、ビット演算を用いて高速に解けることも知られています。 並列化は数行を全探索し、その結果を各スレッドに分配して計算させることでできます。
実装
#include "pzc_builtin.h" #include "../solver.hpp" uint64_t solve(int N, int depth, uint32_t left, uint32_t mid, uint32_t right) { if (depth == N) return 1; uint64_t sum = 0; for (uint32_t pos = (((uint32_t)1 << N) - 1) & ~(left | mid | right); pos; pos &= pos-1) { uint32_t bit = pos & -pos; sum += solve(N, depth+1, (left | bit) << 1, mid | bit, (right | bit) >> 1); } return sum; } void pzc_Solve(const Problem * const probs, uint64_t * const result, size_t count) { size_t offset = get_tid() + get_pid() * get_maxtid(); result[offset] = 0; for (size_t index = offset; index < count; index += get_maxpid() * get_maxtid()) { const Problem &prob = probs[index]; result[offset] += solve(prob.N, prob.depth, prob.left, prob.mid, prob.right); } flush(); }
GPUなどのSIMDアーキテクチャと違って、PEZYプロセッサでは分岐のペナルティーがとても小さいという特徴があるので、このように再帰の中で分岐をするコードを書いても素直かつ高速に動きます。 一方、GPUで高速に動かそうとすると再帰をスタックを用いたループで書くなどコードを変形する必要があり、プログラミングが困難になります。
ホスト側では数行分の全探索をするコードとデバイス側のプログラムを呼び出すコードを書きます。 デバイス側のプログラムを呼び出すのはOpenCLとほぼ同じです。
コード
CPU/GPU版はこちら。
GPU版は再帰をスタックを用いたループで書き直す最適化をしてあります(solve_nonrec関数)。 詳しい最適化の内容については
www.slideshare.net
を見てください(102ページ目から)。
PEZY-SC/SC2版はこちら。
結果
(2020/10/02 追記) N=18、CPUで予め展開する行数は最適なものを選択
マシン | 時間(秒) |
---|---|
Core i7-6700K | 78.857027 |
GeForce GTX 1080 OC | 8.400348 |
Tesla V100 PCIe | 3.936454 |
PEZY-SC | 4.67289 |
PEZY-SC2 | 2.891 |
分岐ペナルティーの小ささなどからPEZY-SC2が最も高速という結果になりました。
さらなる最適化
GPU/PEZY-SC/SC2版ともに、タスク(計算する盤面)によって計算量が異なることにより、後半は遊んでしまうコアが出てきます。 GPUやPEZY-SC2ではアトミック命令が使えるので、アトミック命令を使うなどして動的にタスクを割り振ることでさらに高速化することができます。 アトミック命令の使い方などはNDAを結ばないと教えてもらえませんが、ソースコードは公開OKなので、私の公開しているオセロソルバーのソースコード
などから類推することでおおよそはわかるのではないでしょうか。 実際、PEZY-SC2上でアトミック命令を用いて最適化したところ、2.40199秒まで高速化されました。
GitHub - primenumber/PEZY-nqueen at sc2
オセロソルバーの実装
オセロソルバーとはあるオセロの局面が与えられたときにそこから両者が最善を尽くしたときの試合結果(何石差か)を求めるプログラムのことです。 基本的にはαβ法という再帰的なアルゴリズムを用いて解きます。
すでにCPU/GPU版を実装していたので、簡単にPEZY-SC/SC2版も実装できると思っていました。
CPU版: github.com
しかし、N-Queen問題と同じようにオセロソルバーも実装しようとしたのですが、空きマスが多くなると(=再帰の深さが深くなると)エラーが出て正しく動きませんでした。 原因を探ってみると、どうやらスタックオーバーフローしているようです。 実は、PEZY-SCにはローカルメモリが16KBあり、これを1コア内の8スレッドで分割してスタック領域として使います(ちなみに、スタック領域を縮小して余ったローカルメモリをコア内の共有メモリとして使うこともできます)。 そのため、スタック領域は各スレッドにつき2KBしかありません。PEZY-SC2でローカルメモリが20KBになりましたが、それでもスレッドあたりは2.5KBです。 N-Queenより複雑なオセロソルバーではスタック領域が不足してしまっていたのです。 (再帰で書いているのも一因ですが、スタック消費量を気にせずコードを吐くLLVM側も悪い気はします)
そこで、泣く泣くスタック領域を大量に消費する再帰で書くことを諦めて、自分で実装したスタックをグローバルメモリ上に置いてループで解くことにしました。 これで解けるようにはなりましたが、プログラミングが非常に難解になってしまいました。 このあたりのプログラミングの難しさが解消されるともっと良いと思うのですが…(スタック領域は基本グローバルメモリに置かれて、キャッシュで高速化するなど、難しそうではありますが)
何とかプログラミングの困難さを乗り越えて実装した結果が
です。 根の1段だけFastest-First Heuristic(速さ優先探索)を行っています(複数段に適用しようとするとスタックが足りなかった)。 また、PEZY-SC2向けに最適化を行ったものが
GitHub - primenumber/PEZY-Othello at sc2
にあります。 最適化の内容は
- 速さ優先探索を複数段に適用
- ビット演算命令を使用
- アトミック命令を使用して動的なタスク割り振り
- 置換表の利用
- 葉の近くで速さ優先探索を適用していないところで隅を先に読む静的move ordering
- NegaScout法の適用
になります(まだパラメータ調整をちゃんとしていなくて速度が出ない場面もありますが…)。 これらの改良により、20石空きの局面を219万局読ませて、1ノード(8モジュール)で6時間弱で解けるようになりました(Tesla V100でも測るべきですが、クラウドで8GPUとかぶん回すと結構高いのでやってません、あとGPU版はいくつかの改良を取り入れてないです)。 菖蒲SystemBでは最大1週間のジョブを投げられるので、22~23石空きまでは1ノードで読ませることができそうです(1石空きマスが増えるごとにだいたい3倍ぐらいかかる時間が増える)(ちなみにノード数を増やしてもゲーム木の大きい局面が律速になるのでそんなに高速化しない)。
今後は統計評価関数を利用したmove orderingを実装して更なる高速化を果たして、24石空き程度まで読ませたいと考えています。
感想
PEZY-SC/SC2はメニーコアのMIMDプロセッサという他にはあまりないアクセラレータです。 今回は再帰の中で分岐するようなMIMDに有利なワークロードではPEZY-SC2がTesla V100を超えるパフォーマンスを出せることを示せました。次世代のPEZY-SC3も近いうちに登場するみたいなので楽しみですね。
一方で、スタック領域の制約にはかなり苦しめられました。せっかく再帰的なタスクに有利なので、たくさんスタック領域が使えるとプログラミングがしやすくてよいと思うので、この辺が改善されてほしいです。
最後に
KMCアドベントカレンダー 3日目の記事は id:nojima718 (KMC-ID: nojima) さんの「ぷよぷよのプレイ動画を解析して棋譜を生成する」です。お楽しみに!
ビット演算マニアのためのAVX512入門 【KMCアドベントカレンダー 22日目】
この記事はKMCアドベントカレンダー22日目の記事です。
この記事では全国一億三千万のビット演算マニアのために、AVX512命令セットからビット演算に使えそうなものを紹介します。
AVX512の基本
AVX512とは
Advanced Vector Extensions(AVX)というx86/x64のベクトル演算拡張命令を512ビット幅に拡張したものです。 普通のCPUではSkylake-SP/X/Wから使えるようになりました。 512bitのデータを一括で処理できます。ビット演算なら512並列ですよ512並列。
AVX512になって変わった主なこと
AVX512の命令セット
AVX512は複数の要素から構成されており、プロセッサによってサポートの状況が異なります。
以下に書くサブセットとその説明を書きます。
- AVX512-F (Foundation) 基本命令セット。AVX512をサポートするプロセッサは必ずこれをサポートする。
- AVX512-CD (Conflict Detection) データの衝突検出。Hist処理とかがベクトル化できる。 vpconflictd - Qiita
- AVX512-DQ (Doubleword and Quadword) 32bit/64bit単位での処理。
- AVX512-BW (Byte and Word) 8bit/16bit単位での処理。
- AVX512-VL (Vector Length) 128bit(XMM)/256bit(YMM)レジスタを操作する。
最新のCPU(Skylake-SP/X/W)でもまだ使えない拡張もたくさんあります。
- AVX512-PF (Pre Fetch) プリフェッチ命令。
- AVX512-ER (Exponential and Reciprocal) 指数関数や逆数関数。
- AVX512-IFMA52 (Integer Fused Multiply Add) 整数の積和算を倍精度の積和算演算器を借りて52ビット精度で実行する。
- AVX512-VBMI (Vector Byte Manipulation Instructions) BWに入っていないbyte操作の命令。
- AVX512-VNNI (Vector Neural Network Instructions) ニューラルネットワーク処理用の整数演算命令。16bit整数同士の積を32bit整数と足す。
- AVX512-4VNNIW (Vector Neural Network Instructions Word variable precision) ニューラルネットワーク処理用の整数演算命令。16bit整数同士の積を32bit整数と足すのを4回繰り返す。
- AVX512-4FMAPS (Fused Multiply Accumulation Packed Single precision) 単精度の内積を4回繰り返す。これもニューラルネットワーク処理用。
- AVX512-VPOPCNTDQ (Vector POP CouNT) ベクトル化されたpopcount。
- AVX512-VBMI2 (Vector Byte Manipulation Instructions2) BWにもVBMIにも入っていないbyte操作の命令。
- AVX512-BITARG (不明) ビット操作関連の命令。
マスク演算
AVX512では比較演算の結果などはマスクレジスタと呼ばれる特殊なレジスタに格納されます。
AVX512では命令にマスクとしてマスクレジスタを指定することができます。
VPANDQ (and演算, 64bit単位のマスク)を例にして疑似コードで表すと、
for j in range(0, 8): i = j*64 if mask[j]: dest[i+63:i] = src[i+63:i] else: dest[i+63:i] = a[i+63:i] AND b[i+63:i]
このようにマスクでビットが立っている場所を演算せずに値を残すことが可能です。また、値を残す代わりに0にすることもできます。
intrinsicでは
dest = _mm512_and_epi64(a, b); // マスク無し dest = _mm512_mask_and_epi64(src, mask, a, b); // マスクあり(マスクされたらsrcを残す) dest = _mm512_maskz_and_epi64(mask, a, b); // マスクあり(マスクされたら0)
のように使います。
マスク同士の演算
このマスクレジスタはSIMDレジスタ(XMM/YMM/ZMMレジスタ)とも通常のレジスタとも違うので、通常のレジスタとマスクレジスタとの間で値をやりとりするのにもコストがかかります。 そこで、いくつかの演算(and, or, addなど)はマスクレジスタ同士で行えるようになっています。
命令の解説
AVX512-F
VPTERNLOG{D,Q}
ビット単位の任意の3項演算ができます。演算内容は即値で与えます。
A & (B | ~C)
なら
A | B | C | A&(B|~C) |
---|---|---|---|
1 | 1 | 1 | 1 |
1 | 1 | 0 | 1 |
1 | 0 | 1 | 0 |
1 | 0 | 0 | 1 |
0 | 1 | 1 | 0 |
0 | 1 | 0 | 0 |
0 | 0 | 1 | 0 |
0 | 0 | 0 | 0 |
なので即値として0b11010000 = 0xD0を与えてやれば良いです。
intrinsicは_mm512_ternarylogic_epi64(A, B, C, 0xD0)
のように使います。
これを使えば他のビット論理演算命令いらないじゃんとなりそうですが、怪しいページによるとスループットが2命令/サイクルとandやorなどの3命令/サイクルに比べると劣るので適材適所で使うことになりそうです。
VP{AND,OR,XOR,ANDNOT}{D,Q}
and, orなどの演算です。マスクが32bit/64bit単位で指定できます。
VPS{L,R}{L,A}{,V}{W,D,Q}
16bit/32bit/64bit単位のビット{論理,算術}シフトです。V付きはシフト量もベクトルで指定することで、各データごとに異なるシフト量を指定できます。
VPRO{L,R}{,V}{D,Q}
待望のローテート命令です。でも32bit/64bitだけ。
AVX512-CD
VPLZCNT{D,Q}
leading zero countです。上位ビットから見て最初に1が立っているビットの位置を調べる。整数→浮動小数の変換で前からほぼ同じことができましたが、キャストとか0の特別扱いとか諸々せずに済むようになりました。 ところでなんでConflict Detectionに入っているんでしょうか。
AVX512-VBMI
VPERMI2B
各バイトの下位7bitをインデックスにして、2つのZMMレジスタをくっつけた128bytesの中から表引きします。 これまでの128bit境界に捕われていたり、即値で指定しないといけなかったりしたシャッフル命令とはおさらばできます。 とはいえスループットやレイテンシが大きかったら従来のシャッフル命令も併用しないといけなさそうですね。
VPMULTISHIFTQB
各64bit要素に対して各バイトの分だけ右シフトしたデータを返します。
for i in range(0, 8): q = i*64 for j in range(0, 8): tmp8 = 0 ctrl = src1[q+j*8+7:q+j*8] & 0x3F for l in range(0, 8): tmp8[l] = src2[q+((ctrl+l)&0x3F)] dst[q+j*8+7:q*j*8] = tmp[7:0]
AVX512-VBMI2
VPSH{L,R}D{,V}{W,D,Q}
2つの値をつなげてシフトします。多倍長シフトとかに使えますね。V付きの命令はシフト量もベクトルで指定することで、各データごとに異なるシフト量を指定できます。
AVX512-POPCNT
VPOPCNT{D,Q}
32bit/64bit単位でのpopulation count(1になっているビットの数を数える)です。
AVX512-BITARG
VPOPCNT{B,W}
8bit/16bit単位でのpopcountです。なぜかVPOPCNT{D,Q}とは別の拡張になっています。
VPSHUFBITQMB
謎命令。各64bit要素に対して、各バイトをインデックスとして64bitデータの指定したビットを取ってきてマスクにします。疑似コードを見るのが早いと思います。
for i in range(0, 8): # qword for j in range(0, 8): # byte m = src2[i*64+j*8+7:i*64+j*8] & 0x3F res[i*8+j] = src1[i*64+m] return res
AVX512以外の拡張
VPCLMULQDQ
繰り上がりなしの64bit乗算を並列で行います。ビット演算的にはCRCの計算とかpclmulqdqを用いたビット単位のunpack【ビット演算テクニック Advent Calendar 2016 11日目】 - prime's diaryとかに使えそうですね。
便利Intrinsic
_mm512_reduce_{or,and}_epi{32,64}
32bit/64bitデータ16/8個のOR/ANDを取ります。
_mm512_reduce_and_epi64
なら、
reduced[63:0] = 0xFFFFFFFFFFFFFFFF for j in range(0, 8): i = j*64 reduced[63:0] = reduced[63:0] AND a[i+63:i] return reduced[63:0]
まとめ
AVX512にはどう使うかよくわからない便利そうな命令がたくさん追加されています。
使えるCPUを手に入れた暁にはうまく使ってビット演算ライフを満喫しましょう!
明日のKMCアドベントカレンダーはwass88さんによる「役に立たないこと」の予定です。お楽しみに!
ICPC 2017 Tsukuba Asia Regional 参加記
2017/12/16-18に行われたICPC 2017 Tsukuba Asia RegionalにチームPrimeDragonとしてamano, nikuttoと参加してきました。
コンテスト前日まで
15日
生活リズム的に16日の朝早くから出られそうにないので徹夜することにした。Redbull飲んで徹夜する。
16日
徹夜のおかげで無事新幹線に乗れる。
初つくばエクスプレス乗車。(京都大阪間に比べると)値段高いですね…
中学校の知人とばったり出会う。偶然ではなくてICPCに出てるらしい。
JAXA見学。管制システムかっこいい…
リハーサルではコンテスト時の段取り(誰が最初にテンプレートとか書くのかとか)を確認、コンテストシステムに慣れる。 プログラムが実行されるサーバーと手元の速度差を確認。(ほとんど違いはなかった) 印刷やClarを投げる練習もした。
筑波学院大学でWelcome Party。 京大の伝統(?)でチーム紹介は短めに。 料理がだいたい茶色系だったのはなんでだろう…おいしかったけど
歩いて会場から宿まで歩いた。さむい。
ちょうどARCがあったので参加。3完できたのでまあまあ良い結果。
寝る。
コンテスト当日(17日)
朝
緊張してあんまり寝られなかったけどおかげで寝坊せず起きられた。
朝食はつくばカピオの近くのCASAってところで食べた。natsugiriさんとかがいた。
コンテスト本番
記憶が曖昧なところがあるので時系列おかしいかも知れないです。
コンテスト開始。
まず、Xmodmapの設定(Caps LockをCtrlにする)、.vimrcやテンプレートを書く。 その間にAやBを読んでもらう。
書けた時点でAが読めて解法もわかったらしいので教えてもらって書く。AC。 Bが解けるらしいので交代する。
Cを読む。特に難しいところはなさそう。Bを解いている間に実装を詰める。 BがACしたので交代してCの実装をする。サンプルが通らなくて焦るが、インデックスがずれているだけだった。 2個もずれててこんなずれるはずないだろと思っていたが、サンプルは通るので投げる。AC。
この時点で10位ぐらいだった気がする。出だしとしてはいい感じ。
ここからは難易度順に並んでいないので簡単そうな問題をさがす。 Fが始点と終点からDijkstraすれば良さそうなので解けると主張するが、それだけだと経路が長くなるか判定できないと言われてたしかにとなる もうちょっと考えると最短経路上の辺で作った部分グラフ(DAGになる)で橋になっているかがわかれば良いことがわかる。 コードを書いてサンプルも通ったので投げる。WA。コード印刷するが原因わからず。
その間にGで実装軽そうな方針が立ったので書いてもらう。AC。
まだわからないので他の問題を考えていた。Iが解けるらしいので書いてもらう。AC。
ようやくFで橋の判定アルゴリズムが間違っていたことに気づき(DAGという条件で軽く書けると思っていたが、嘘解法だった)、ライブラリの橋検出を写す。AC。
この時点で5位くらいだったかな?かなり良い順位。このまま比較的解いている人の多いEを通せば良さそう。
しかし、EでDP解(これが正解)ではなく貪欲法解をずっと考察していていくつも嘘解法を作ってはだめとなっていた。
結局Eは通せずにコンテスト終了。
コンテスト後
解説&順位発表&表彰。D実装つらそう…こんなんコンテスト中に通せる気がしない… 結果は6完で7位。
夕食&企業ブース見学。料理うまい。前日に比べて彩りもあってよかった。 企業ブースを一通り回る。 レコチョクのパズルを解いて景品(モバブ、アメなど)をもらう。このことが後々生死を分かつことになる。 あとドワンゴのWA投げでトランプをもらった。
企業賞みたいなのでレコチョクからワイヤレスイヤホンをもらう。結構良いやつっぽい。
E通してればWF行けたんでは?とか一応まだWFの可能性ありそう、とかのコメントももらった。まさかWFの可能性がある順位取れると思ってなかったので嬉しいが、余計Eを解けなかったことが辛くなる。
昨日と同じく歩いて宿に帰る。つくばカピオからは近い。
この日はなかなか眠れなかったのでオセロでモンテカルロ木探索を実装した。合法手生成でバグってうまく動かなかった。
ホテルの自販機でDrPepperが売っていたので人生で初めて飲む。謎の味がした。 よく考えたら寝れなくて困ってるのにカフェイン入り飲料を飲むのは間違っていた。
結局朝6時頃に就寝。
コンテスト翌日(18日)
企業見学のためつくばから東京に移動しようとするが、つくばエクスプレスが30分ぐらい遅れて遅刻する。
企業見学前半はLINE。以前ISUCONで行ったことがあるが、当時は渋谷にあったので新しい新宿のオフィスに来るのは初めて。 窓からの眺めが良かった。あと叙々苑の焼肉弁当が美味しかった。
企業見学後半はPFN。ICPC/JOI出身者がすごく多いらしい。ロボットの制御もやっていたのは知らなかった。
解散。このあとせっかく東京に着たのでオフ会して終電で帰った。途中で持ってきていたモバブの充電が切れたが、レコチョクからもらったモバブのおかげで生きながらえた。
感想
E通したかった…
オセロの着手可能位置生成 【ビット演算テクニック Advent Calendar 2016 23日目】
はじめに
この記事はビット演算テクニック Advent Calendar 2016の23日目の記事です。
オセロとビット演算
オセロは盤面サイズが8x8=64マスなので、64bit変数を2つ用意し、1つめを自分の石(あるいは黒石)の位置を、2つめを相手の石(あるいは白石)を表すために使うと合計128ビットで自然に表現できます。*1
自然に表現できるだけでなく、盤面の回転や反転などもDelta swapとその応用 【ビット演算テクニック Advent Calendar 2016 3日目】 - prime's diaryで紹介したような方法で効率的に行なえます。
着手可能位置の生成はタテ・ヨコ・ナナメの一直線上をpextなどで取ってきて表引きをすることにより、ある程度の速度で可能ですが、(38=6561なのでL1キャッシュに乗る程度の大きさとはいえ)メモリアクセスが必要になってしまいそれなりの時間がかかります。
以下、ビットボードには以下のようにマス目とビット列の対応が付いているとします(0が最下位ビット)。
A | B | C | D | E | F | G | H | |
---|---|---|---|---|---|---|---|---|
1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
2 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
3 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 |
4 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 |
5 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 |
6 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 |
7 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 |
8 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 |
ひっくり返る石の計算
参考文献の方法を使っています。 参考文献のOpenCL風のコードがわかりやすいのでそれを元に動作を解説します。
ulong flip(int pos, ulong P, ulong O) { ulong4 flipped, OM, outflank, mask; OM.x = O; OM.yzw = O & 0x7e7e7e7e7e7e7e7eUL; // 1 mask = (ulong4) (0x0080808080808080UL, 0x7f00000000000000UL, 0x0102040810204000UL, 0x0040201008040201UL) >> (63 - pos); // 2 outflank = (0x8000000000000000UL >> clz(~OM & mask)) & P; // 3 flipped = (-outflank * 2) & mask; // 4 mask = (ulong4) (0x0101010101010100UL, 0x00000000000000feUL, 0x0002040810204080UL, 0x8040201008040200UL) << pos; outflank = mask & ((OM | ~mask) + 1) & P; // 5 flipped |= (outflank - (outflank != 0)) & mask; return flipped.x | flipped.y | flipped.z | flipped.w; }
- 盤面の端を超えて石が繋がっていても無視する
- posから見て上、左、右上、左上のマス目をマスクするビット列を計算する
- posからみて相手の石が途切れたところをclzで求め、自分の石の位置とandを取ることでひっくり返せるかどうかがわかる。outflankの各要素は高々1ビットが立つ。
- -outflank * 2でoutflankより上位のビットが全部1になる。これとマスクを取ることでひっくり返る石の位置がわかる。
- 後半部分は下位ビットから見ていって相手の石が途切れたところを探します。これは関係ない場所を1埋めしてから+1すればよいです。あとは単に前半部分の逆です。
実際のAVX2をゴリゴリ使った実装はissen/move_generator.cpp at 72f450256878094ffe90b75f8674599e6869c238 · primenumber/issen · GitHubにあります。
このコードでは各方向に対するひっくり返る石をビット演算で求めることができ、AVX2を使うと256bit一気に操作でき高速です。ただ、clzの並列版はないので、立っている一番上のビット位置を頑張って求めます。
Bit scan reverseをAVXで並列化する - prime's diaryと違い、1<<bsr
を求める関係で浮動小数点数への変換より下位ビット埋め+vpshufb+最下位ビット抽出のほうが速いです。
立っている一番上のビットを求める実装例: issen/move_generator.cpp at 72f450256878094ffe90b75f8674599e6869c238 · primenumber/issen · GitHub
着手可能位置の計算
着手可能位置(mobility)の計算をビットボードで行う方法は
で言及されているビットシフトを繰り返す方法が知られているようです。今回はそれより高速な方法を解説します。
この方法では、石を置くと左に向かって相手の石をひっくり返せるような位置を一気に列挙します。
これを盤面を回転・反転させて8方向に対して行うことで着手可能位置を計算することができます。
1列分の着手可能位置を求めるコードは以下の通り。
// p: 自分の石, o: 相手の石 uint8_t mobility_line(uint8_t p, uint8_t o) { uint8_t p1 = p << 1; return ~(p1 | o) & (p1 + o); }
このプログラムの動作を具体例を元に解説します。
例えば-oop-opo
(p: 自分の石, o: 相手の石, -: 空きマス)のとき、
7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 | |
---|---|---|---|---|---|---|---|---|
p | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 |
o | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 1 |
pを1ビット左シフトして、
7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 | |
---|---|---|---|---|---|---|---|---|
p1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 |
o | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 1 |
2つを足し算すると、自分の石より上位にある相手の石が途切れる位置(つまり着手可能位置) または 孤立した相手の石 または 孤立した自分の石の一つ上位 のビットが立ちます。
7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 | |
---|---|---|---|---|---|---|---|---|
p1+o | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 1 |
~(p1|o)でマスクを取って、
7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 | |
---|---|---|---|---|---|---|---|---|
res | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
となり着手可能位置が計算できました。 これを64ビット*2の盤面全体に対して一気にやろうとすると、8ビット単位のシフトや足し算が必要になります。
うまくマスクを取って64ビット演算に帰着してもいいですが、SSE/AVX2には8ビット単位の加算が用意されているので、それを使うと良いでしょう。1ビットの左シフトはp+pで実現できます。 SSEなら2並列、AVX2なら4並列で着手可能位置を求められるので、8方向に対する計算をかなり高速化することができます。
全体をSSE/AVX2で書いたものがこちらになります。
issen/movable_generator.cpp at master · primenumber/issen · GitHub
参考文献
ひっくり返る石の計算方法 Reversi - bitboard and GPU computing
すでに知られている着手可能位置の計算方法 リバーシプログラム - bitboard による合法手生成 - プログラミング練習帳&備忘録
*1:1つ目で両方の石の位置を、2つ目で相手の石の位置を表す流儀もある
GPGPUで爆速でオセロを解く 【KMC Advent Calendar 2016 19日目】
はじめに
この記事はKMC Advent Calendar 2016の19日目の記事です。
GPGPUとは
GPGPU(General-purpose computing on graphics processing units; GPUによる汎用計算)とは、GPUの演算資源を画像処理以外の目的に応用する技術のことである。
(出典: GPGPU - Wikipedia )
GPUは並列演算に特化しており、大量のスレッドをハードウェアで処理でき、汎用CPUに比べて演算処理性能が高いという特徴があるので、これを用いて演算を高速に行おうというのがGPGPUのアイデアです。
ただ、良いことばかりではなく、汎用CPUに比べてプログラミング上の制約が大きいなど、性能を出し切るのはかなり難しいです。
GPGPUでオセロを解く
今回はGPGPUを用いてオセロの終盤の完全読み*1を高速化します。
オセロを解く際に基本となるのがアルファベータ法です。
詳しい解説は他に譲るとして、このアルゴリズムの肝はすでに計算した部分の結果を利用して枝刈りを行うことです。そのため、並列化が難しいアルゴリズムになっています。 今回はアルゴリズムレベルでの並列化は後回しにして、たくさんの局面を並列で完全読みすることを考えます。
各スレッドに一つの局面を割り当てて完全読みを行うのが素直な実装ですね。
アルファベータ法は再帰を用いると自然に書けるので再帰で実装してみます。 今回はCUDAで実装します(したがって、この記事では以後GPUはNVIDIAのGPUと仮定します)。CUDAで再帰するときは必要なスタック数が静的に見積もれないので必要に応じて動的にスタックサイズを指定してやります(今回は4096バイト/スレッドにしました)。*2 高速化のためビット演算を用いてひっくり返る石の位置を計算しています。詳しくはReversi - bitboard and GPU computingを参照してください。
CUDAカーネル部分だけ示します(CPU側のコードはスタックサイズの指定を除いて後で示す改良バージョンと同じです)。
constexpr int threadsPerBlock = 128; using ull = unsigned long long; struct Board { ull player; ull opponent; __host__ __device__ Board() = default; __host__ __device__ Board(ull p, ull o) : player(p), opponent(o) {} __host__ __device__ Board(const Board &) = default; __host__ __device__ Board(Board &&) = default; __host__ __device__ Board &operator=(const Board &) = default; __host__ __device__ Board &operator=(Board &&) = default; }; __shared__ int count[threadsPerBlock]; __device__ char score(const Board &bd) { char pnum = __popcll(bd.player); char onum = __popcll(bd.opponent); if (pnum == onum) return 0; if (pnum > onum) return 64 - 2*onum; else return 2*pnum - 64; } __constant__ ull mask1[4] = { 0x0080808080808080ULL, 0x7f00000000000000ULL, 0x0102040810204000ULL, 0x0040201008040201ULL }; __constant__ ull mask2[4] = { 0x0101010101010100ULL, 0x00000000000000feULL, 0x0002040810204080ULL, 0x8040201008040200ULL }; __device__ ull flip(const Board &bd, int pos, int index) { ull om = bd.opponent; if (index) om &= 0x7E7E7E7E7E7E7E7EULL; ull mask = mask1[index] >> (63 - pos); ull outflank = (0x8000000000000000ULL >> __clzll(~om & mask)) & bd.player; ull flipped = (-outflank * 2) & mask; mask = mask2[index] << pos; outflank = mask & ((om | ~mask) + 1) & bd.player; flipped |= (outflank - (outflank != 0)) & mask; return flipped; } __device__ ull flip_all(const Board &bd, int pos) { return flip(bd, pos, 0) | flip(bd, pos, 1) | flip(bd, pos, 2) | flip(bd, pos, 3); } __device__ Board move_pass(const Board &bd) { return Board(bd.opponent, bd.player); } __device__ char alpha_beta(const Board &bd, char alpha, char beta, bool passed) { ++count[threadIdx.x]; ull puttable = ~(bd.player | bd.opponent); if (puttable == 0) return score(bd); bool pass = true; for (; puttable;) { ull bit = puttable & -puttable; puttable ^= bit; int pos = __popcll(bit-1); ull flipped = flip_all(bd, pos); if (flipped) { pass = false; Board next(bd.opponent ^ flipped, (bd.player ^ flipped) | bit); alpha = max(alpha, -alpha_beta(next, -beta, -alpha, false)); if (alpha >= beta) return alpha; } } if (pass) { if (passed) { return score(bd); } else { return -alpha_beta(move_pass(bd), -beta, -alpha, true); } } return alpha; } __global__ void search_noordering(const Board *bd_ary, int *res_ary, int *nodes_count, const size_t size) { size_t index = threadIdx.x + blockIdx.x * blockDim.x; while(index < size) { count[threadIdx.x] = 0; res_ary[index] = alpha_beta(bd_ary[index], -64, 64, false); nodes_count[index] = count[threadIdx.x]; index += blockDim.x * gridDim.x; } }
GTX1080(ZOTAC Geforce GTX 1080 AMP Edition, Base 1683MHz/Boost 1822MHz)でGGSの棋譜から生成した10石空きの問題1091780局面の完全読みを実行してみると5分48秒かかりました。
一方、CPUで各種高速化を施したソルバー*3ではCore i7-6700K上で同じ処理を行うのに1分13秒しかかかりませんでした。これではGPGPUにした意味がほとんどありません。
高速化する
2018/03/31追記
このプログラムが遅い最も大きな要因は、以下で説明するレイテンシの隠蔽失敗によるものではなく、再帰的な分岐により演算器が遊んでしまうことによるものです。 実際、Shared Memoryを少し節約しつつ、1スレッドで1局面読むようにすると下記の結果の9秒より速く計算することが出来ます。
追記終わり
GPUはSM(Streaming Multiprocessor)という単位でワープを管理していて、GTX 1080(や他のCompute Capability 6.1のGPU)では同時に4つのワープを実行できます。しかし、実際にはSMは4つよりも多くのワープを管理でき、実行可能なワープを選んで実行することで、命令やメモリの読み書きのレイテンシを隠蔽しています。逆に言うと、同時にたくさんのワープを管理できなければレイテンシを隠蔽できず、多くの時間演算器は遊んでしまうことになります。
今回の場合、再帰でレジスタをかなりたくさん消費してしまっているので、一つのSMで持てるワープは数個しかありません。
そこで、Shared Memoryに再帰のスタックの内容を移し、再帰をループに直すことにします。 こうした場合今度はShared Memoryの使用量がSMの持てるワープ数を制限してしまうので、4スレッドで1局面を読むことで使用量を減らし、その代わり4スレッドをSIMD的に使うことでひっくり返る石の場所の計算を高速化します。
また、一つの局面を読み終わったらすぐ次の局面を読めるように変更しました。これにより、ワープの中で最も遅いスレッドが律速になるのを軽減できます。
その他の細かい改良を含めた結果、実行時間は最終的に9秒になりました。CPUの8倍高速です。
ソースコード及びベンチデータ
ソースコードは以下の通り。 問題ファイルの読み込みはBase81という1文字で4マス分表すフォーマットです。詳しくはGitHub - primenumber/issen: オセロAIの高速な実装を参照。
ベンチマークに使用した入力ファイルは次のリンクからダウンロードできます
今後の課題
枝刈りの成功率を上げるため、次の手を読む順番を良さそうな順番になるように並べ替える。 並べ替えのためにメモリをかなり消費するので、グローバルメモリ上に置かないと実行効率が上がらないと予想される。
多数の局面を並列に読むのではなく、アルゴリズムを工夫して一つの局面を並列に読めるようにする。 あわよくばCPUで読んだときより速くならないだろうか…
参考文献
Reversi - bitboard and GPU computing 五石空きをminimax法で完全読みするのをGPUにやらせてみたがあまり速くなかったという結果。 ビット演算部分はここを参考にした。
GPU Acceleration for Board Games | NVIDIA Developer この記事と似ていて、並列にオセロのAIを走らせているらしい。Windowsを持ってないので動作を確かめられない。
*1:その局面から両者が最善を尽くしたときに、勝敗及び石差がどうなるかを判定すること。
*2:http://primenumber.hatenadiary.jp/entry/2016/11/27/031429
Bit scan reverseをAVXで並列化する
Bit scan reverse (bsr)
00000001010100010100000100101101 <-clz->*<---bit scan reverse--->
1になっている最も高位のビットのインデックスを0-indexedで計算します。0に対する結果は未定義(不定)です。 似た処理にcount leading zero(clz)があり、こちらは最上位ビット側から数えます。(0に対する結果は整数型のビット数num_bitsになる) したがって、0以外の値に対してはbsr + 1 + clz = num_bitsが成り立ちます。
bsr/lzcnt命令
x86にはズバリbsr命令があり、32ビット/64ビットの値に対してBit scan reverseを計算できます。 また、BMI拡張命令セットにはlzcntという命令があり、clzも1命令で計算できます。
しかし、bsr、lzcntともにSSEやAVX拡張命令セットの中に対応する命令がないため、並列に計算することができません。
今回はAVX命令を使って256ビット幅で並列にBit scan reverseを計算する方法を考えます。
また、8/16/32/64ビットそれぞれについて、どのアルゴリズムが最適なのかベンチマークを取ります。
アルゴリズム
スカラーに変換してbsrで処理 (naive)
まずは単純にスカラー値に変換してbsr命令で処理するのが一番自然でしょう。 pextr{b,w,d,q}は意外にコストが高いのでvmovdqaで一括でメモリ上(実際にはL1キャッシュに乗るはず)にコピーして処理することにします。 計算結果をYMMレジスタに格納するのは_mm256_set_epi* intrinsicsに任せました(メモリ上に書き戻してvmovdqaで格納するより速かった)
inline __m256i bsr_256_8_naive(__m256i x) { alignas(32) uint8_t b[32]; _mm256_store_si256((__m256i*)b, x); return _mm256_setr_epi8( __bsrd(b[ 0]), __bsrd(b[ 1]), __bsrd(b[ 2]), __bsrd(b[ 3]), __bsrd(b[ 4]), __bsrd(b[ 5]), __bsrd(b[ 6]), __bsrd(b[ 7]), __bsrd(b[ 8]), __bsrd(b[ 9]), __bsrd(b[10]), __bsrd(b[11]), __bsrd(b[12]), __bsrd(b[13]), __bsrd(b[14]), __bsrd(b[15]), __bsrd(b[16]), __bsrd(b[17]), __bsrd(b[18]), __bsrd(b[19]), __bsrd(b[20]), __bsrd(b[21]), __bsrd(b[22]), __bsrd(b[23]), __bsrd(b[24]), __bsrd(b[25]), __bsrd(b[26]), __bsrd(b[27]), __bsrd(b[28]), __bsrd(b[29]), __bsrd(b[30]), __bsrd(b[31])); } inline __m256i bsr_256_16_naive(__m256i x) { alignas(32) uint16_t b[16]; _mm256_store_si256((__m256i*)b, x); return _mm256_setr_epi16( __bsrd(b[ 0]), __bsrd(b[ 1]), __bsrd(b[ 2]), __bsrd(b[ 3]), __bsrd(b[ 4]), __bsrd(b[ 5]), __bsrd(b[ 6]), __bsrd(b[ 7]), __bsrd(b[ 8]), __bsrd(b[ 9]), __bsrd(b[10]), __bsrd(b[11]), __bsrd(b[12]), __bsrd(b[13]), __bsrd(b[14]), __bsrd(b[15])); } inline __m256i bsr_256_32_naive(__m256i x) { alignas(32) uint32_t b[8]; _mm256_store_si256((__m256i*)b, x); return _mm256_setr_epi32( __bsrd(b[0]), __bsrd(b[1]), __bsrd(b[2]), __bsrd(b[3]), __bsrd(b[4]), __bsrd(b[5]), __bsrd(b[6]), __bsrd(b[7])); } inline __m256i bsr_256_64_naive(__m256i x) { alignas(32) uint64_t b[4]; _mm256_store_si256((__m256i*)b, x); return _mm256_setr_epi64x( __bsrq(b[0]), __bsrq(b[1]), __bsrq(b[2]), __bsrq(b[3])); }
popcountに帰着 (popcount)
uint32_t bsr32(uint32_t x) { x |= x >> 1; x |= x >> 2; x |= x >> 4; x |= x >> 8; x |= x >> 16; return popcount32(x>>1); }
右シフトとORを組み合わせることで1の立っている一番高位のビットより下位のビットを1で埋めることができます。
これをpopcountと組み合わせればbsrが求められます。
inline __m256i popcount_256_8(__m256i x) { x = _mm256_sub_epi8(x, _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi8(0x55))); x = _mm256_add_epi8(_mm256_and_si256(x, _mm256_set1_epi8(0x33)), _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi8(0x33))); return _mm256_and_si256(_mm256_add_epi8(x, _mm256_srli_epi16(x, 4)), _mm256_set1_epi8(0x0F)); } inline __m256i bsr_256_8_popcnt(__m256i x) { x = _mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi8(0x7F))); x = _mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi8(0x3F))); x = _mm256_and_si256(_mm256_srli_epi16(_mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi8(0x0F))), 1), _mm256_set1_epi8(0x7F)); return popcount_256_8(x); } inline __m256i popcount_256_16(__m256i x) { x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x5555)), _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi16(0x5555))); x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x3333)), _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi16(0x3333))); x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x0F0F)), _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi16(0x0F0F))); return _mm256_maddubs_epi16(x, _mm256_set1_epi16(0x0101)); } inline __m256i bsr_256_16_popcnt(__m256i x) { x = _mm256_or_si256(x, _mm256_srli_epi16(x, 1)); x = _mm256_or_si256(x, _mm256_srli_epi16(x, 2)); x = _mm256_or_si256(x, _mm256_srli_epi16(x, 4)); x = _mm256_srli_epi16(_mm256_or_si256(x, _mm256_srli_epi16(x, 8)), 1); return popcount_256_16(x); } inline __m256i popcount_256_32(__m256i x) { x = _mm256_sub_epi32(x, _mm256_and_si256(_mm256_srli_epi32(x, 1), _mm256_set1_epi8(0x55))); x = _mm256_add_epi32(_mm256_and_si256(x, _mm256_set1_epi8(0x33)), _mm256_and_si256(_mm256_srli_epi32(x, 2), _mm256_set1_epi8(0x33))); x = _mm256_and_si256(_mm256_add_epi8(x, _mm256_srli_epi32(x, 4)), _mm256_set1_epi8(0x0F)); return _mm256_srli_epi16(_mm256_madd_epi16(x, _mm256_set1_epi16(0x0101)), 8); } inline __m256i bsr_256_32_popcnt(__m256i x) { x = _mm256_or_si256(x, _mm256_srli_epi32(x, 1)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 2)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 4)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 8)); x = _mm256_srli_epi32(_mm256_or_si256(x, _mm256_srli_epi32(x, 16)), 1); return popcount_256_32(x); } inline __m256i bsr_256_64_popcnt(__m256i x) { x = _mm256_or_si256(x, _mm256_srli_epi64(x, 1)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 2)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 4)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 8)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 16)); x = _mm256_srli_epi64(_mm256_or_si256(x, _mm256_srli_epi64(x, 32)), 1); return popcount_256_64(x); }
シャッフル命令vpshufbを使った最適化 (rev popcount)
64ビット単位のときのみ、しかも若干効果がある程度ですが、vpshufbを使ってバイト順を反転させ、立っている一番下位のビットを算術演算で求めるテクニックを使うことができます。
inline __m256i bswap_256_64(__m256i x) { __m128i swap_table_128 = _mm_setr_epi8(7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8); __m256i swap_table = _mm256_broadcastsi128_si256(swap_table_128); return _mm256_shuffle_epi8(x, swap_table); } inline __m256i bsr_256_64_rev_popcnt(__m256i x) { x = _mm256_or_si256(x, _mm256_srli_epi64(x, 1)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 2)); x = _mm256_or_si256(x, _mm256_srli_epi64(x, 4)); x = _mm256_andnot_si256(_mm256_srli_epi64(x, 1), x); x = bswap_256_64(x); x = _mm256_and_si256(x, _mm256_sub_epi64(_mm256_setzero_si256(), x)); x = bswap_256_64(x); return popcount_256_64(_mm256_sub_epi64(x, _mm256_set1_epi64x(1))); }
ギャザー命令を使ったテーブル参照 (table gather)
32/64ビット単位で表引きができるvpgather{d,q}{d,q}を使ってpopcountを表引きします。32ビット単位なら4byte*231byte = 8GBに収まります。 表引きは000111...111の形のものしかないので実は32通りしかないので、L1キャッシュに乗る可能性が高いです。 しかもmallocでメモリを確保してもアクセスしていない領域には実際にはページ割当がされないので物理メモリ使用量はそこまで大きくなりません。
int32_t *table_32; void init_table_32() { table_32 = (int32_t*)malloc(sizeof(int32_t)*((size_t)1 << 31)); table_32[0] = 0; int64_t j = 1; for (int i = 1; i < 31; ++i) { table_32[j] = i; j = 2*j + 1; } } inline __m256i bsr_256_32_table_gather(__m256i x) { x = _mm256_or_si256(x, _mm256_srli_epi32(x, 1)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 2)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 4)); x = _mm256_or_si256(x, _mm256_srli_epi32(x, 8)); x = _mm256_srli_epi32(_mm256_or_si256(x, _mm256_srli_epi32(x, 16)), 1); return _mm256_i32gather_epi32(table_32, x, 4); }
vcvtdq2ps命令による浮動小数点数への変換 (cvtfloat)
浮動小数点数のビット表現は
などを参照してください。
vcvtdq2ps命令で浮動小数点数へ変換して、指数部だけになるようにシフトして、下駄履き(単精度なら127)を引いてやることでbsrが得られます。 しかしこれには罠があります。具体的には、
- ビットが浮動小数点数の精度以上に連続していると丸めによって間違った結果になる(1大きくなってしまう)。
- 最上位ビットは符号として扱われるので間違った結果になる(負数扱いになり全く違う値になる)。
これらへの対策については
或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化
或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化(その2)
を参照してください。
16ビット以下では対策は不要になります。 また、64ビット版では64ビットを3分割することで問題を回避したバージョンも示します。
inline __m256i bsr_256_32_cvtfloat_impl(__m256i x, int32_t sub) { __m256i cvt_fl = _mm256_castps_si256(_mm256_cvtepi32_ps(x)); __m256i shifted = _mm256_srli_epi32(cvt_fl, 23); return _mm256_sub_epi32(shifted, _mm256_set1_epi32(sub)); } inline __m256i bsr_256_8_cvtfloat(__m256i x) { __m256i r0 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi32(0x000000FF)), 127); __m256i r1 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi32(x, 8), _mm256_set1_epi32(0x000000FF)), 127); __m256i r2 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi32(x, 16), _mm256_set1_epi32(0x000000FF)), 127); __m256i r3 = bsr_256_32_cvtfloat_impl(_mm256_srli_epi32(x, 24), 127); __m256i r02 = _mm256_blend_epi16(r0, _mm256_slli_epi32(r2, 16), 0xAA); __m256i r13 = _mm256_blend_epi16(r1, _mm256_slli_epi32(r3, 16), 0xAA); return _mm256_blendv_epi8(r02, _mm256_slli_epi16(r13, 8), _mm256_set1_epi16(0xFF00)); } inline __m256i bsr_256_16_cvtfloat(__m256i x) { __m256i lo = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi32(0x0000FFFF)), 127); __m256i hi = bsr_256_32_cvtfloat_impl(_mm256_srli_epi32(x, 16), 127); return _mm256_blend_epi16(lo, _mm256_slli_epi32(hi, 16), 0xAA); } inline __m256i bsr_256_32_cvtfloat(__m256i x) { x = _mm256_andnot_si256(_mm256_srli_epi32(x, 1), x); // 連続するビット対策 __m256i result = bsr_256_32_cvtfloat_impl(x, 127); result = _mm256_or_si256(result, _mm256_srai_epi32(x, 31)); return _mm256_and_si256(result, _mm256_set1_epi32(0x0000001F)); } inline __m256i bsr_256_64_cvtfloat(__m256i x) { __m256i bsr32 = bsr_256_32_cvtfloat(x); __m256i higher = _mm256_add_epi32(_mm256_srli_epi64(bsr32, 32), _mm256_set1_epi64x(0x0000000000000020)); __m256i mask = _mm256_shuffle_epi32(_mm256_cmpeq_epi32(x, _mm256_setzero_si256()), 0xF5); return _mm256_blendv_epi8(higher, _mm256_and_si256(bsr32, _mm256_set1_epi64x(0xFFFFFFFF)), mask); } inline __m256i bsr_256_64_cvtfloat2(__m256i x) { __m256i bsr0 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi64x(0x3FFFFF)), 127); __m256i bsr1 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi64(x, 22), _mm256_set1_epi64x(0x3FFFFF)), 105); __m256i bsr2 = bsr_256_32_cvtfloat_impl(_mm256_srli_epi64(x, 44), 83); __m256i result = _mm256_max_epi32(_mm256_max_epi32(bsr0, bsr1), bsr2); return _mm256_and_si256(result, _mm256_set1_epi64x(0xFFFFFFFF)); }
vpshufbを使ったテーブル参照 (table)
バイト単位のシャッフル命令vpshufbは、バイト列の並べ替え以外に、4ビット入力、8ビット出力のルックアップテーブルとして用いることができます。 これを用いて、4ビット単位でbsrを求めて、それを組み合わせることで8ビット以上のbsrを求めることができます。
inline __m256i bsr_256_8_table(__m256i x) { __m128i table128_lo = _mm_setr_epi8(0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3); __m128i table128_hi = _mm_setr_epi8(0, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7); __m256i table_lo = _mm256_broadcastsi128_si256(table128_lo); __m256i table_hi = _mm256_broadcastsi128_si256(table128_hi); return _mm256_max_epi8( _mm256_shuffle_epi8(table_lo, _mm256_and_si256(x, _mm256_set1_epi8(0x0F))), _mm256_shuffle_epi8(table_hi, _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi8(0x0F)))); } inline __m256i bsr_256_16_table(__m256i x) { __m256i half = bsr_256_8_table(x); __m256i mask = _mm256_cmpgt_epi8(x, _mm256_setzero_si256()); __m256i res8 = _mm256_add_epi8(half, _mm256_and_si256(mask, _mm256_set1_epi16(0x0800))); return _mm256_max_epi16(_mm256_and_si256(res8, _mm256_set1_epi16(0x00FF)), _mm256_srli_epi16(res8, 8)); }
ベンチマーク
Haswell/Broadwell/Skylakeの各アーキテクチャのCPUで各アルゴリズムを並列に実行した時のの実行時間を計測しました。
bsrが合計231回実行されるように、8ビット幅なら226回、64ビット幅なら229回ループを回します。
環境
Haswell:
- Google compute engine n1-highmem-2 (2vCPU, 13GB mem)
- Xeon 2.3GHz, Turbo Boost時は不明
- OS: Ubuntu 16.10
- Compiler: GCC 6.1.1
Broadwell:
- Google compute engine n1-highmem-2 (2vCPU, 13GB mem)
- Xeon 2.2GHz, Turbo Boost時は不明
- OS: Ubuntu 16.10
- Compiler: GCC 6.1.1
Skylake:
コンパイルオプションはいずれも
-O3 -flto -mtune=(アーキテクチャ名) -march=(アーキテクチャ名) -lboost_system -lboost_timer
結果
(異なるアーキテクチャ同士で実行時間を比較するのはTurbo boostなどの影響があるのであまり意味がないと思います。)
8ビット
アルゴリズム | Haswell | Broadwell | Skylake |
---|---|---|---|
naive | 1.323s | 0.993s | 0.832s |
cvtfloat | 0.284s | 0.209s | 0.165s |
popcnt | 0.264s | 0.201s | 0.169s |
table | 0.079s | 0.058s | 0.049s |
いずれも表引きが一番高速なようです。
16ビット
アルゴリズム | Haswell | Broadwell | Skylake |
---|---|---|---|
naive | 1.562s | 1.192s | 0.948s |
cvtfloat | 0.257s | 0.200s | 0.146s |
popcnt | 0.550s | 0.423s | 0.337s |
table | 0.262s | 0.199s | 0.179s |
表引きあるいは浮動小数点数への変換が一番高速な結果になりました。
32ビット
アルゴリズム | Haswell | Broadwell | Skylake |
---|---|---|---|
naive | 1.581s | 1.220s | 0.993s |
cvtfloat | 0.416s | 0.318s | 0.237s |
popcnt | 1.278s | 0.966s | 0.701s |
table gather | 1.202s | 0.915s | 0.523s |
浮動小数点数への変換が圧倒的に速い結果になりました。 gather命令を使った表引きはSkylakeで若干高速になったようです。
64ビット
アルゴリズム | Haswell | Broadwell | Skylake |
---|---|---|---|
naive | 1.574s | 1.220s | 0.911s |
cvtfloat | 1.369s | 1.037s | 0.804s |
cvtfloat2 | 1.595s | 1.225s | 0.897s |
popcnt | 2.521s | 1.924s | 1.498s |
rev popcnt | 2.448s | 1.898s | 1.476s |
ここでも浮動小数点数への変換が一番速い結果となりました。3分割してコーナーケースの処理を不要にしたバージョンは思ったほど性能が出ませんでした。 64ビット幅になるとpopcntは命令数が増えてナイーブな実装より遅くなってしまいました。
来るAVX512時代には
少なくとも32ビット幅以上なら並列でclzを求めるvplzcnt{d,q}命令が追加されるのでそれを使えば良さそうです。 8ビット幅の場合vpermi2bを使って7ビットの表引きをして求めるのがいいかもしれません。 16ビット幅の場合はやってみないとわからなさそうです。
参考文献
或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化
或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化(その2)
bsrを並列処理する方法について考察した記事。当時はAVX2が使えるプロセッサはまだ発売されていなかった。
pclmulqdqを用いたビット単位のunpack【ビット演算テクニック Advent Calendar 2016 11日目】
この記事はビット演算テクニック Advent Calendar 2016
の11日目の記事です。
ビット単位のunpack
C++の固定長ビット配列を扱うクラスstd::bitset
を例として解説します。
std::bitset<128> unpacklo(const std::bitset<128> &a, const std::bitset<128> &b) { std::bitset<128> result; for (int i = 0; i < 64; ++i) { result[2*i] = a[i]; result[2*i+1] = b[i]; } return result; } std::bitset<128> unpackhi(const std::bitset<128> &a, const std::bitset<128> &b) { std::bitset<128> result; for (int i = 0; i < 64; ++i) { result[2*i] = a[i+64]; result[2*i+1] = b[i+64]; } return result; }
要は2つのビット列から交互に1ビットずつ取ってきたデータを生成します。 今回は簡単のため128ビット固定とします。
SSE
普通にstd::bitsetを使っていては遅いのでintelのx86の拡張命令セットであるSSE(Streaming SIMD Extensions)を利用します。 (追加された時期によりSSE, SSE2, ..., SSE4.2がありますがここではまとめてSSEと呼ぶことにします)
SSEを用いると128ビット単位で算術演算、ビット演算、シャッフルなど様々な操作を行うことができます。 SSEには8ビット単位のunpackを行うpunpacklbw/punpackhbwという命令があるので、まずこれを用いて8ビットごとにunpackし、残りを3日目に解説したDelta Swapを使ってunpackするのがぱっと思いつく実装でしょうか。 intrinsicsを使って実装してみます。
#include <x86intrin.h> __m128i mm_unpacklo_epb_unpack_dswap(__m128i a, __m128i b) { __m128i unpack8 = _mm_unpacklo_epi8(a, b); __m128i unpack4 = mm_delta_swap_epi64(unpack8, _mm_set1_epi16(0x00F0), 4); __m128i unpack2 = mm_delta_swap_epi64(unpack4, _mm_set1_epi8(0x0C), 2); return mm_delta_swap_epi64(unpack2, _mm_set1_epi8(0x22), 1); } __m128i mm_unpackhi_epb_unpack_dswap(__m128i a, __m128i b) { __m128i unpack8 = _mm_unpackhi_epi8(a, b); __m128i unpack4 = mm_delta_swap_epi64(unpack8, _mm_set1_epi16(0x00F0), 4); __m128i unpack2 = mm_delta_swap_epi64(unpack4, _mm_set1_epi8(0x0C), 2); return mm_delta_swap_epi64(unpack2, _mm_set1_epi8(0x22), 1); }
pdep
ビット演算に詳しい人ならBMIという拡張命令セットにpdepという命令があってこれが使えるということを思いつくかもしれません。
#include <x86intrin.h> __m128i mm_unpacklo_epb_pdep(__m128i a, __m128i b) { uint64_t alo = _mm_cvtsi128_si64(a); uint64_t blo = _mm_cvtsi128_si64(b); return _mm_set_epi64x( _pdep_u64(alo >> 32, UINT64_C(0x5555555555555555)) | _pdep_u64(blo >> 32, UINT64_C(0xAAAAAAAAAAAAAAAA)), _pdep_u64(alo, UINT64_C(0x5555555555555555)) | _pdep_u64(blo, UINT64_C(0xAAAAAAAAAAAAAAAA))); } __m128i mm_unpackhi_epb_pdep(__m128i a, __m128i b) { uint64_t ahi = _mm_extract_epi64(a, 1); uint64_t bhi = _mm_extract_epi64(b, 1); return _mm_set_epi64x( _pdep_u64(ahi >> 32, UINT64_C(0x5555555555555555)) | _pdep_u64(bhi >> 32, UINT64_C(0xAAAAAAAAAAAAAAAA)), _pdep_u64(ahi, UINT64_C(0x5555555555555555)) | _pdep_u64(bhi, UINT64_C(0xAAAAAAAAAAAAAAAA))); }
pclmulqdq
さて、また別の拡張命令セットCLMULにはpclmulqdqという命令があります。これは繰り上がりを無視しての2進数の掛け算を計算する命令で、CRCの計算やガロア体上の乗算を高速に行うために使われます。
pclmulqdqの動作をstd::bitset
を用いた疑似コードで書くとこんな感じでしょうか。
std::bitset<128> pseudo_pclmulqdq(const std::bitset<128> &a, const std::bitset<128> &b, const int index) { std::bitset<64> x, y; if (index & 0x01) { for (int i = 0; i < 64; ++i) { x[i] = a[i+64]; } } else { for (int i = 0; i < 64; ++i) { x[i] = a[i]; } } if (index & 0x10) { for (int i = 0; i < 64; ++i) { y[i] = b[i+64]; } } else { for (int i = 0; i < 64; ++i) { y[i] = b[i]; } } std::bitset<128> result; for (int i = 0; i < 64; ++i) { for (int j = 0; j < 64; ++j) { if (x[i] && y[j]) { result.flip(i+j); } } } return result; }
最初のindexにより分岐している部分で128ビットのうちどの64ビットを使って演算するかを選択し、64ビットx64ビットの繰り上がりなしの掛け算を行います。 繰り上がりなしの掛け算は有限体係数の多項式の乗算だと思うことができます。(足し算はXOR、掛け算はANDに相当するので)
ここで、の多項式の自乗を考えてみます。
これはつまり
for (int i = 0; i < 64; ++i) { result[2*i] = a[i]; }
と同じことになります。これを用いるとかなり簡潔にunpackを実装できます。
__m128i mm_unpacklo_epb_pclmulqdq(__m128i a, __m128i b) { return _mm_or_si128(_mm_clmulepi64_si128(a, a, 0x00), _mm_slli_epi32(_mm_clmulepi64_si128(b, b, 0x00), 1)); } __m128i mm_unpackhi_epb_pclmulqdq(__m128i a, __m128i b) { return _mm_or_si128(_mm_clmulepi64_si128(a, a, 0x11), _mm_slli_epi32(_mm_clmulepi64_si128(b, b, 0x11), 1)); }
ベンチマーク
以上3種類のアルゴリズムの実行コストを実際に計測してみました。
Hi/Loそれぞれ230回(≒10億回)unpackを計算しその時間を計測します。
1回の処理が軽いので、関数呼び出しやデータ生成が遅いと差が出る比較ができないので、(もともと最適化したら関数呼び出しは消えますが)関数にインライン化指定をして、データ生成はxorshiftを用いています。
ベンチマーク結果は以下のとおりです。
- CPU: Core i7-6700K 4.0GHz Turbo Boost 4.2GHz
- Memory: DDR-2133 64GB
- OS: Arch Linux (Linux 4.8.10-1-ARCH x86_64)
- Compiler: GCC 6.2.1
- Compiler Options: -O3 -flto -mtune=skylake -march=skylake -lboost_system -lboost_timer
アルゴリズム | unpacklo | unpackhi |
---|---|---|
8bitunpack+Delta Swap | 3.767s | 3.792s |
pdep | 2.624s | 2.938s |
pclmulqdq | 2.252s | 2.265s |
Hi/Loともにpclmulqdqを用いた場合が一番高速という結果になりました。 ちなみにpdepバージョンがHiとLoで実行時間に開きがあるのはmovq命令とpextrq命令の違いによると思われます。
pclmulqdq命令は使いどころが難しいですがスポッとハマると気持ちいいですね。
さいごに
ビット演算テクニックアドベントカレンダー 明日の12日目はocxtalさんによる「たぶんbitDPの話」です。お楽しみに!