できるXeonPhi(2) ！
はじめに
XeonPhiを導入したりvecAddしたりしたけど、それでも使いこなすにはまだまだだよね！
ということで、最適化について書いていきたいと思います
…実際、
http://software.intel.com/en-us/articles/intel-cilk-plus-aobench-sample
のコード解説以上の何かではないので、あらかじめお断りしておきます
題材
取り上げる題材は、みんな大好きAobenchです
https://code.google.com/p/aobench/
ふぁーすとすてっぷ
まずはCPUでどれぐらいの速度が出ているのかを見ます
計測に使うPCは、
あと、サイトに乗っている条件だと軽そうなので、条件をもう少し厳しくして
#define WIDTH 512 #define HEIGHT 512 #define NSUBSAMPLES 2 #define NAO_SAMPLES 16
にします。
あと適当に時間計測用関数を仕込んで、とりあえずシングルコアで実行してみましょう。そーれー
$ icc Aobench.c -O3 -lm -o aobench $ ./aobench time : 12.181054 sec.
何一つ並列化しないで12秒でした
あとでHaswellで測ってみるか…
一応XeonPhi上のコア1コアでも測ってみると
$ icc Aobench.c -O3 -lm -mmic -o aobench_mic $ scp aobench_mic mic0:~/ <- SCPでmic0にデータ送って $ ssh mic0 <- mic0にログインして $ ./aobench_mic <- 実行 time : 126.206486 sec.
あぁ、なんかこういう数字、Cell/B.E.のコード書いてる時によく見た気がする…
さて、ではまずどこが一番重いのかを特定します
アムダールの法則的に、そこを重点的に潰していくのが一番いいのでス
とりあえず適当にVtuneで速度を見ました
ambient_occlusion関数が重いです
どんな時に呼ばれるかというと
void render(unsigned char *img, int w, int h, int nsubsamples) { ... for (y = 0; y < h; y++) { for (x = 0; x < w; x++) { for (v = 0; v < nsubsamples; v++) { for (u = 0; u < nsubsamples; u++) { ... Isect isect; isect.t = 1.0e+17; isect.hit = 0; ray_sphere_intersect(&isect, &ray, &spheres[0]); ray_sphere_intersect(&isect, &ray, &spheres[1]); ray_sphere_intersect(&isect, &ray, &spheres[2]); ray_plane_intersect (&isect, &ray, &plane); if (isect.hit) { ambient_occlusion(&col, &isect); ... } } } } } }
width * height * sample * sampleの中で、飛ばしたレイがどこかにあたったときに計算するとそういうわけですね
# 言ってることが謎みたいな方は
# http://kagamin.net/hole/edupt/index.htm
# の物理ベースレンダラの資料を読んだりアンビエントオクルージョンについてググられるとそれなりに意味がわかるかと
まずambient_occlusion関数やるとか言った舌の根も乾かぬうちに、とりあえずこのループを何とかしてみましょう
ここで、上のループは全て独立であるということが一つ大きなポイントです
独立であるということは、それは別に並列で動作しても何も問題ないよね？ ということで並列に動かしてみます
まぁ言うまでもないことですが、OpenMPを使うときのループ変数はきちんとスレッドローカルになるようにしてあげないと大変なことになるので、頑張って変えてください
そして、変えた結果がこちら
$ icc Aobench.c -O3 -lm -openmp -std=gnu99 -mmic -o aobench_mic $ scp aobench_mic mic0:~/ $ ssh mic0 $ ./aobench_mic time : 44.997415 sec.
さっくり3倍ぐらい
とはいえ、CPUのシングルスレッドにまだ4倍近い差をつけられて負けているので、このままだとXeonPhiの存在意義に関わるわけで、さらに最適化を推し進めましょう
せかんどすてっぷ
今度こそambient_occlusion関数を見ます
void ambient_occlusion(vec *col, const Isect *isect) { ... int ntheta = NAO_SAMPLES; int nphi = NAO_SAMPLES; ... real occlusion = 0.0; for (j = 0; j < ntheta; j++) { for (i = 0; i < nphi; i++) { real theta = sqrt(drand48()); real phi = 2.0 * M_PI * drand48(); real x = cos(phi) * theta; real y = sin(phi) * theta; real z = sqrt(1.0 - theta * theta); real rx = x * basis[0].x + y * basis[1].x + z * basis[2].x; real ry = x * basis[0].y + y * basis[1].y + z * basis[2].y; real rz = x * basis[0].z + y * basis[1].z + z * basis[2].z; Ray ray; ray.org = p; ray.dir.x = rx; ray.dir.y = ry; ray.dir.z = rz; Isect occIsect; occIsect.t = 1.0e+17; occIsect.hit = 0; ray_sphere_intersect(&occIsect, &ray, &spheres[0]); ray_sphere_intersect(&occIsect, &ray, &spheres[1]); ray_sphere_intersect(&occIsect, &ray, &spheres[2]); ray_plane_intersect (&occIsect, &ray, &plane); if (occIsect.hit) occlusion += 1.0; } } occlusion = (ntheta * nphi - occlusion) / (real)(ntheta * nphi); ... }
一度あたったところが、どれぐらい遮蔽されているか？ というのを求めるのがこのコード
NAO_SAMPLE * NAO_SAMPLEの二乗のループになっています
今の条件だと、NAO_SAMPLE * NAO_SAMPLE = 16 * 16 = 256ですね
16とかすごくキリのいい数字なので、なんか並列化できそうですね。4の倍数とか2の倍数とかを見ると割とラッキーとかなります
でもIntrinsicsを直で書くのは面倒です…
そんな時はCilk+を使いましょう！
Cilk+を使うと何がイイか、というのは、下の記事を読めばよいかと思います
配列表記 (アレイ・ノーテーション) による SIMD 並列処理
http://www.isus.jp/article/parallel-special/simd-parallelism-using-array-notation/
コンパイラによる自動ベクトル化を行う上で必要になるあれこれをうまいこと何とかしてくれるいい感じのアノテーションです
どんなコードになったかは、ここに乗せるのは結構長くなってしまったので一部抜粋でお送りします。
詳細は後述のリポジトリのほうを見てください
void ambient_occlusion(vec *col, const Isect *isect) { ... real rand1[NAO_SAMPLES] __attribute__((aligned(64))); real rand2[NAO_SAMPLES] __attribute__((aligned(64))); real theta[NAO_SAMPLES] __attribute__((aligned(64))); real phi[NAO_SAMPLES] __attribute__((aligned(64))); real x[NAO_SAMPLES] __attribute__((aligned(64))); real y[NAO_SAMPLES] __attribute__((aligned(64))); real z[NAO_SAMPLES] __attribute__((aligned(64))); ... for(int j = 0; j < ntheta; ++j) { for(int k = 0; k < nphi; ++k) { rand1[k] = (real)drand48(); rand2[k] = (real)drand48(); } theta[:] = sqrt(rand1[:]); phi[:] = (real)2.0 * M_PI * rand2[:]; x[:] = cos(phi[:]) * theta[:]; y[:] = sin(phi[:]) * theta[:]; z[:] = sqrtf((real)1.0 - theta[:] * theta[:]); rx[:] = x[:] * basis[0].x + y[:] * basis[1].x + z[:] * basis[2].x; ry[:] = x[:] * basis[0].y + y[:] * basis[1].y + z[:] * basis[2].y; rz[:] = x[:] * basis[0].z + y[:] * basis[1].z + z[:] * basis[2].z; } .... }
real型はただfloatをtypedefしただけです
こうすると、確かにSIMD化されたコードが出力されます
確認、一応-Sオプションつけてコンパイルしてアセンブリ出力すると見えるんですが、ここに載せようとすると超長くなるのでやめます
さて、この状態だと？
$ ./aobench_cilkplus time = 23.9387 sec.
それでもまだCPUと比べると倍近いわけですね。さて、どうしたものか。
さーどすてっぷ
何が遅いんですかねぇ？
入念な調査の結果、遅い箇所がどこなのかを特定しました
# 実は関東GPGPU勉強会#2後、Haswell買いに行って始発を待つまでのファミレスの中で気が付いてその場でコードを書くというエクストリームなことをやりました
遅いのは、ここでした
for(int k = 0; k < nphi; ++k) { rand1[k] = (real)drand48(); rand2[k] = (real)drand48(); }
もっというと、drand48()でした
標準ライブラリのrand()とかは、XeonPhiではベクタライズされて速くなる、みたいなことを書いてあるのを見るのですが、そうでもない？ような感じです
ここで高速な乱数生成器であるXorShiftを導入してみます
XorShiftについてはググっていただきたいところですが、XORとビットシフトだけで構成された乱数生成器で、そこそこの周期でなかなか高速といった特徴を持ちます
コードについては、リポジトリのほうを参照してください
というわけで、XorShiftを導入した結果というのは
$ ./aobench_cilkplus time = 1.05928 sec.
以上のようになりました
まとめ
Aobenchを対象に、XeonPhiでの最適化についてちょっとだけ触れました
今回はIntelさんのAobenchのCilk+サンプルをちょっと改変してXeonPhiに載せて高速化していますが、Cilk+のアノテーションによるSIMD化は、実はそこまで効率が良いわけではないっぽいことを感じています
あと鬼門は標準ライブラリのrand()でした。最初全然気づきませんでした
XeonPhiではシングルスレッド比で約12倍の性能を達成しましたが、シングルスレッド比なのでフェアじゃないです
これをマルチスレッド化してSIMD化したらどうなるのか？ それについては、今ここで触れると泣きたくなるのでやめましょう
というわけで、リポジトリは以下です
https://bitbucket.org/telmin/aobench_cilkplus/
お試しになられたい方はどうぞ
ちなみにXeonPhi固有のコードというのはないので、普通のマルチコアCPUでも動作します
ただしコンパイルにiccが必要です
# Linuxだと非商用であればフリーのがあります
