プログラムを高速化する話Ⅱ ～GPGPU編～

GPUを利用して汎用演算を行う技術であるGPGPUを用いて、プログラムを高速化する技法についてまとめました。高速化の具体例も適宜用いて解説しています。
注: 京大マイコンクラブの春合宿2015で使ったスライドをそのまま掲載しているため、サンプルコードのリンクは利用できません。
CPU編はこちら https://speakerdeck.com/primenumber/puroguramuwogao-su-hua-suruhua

1. 〜GPGPU編〜 KMC5回生 prime プログラムを高速化する話Ⅱ
2. 自己紹介 KMC-ID: prime (KMC 5 回生 ) 京都大学理学部数学系 5 回生 ( 卒業確定 ) KMC での活動 : 競技プログラミング ゲーム AI 作成 難解プログラミング言語 電子工作 37 代会計 root, 電子錠の管理 etc.
3. この講座の内容 GPU: コンピューターの画像処理に使われる回路 graphics processing units の略 GPGPU: GPU を用いて様々な処理を行う General-purpose computing on GPU の略 GPU は CPU に比べて演算能力が高い CPU に比べ数倍〜数十倍高速なことも その性能を活かすためのテクニックを学ぶ 高速化の具体例もいくつか挙げていきます
4. もくじ はじめに その高速化、必要ですか？ GPGPU について CUDA 入門 CPU 編の復習 GPU 特有の高速化テクニック CPU とも共通する高速化テクニック
5. お断り このスライドでは最適化 = 高速化として扱います 高速なプログラムこそ最適 サンプルは CUDA で示します デファクトスタンダードなので サンプルコードは http://192.168.220.117:8000/ から見ることができます
6. はじめに その高速化、必要ですか？
7. その高速化、必要ですか？ 高速な既存のソフトウェア、ライブラリはないか クラウドコンピューティング等金の力で殴れないか より高速なアルゴリズムはないか これらのことをまず検討すべき
8. 最適化について 「細かい効率のことは忘れて、時間の 97% につ いて考え よう。時期尚早な最適化は諸悪の根源だ。 それでも残り 3% についても機会を逃すべきでは　 ない」 - Donald E. Knuth 「プログラム最適化の第一法則 : 最適化するな。 プログラム最適化の第二法則 ( 上級者限定 ): まだ するな。」 - Michael A. Jackson
9. 最適化について 最適化は、コードを複雑にすることが多いので、 コード の変更やデバッグを困難にする そのうえ、パフォーマンスに重大な影響を与える コード は全体のうちのほんの僅かなことが多い パフォーマンスに大きな影響を与えないコードを 最適化 してもほとんど意味がない
10. それでも どうしても限界まで高速化したい場合もある 莫大な計算量で金も時間も莫大にかかるとき ゲーム AI など高速であればあるほど有利な時 消費電力を減らしたいとき 利用できる計算資源が限られているとき
11. 前回の目標 コンパイラ・ライブラリ・ CPU 等の高速化技術 を上手く 活用して楽に高速化する アルゴリズムの改良 既存ライブラリ等を使用 など可能なことすべてをしても、必要な パフォーマンスが得られない時に、最 後の手段と して、手作業でプログラムを高速化するため に必 要なテクニックを学ぶ
12. 今回の目標 アルゴリズムの改良 既存ライブラリ等を使用 CPU 上のプログラムの最適化 など可能なことすべてをしても、必要なパフォーマ ンスが得られない時に、最後の最 後の手段として 、手作業で GPU プログラムを高速化するため に必 要なテクニックを学ぶ
13. CPU と GPU の性能比較 名前 Core i7- 8700K Core i9- 7890XE Ryzen Thread- Ripper 1950X Geforce GTX　 1080Ti RX VEGA 64 単精度 (GFLOPS) 710.4 2995.2 870.4 10596.8 10436.6 倍精度 (GFLOPS) 355.2 1497.6 435.2 331.15 652.3 理論値/いずれもベースクロックでの値
14. CPU と GPU の性能比較 名前 Xeon Platinum 8180 EPYC 7601 Tesla V100 単精度 (GFLOPS) 4480 1126.4 14899.2 倍精度 (GFLOPS) 2240 563.2 7449.6 CPU/GPUとも1個あたりの値 理論値/CPUはベースクロックでの値 Tesla V100はベースクロック不明のため ブーストクロックでの値
15. GPU はなぜ性能が高いのか？ 答え：演算器をたくさん積んでいるから CPU ではシングルスレッド性能向上等のために、 演算器以外の部分に大きく回路を割いている GPU では演算器をたくさん積むために、それ以外 の回路は控えめになっている
16. このスライドのベンチマークについて CPU は Core i7-6700K (Skylake) GPU は GeForce GTX 1080 (Pascal) にてそれぞれ計測を行った 備考 メンテの関係で CPU メモリはシングルチャネルで動作 GPU を接続している PCIe レーンは x8
17. GPU メーカーについて 主な GPU メーカーは 3 社 Intel CPU に内蔵された GPU(iGPU) のみ製造 AMD iGPU ・外付け (dGPU) ともに製造 NVIDIA dGPU のみ製造 CUDA は基本的に NVIDIA の GPU でのみ使える
18. GPGPU の基本 CPU と GPU は物理的に異なる 別々のプログラムが走る CPU GPU メモリ メモリ PCI Express/NVLink等
19. GPGPU の基本 CPU 用のメモリと GPU 用のメモリは別 データをやり取りするのに通信する必要がある CPU GPU メモリ メモリ PCI Express/NVLink等
20. GPGPU の基本 GPU メモリの特徴： CPU メモリの数倍帯域が広い GPU が利用される理由の１つ CPU GPU メモリ メモリ PCI Express/NVLink等
21. GPGPU の基本 CPU 上のプログラムがそのまま動くわけではない GPU 上で動くプログラムを作る必要 とはいえ、 CPU 上のプログラムとほぼ同じように書ける
22. GPGPU の基本 プログラムから GPU を扱う方法はいろいろある CUDA OpenCL OpenACC OpenMP (>4.0) DirectX DirectCompute OpenGL Compute Shader 今回はサンプルには CUDA を使う
23. CUDA とは NVIDIA が開発・提供している、 GPU 向けの汎用 並列コンピューティングプラットフォーム（並列コ ンピューティングアーキテクチャ）およびプログラ ミングモデルである。 - Wikipedia より CUDA C という言語を書くことで GPU 上で動く プログラムを書くことが出来る C++ とほぼ互換の言語 printf 等一部を除いて C++ 標準ライブラリは使えない
24. CUDA 入門 main.cu #include <cstdio> __global__ void kernel() { printf("Hello, World!\n"); } int main() { kernel<<<1, 1>>>(); cudaDeviceSynchronize(); return 0; }
25. CUDA 入門 コンパイルするには nvcc コンパイラを使う $ nvcc -o prog main.cu 実行 $ ./prog 実行結果 Hello, World!
26. CUDA 入門 CPU から呼ばれる GPU 側の関数（カーネル）には __global__ をつける 戻り値の型は void でなければならない 計算結果を CPU 側に伝えたい時はポインタを使う CPU 側から呼び出す時は kernel<<<1, 1>>>(); のように 三重山括弧を付けて呼び出す（数字の意味は後で） CPU から呼ばれる CPU 側の関数には __host__ 省略可能 GPU から呼ばれる GPU 側の関数には __device__
27. CUDA 入門 カーネルは呼び出したら計算が終わるのを待たずに 制御が戻ってくる GPU 側の計算がすべて終わるまで待つのに cudaDeviceSynchronize 関数を使う
28. CUDA 入門 変数 a, b の和を計算するコードを書く デバイス (GPU) 側コード __global__ void kernel(int * a, int * b, int * c) { *c = *a + *b; }
29. CUDA 入門 ホスト (CPU) 側コード int a = 1, b = 2; int *a_d, *b_d, *c_d; cudaMalloc((void**)&a_d, sizeof(int)); cudaMalloc((void**)&b_d, sizeof(int)); cudaMalloc((void**)&c_d, sizeof(int)); cudaMemcpy(a_d, &a, sizeof(int), cudaMemcpyHostToDevice); cudaMemcpy(b_d, &b, sizeof(int), cudaMemcpyHostToDevice); kernel<<<1, 1>>>(a_d, b_d, c_d); int c; cudaMemcpy(&c, c_d, sizeof(int), cudaMemcpyDeviceToHost); std::cout << c << std::endl;
30. CUDA 入門 GPU と CPU のメモリは違うので、 GPU 側メモリ を使いたいときには明示的に確保する必要がある CPU 側から cudaMalloc 関数を使う CPU-GPU 間でデータをコピーする必要がある CPU 側から cudaMemcpy 関数を使う CPU→GPU の転送には cudaMemcpyHostToDevice GPU→CPU の転送には cudaMemcpyDeviceToHost を指定する cudaMemcpy は同期的に実行されるので cudaDeviceSynchronize は必要ない
31. 例題：行列同士の積 N 行 M 列の行列 A と M 行 K 列の行列 B の積 を計算する 備考：既に高速なライブラリが存在 あくまで例題 (a ij )(b ij )=(∑ k=0 M a ik b kj )=(c ij )
32. 例題：行列同士の積 簡単のため、行列は 1 次元配列として確保 A, B, C をそれぞれ N*pitch1, M*pitch2, N*pitch3 個の長さを持つ配列とする pitch はコンピューターにとってキリの良い数 cudaMallocPitch 関数がやってくれる 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 A pitch1
33. 例題：行列同士の積 CPU 上で素直に動かすコード for (size_t i = 0; i < N; ++i) { for (size_t j = 0; j < K; ++j) { for (size_t k = 0; k < M; ++k) { mat3[i*pitch2 + j] += 　 mat1[i*pitch1 + k] * mat2[k*pitch2 + j]; } } }
34. 例題：行列同士の積 サンプルコード　デバイス (GPU) 側 __global__ void matmul_kernel_ver1( const float * const mat1, const float * const mat2, float * const mat3, const size_t pitch1, const size_t pitch2, const size_t pitch3, const size_t size1, const size_t size2, const size_t size3) { for (size_t i = 0; i < size1; ++i) { for (size_t j = 0; j < size3; ++j) { for (size_t k = 0; k < size2; ++k) { mat3[i*pitch3 + j] += mat1[i*pitch1 + k] * mat2[k*pitch2 + j]; } } } }
35. 例題：行列同士の積 サンプルコード　ホスト (CPU) 側 mat1_d 等は GPU メモリのアドレス matmul_kernel_ver1<<<1, 1>>>( mat1_d, mat2_d, mat3_d, pitch1, pitch2, pitch3, size1, size2, size3);
36. 例題：行列同士の積 10 100 1000 10000 0.01 0.1 1 10 100 1000 10000 100000 1000000 CPU(naive) CPU(optimized) GPU(naive) 行列のサイズ 実行時間 (ms)
37. 例題：行列同士の積 55 倍になった！！！
38. 例題：行列同士の積 55 倍になった！！！（実行時間が） GPU を使えば速くなるとはなんだったのか
39. GPGPU の高速化の大前提 並列化しろ！
40. GPGPU の高速化の掟① GPU は並列演算に特化したハードウェア 並列化しないと性能が出ない
41. 例題：行列同士の積 行列積の定義 結果の行列 C の i, j に関して依存なく並列化出来る 各スレッドに固有の i, j を割り当て、要素数分並列化 GPU では 100 万スレッドでも容易に走らせられる ハードウェアスレッドのおかげ (a ij )(b ij )=(∑ k=0 M a ik b kj )=(c ij )
42. 例題：行列同士の積 デバイス (GPU) 側 ( 引数部分は省略 ) const size_t i = blockIdx.y * blockDim.y + threadIdx.y; const size_t j = blockIdx.x * blockDim.x + threadIdx.x; if (i < size1 && j < size3) { for (size_t k = 0; k < size2; ++k) { mat3[i*pitch3 + j] += mat1[i*pitch1 + k] * mat2[k*pitch2 + j]; } }
43. 例題：行列同士の積 ホスト (CPU) 側 size_t bv = (size1+threadsPerBlock-1) / threadsPerBlock; size_t bh = (size3+threadsPerBlock-1) / threadsPerBlock; dim3 block(threadsPerBlock, threadsPerBlock); dim3 grid(bh, bv); matmul_kernel_ver2<<<grid, block>>>( mat1_d, mat2_d, mat3_d, pitch1, pitch2, pitch3, size1, size2, size3);
44. 例題：行列同士の積 ホスト (CPU) 側 size_t bv = (size1+threadsPerBlock-1) / threadsPerBlock; size_t bh = (size3+threadsPerBlock-1) / threadsPerBlock; dim3 block(threadsPerBlock, threadsPerBlock); dim3 grid(bh, bv); matmul_kernel_ver2<<<grid, block>>>( mat1_d, mat2_d, mat3_d, pitch1, pitch2, pitch3, size1, size2, size3);
45. 例題：行列同士の積 カーネル呼び出しの三重山括弧の中で並列度を指定 する 1 ブロックあたり 8*8 スレッドを、 bv*bh ブロック同時に走らせる 各スレッドで同じプログラムが実行される SIMT(Single Instruction Multiple Thread) という 自分がどのスレッドかは threadIdx, blockIdx で知れる
46. 例題：行列同士の積 10 100 1000 10000 0 0.01 0.1 1 10 100 1000 10000 100000 1000000 CPU(naive) CPU(optimized) GPU(naive) GPU(parallel) 行列のサイズ 実行時間 (ms)
47. 例題：行列同士の積 とりあえず 6000 倍速くなった まだ CPU で最適化したものに比べると遅い 行列のサイズが大きいとき このあといろいろな方法で高速化していく
48. GPGPU で高速化するには 並列化する CPU と共通の高速化テクニック GPU 特有の高速化テクニック これらを組み合わせて高速化する
49. CPU 編の復習 キャッシュの活用 キャッシュに収まるように工夫することで高速化 ビット演算の活用 ビットレベルの並列性を生かして高速化 SIMD(single instruction multiple data) 命令の活用 1 つの命令で複数データを処理して高速化
50. CPU 編の復習：キャッシュの活用 メインメモリは遅い アクセスに数十〜数百 cycle かかる 帯域も CPU の演算能力に対して十分でない CPU コアに近いキャッシュを利用 容量は小さいが高速 L1 〜 L3 など容量と速度により階層化されている
51. CPU 編の復習：キャッシュの活用 具体的なテクニック： 局所的なメモリアクセス SoA なデータ構造 ストリップマイニング ブロック化
52. CPU 編の復習：ビット演算の活用 ビット演算とは 2 進数の 0/1 列を操作する演算の総称 ビット論理積 ビット論理和 ビット排他的論理和 ビット否定 ビットシフト ビット操作用のその他の命令 加減乗算命令を使うことも
53. CPU 編の復習：ビット演算の活用 なぜビット演算は高速なのか？ ビット演算命令自体が高速 一度にビット幅分並列実行できる メモリ使用量を減らしてキャッシュヒット率の向上 ビット配列： 整数等の配列を 0/1 の配列として利用するテクニック
54. CPU 編の復習：ビット演算の活用 ビット列への操作を行う様々なテクニックが存在 特定ビットの操作・マスク popcount ・ハミング距離 立っている一番上・下のビットを求める 立っているビットを走査する ビット列の並びを反転する 部分集合の列挙 ビット列の一部をスワップ magic bitboard
55. CPU 編の復習： SIMD 命令の活用 SIMD 命令とは： 1 つの命令で複数のデータに対し一括で処理を行う命令 A[0] A[1] A[2] A[3] B[0] B[1] B[2] B[3] SIMD幅 + A[0]+B[0] A[1]+B[1] A[2]+B[2] A[3]+B[3]
56. CPU 編の復習： SIMD 命令の活用 1 命令で複数データを処理できるので高速 注意点・テクニック： メモリアラインメントを合わせたほうが良い メモリアクセスが律速になりやすい 条件分岐をマスクを用いてベクトル化する 単純な SIMD にならない場合は水平可算やシャッフル
57. CPU 編の復習： GPGPU 編との関連 これらのテクニックは GPGPU でも有効 キャッシュ： GPU のキャッシュは小さいが有効 ビット演算： CPU 同様速い SIMD 命令： GPU は大体 SIMD アーキテクチャ
58. NVIDIA の GPU の構造 GPU は複数の SM(Streaming Multiprocesser) からなる 各 SM は複数の CUDA Core からなる GPU SM CUDA Core
59. CUDA の実行単位 1 個以上のスレッドからなる「ブロック」 1 個以上のブロックからなる「グリッド」 CPU からはグリッド単位で実行する 1 ブロックあたり 64 スレッド、 1 グリッドあたり 100 ブロックなら 64*100=6400 スレッドが走ることになる
60. SM の構造 CUDA Core の組が複数と、 L1 Cache, Shared Memory などからなる（関係ない部分は省略） L1 Cache Shared Memory
61. GPU の性能と SM GPU の性能は SM の数とクロックとメモリ帯域等 で決まる GPU SMの数 GTX1080Ti 28 GTX1080 20 GTX1070Ti 18 GTX1070 15 GTX1060 10 GTX1050Ti 6 GTX1050 5 GT1030 3
62. SM とブロック 1 つのブロックに入っているスレッドは必ず同じ SM の中で実行される 1 つの SM は 1 つ以上のブロックを同時に実行する
63. Warp ブロックの中では 32 スレッドごとに塊を作る Warp という Warp の中ではスレッドは基本的に同じ命令を実行 SM の中の CUDA Core の組が同時に動いて実行 各スレッドが SIMD の各レーンに対応するようなもの Volta アーキテクチャだとちょっと違う ブロック内のスレッド数が 32 で割り切れないと 余った演算機は休んでしまう スレッド数は 32 の倍数が望ましい
64. Warp スレッド 1 スレッド 2 スレッド 3 スレッド 4 スレッド 5 スレッド 6 … スレッド 32 命令A 命令A 命令A 命令A 命令A 命令A 命令A 命令A 命令B 命令B 命令B 命令B 命令B 命令B 命令B 命令B 命令C 命令C 命令C 命令C 命令C 命令C 命令C 命令C 命令D 命令D 命令D 命令D 命令D 命令D 命令D 命令D
65. GPGPU で性能を出すコツ GPU 特有のテクニック CPU-GPU 間の通信を抑える コアレスアクセス シェアードメモリ・コンスタントメモリの活用 同時マルチスレッディングと Occupancy ワープ内分岐を抑える ストリームによるタスク並列化
66. • CPU-GPU 間の通信を抑える CPU-GPU 間の帯域はそれほど広くない PCIe3.0x16 で 16GB/s CPU-GPU 間の通信をなるべく少なくする CPU GPU メモリ メモリ PCI Express/NVLink等
67. コアレスアクセス メモリアクセスは Warp ごとにまとめられる A B C B A B C A A B C (A, B, Cのアドレスが離れていたら) 3回のメモリアクセスになる
68. コアレスアクセス 近い場所にアドレスが固まっている場合、１つの　 メモリトランザクションにまとめられることがある アクセスがまとめられると、メモリトランザクションが 減り、性能が向上する
69. コアレスアクセス 32 バイト単位等 GPU にとってキリのいい単位の 中に、一つのワープからのメモリアクセスが 複数ある場合、ひとつのトランザクションに まとめられる Warp 64 68 72 76 80 84 88 92 96 60 アドレス
70. 例題：二次元配列のコピー 二次元配列の内容をコピーする cudaMemcpy2D で出来るが、例題ということで… cudaMallocPitch で確保した二次元配列 幅 width, 高さ height, ストライド pitch
71. 例題：二次元配列のコピー 縦にコピーと横にコピーで比較する
72. 例題：二次元配列のコピー 1000 10000 100000 0.01 0.1 1 10 100 1000 GPU(column-major) GPU(row-major) width = height 実行時間 (ms)
73. 例題：二次元配列のコピー 横にコピーするとコアレスアクセスになるので高速 縦にコピーするとコアレスアクセスにならず低速
74. シェアードメモリの活用 各ブロック内でしか読み書きできないメモリ 容量は小さいが低遅延 & 広帯域 うまく使えれば高速化
75. 例題：行列同士の積 シェアードメモリを使って高速化する そもそもなぜ行列積の性能がでなかったか？ メモリ帯域の制限があるから
76. メモリ帯域と演算能力 GTX1080 の演算能力は約 9TFLOPS 1 秒間に 9 兆回計算できる 実際には 4.5 兆回の積和算 (A*B+C) ができる foat は 4byte 、積和算は 3 つの値を読み込む 単純に考えると 4bytes*3*4.5T/s=54TB/s 必要 一方メモリ帯域は 320GB/s 全然足りない！
77. 例題：行列同士の積 行列が N 行 N 列の正方行列とすると、 読み込む必要があるのは A, B それぞれ N^2 個 計算する必要があるのは約 2*N^3 回 うまくキャッシュ出来れば、高速化出来る ここではブロック化という方法を使う
78. 例題：行列同士の積 答えの行列 C の i 行 j 列を計算するとき、 = × c ij a i 1 b 1 j 1回の積和算に2回読み込み
79. 例題：行列同士の積 同時に i+1 行 j 列を計算すると、 = × c ij a i 1 b 1 j 2回の積和算に3回読み込み c i+1 j a i+11
80. 例題：行列同士の積 列に関しても同じことをすると、 = × c ij a i 1 b 1 j 4回の積和算に4回読み込み c i+1 j a i+11 c ij+1 c i+1 j +1 b 1 j +1
81. 例題：行列同士の積 今回は 8x8 のブロック化を行う 64 回の積和算に 16 回の読み込みで済む コアレスアクセスのため、ループ 8 回分の 読み込みを一度に行う = × ループ8回分 8 ループ8回分 8 8 8
82. 例題：行列同士の積 サンプルのため端数処理は省略 __shared__ float localA[threadsPerBlock][threadsPerBlock]; __shared__ float localB[threadsPerBlock][threadsPerBlock]; float tmp = 0.0; for (size_t k = 0; k < size2; k += threadsPerBlock) { __syncthreads(); localA[threadIdx.y][threadIdx.x] = mat1[i*pitch1 + (k + threadIdx.x)]; localB[threadIdx.y][threadIdx.x] = mat2[(k + threadIdx.y)*pitch2 + j]; __syncthreads(); for (size_t k2 = 0; k2 < threadsPerBlock; ++k2) { tmp += localA[threadIdx.y][k2] * localB[k2][threadIdx.x]; } } mat3[i*pitch3 + j] = tmp;
83. 例題：行列同士の積 __syncthreads 関数 ブロック内で同期をとる関数 ブロック内のすべてのスレッドがここに到達する まで待機する すべてのスレッドがシェアードメモリに書き込み 終わるのを待っている
84. 例題：行列同士の積 同期をとらないと シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリから 読み込み シェアードメモリから 読み込み 意図した動作 スレッドA スレッドB
85. 例題：行列同士の積 同期をとらないと シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリから 読み込み シェアードメモリから 読み込み スレッドA スレッドB
86. 例題：行列同士の積 同期をとらないと シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリへの 書き込み シェアードメモリから 読み込み シェアードメモリから 読み込み 意図しない動作 スレッドA スレッドB ←まだ書き込みしていない ところにアクセス
87. 例題：行列同士の積 10 100 1000 10000 0 0.01 0.1 1 10 100 1000 10000 100000 CPU(optimized) GPU(parallel) GPU(shared memory) 行列のサイズ 実行時間 (ms)
88. 例題：行列同士の積 数倍高速化、 CPU よりも速くなった とはいえ、理論性能からするとまだまだ
89. シェアードメモリのバンクコンフリクト シェアードメモリにはバンクという概念がある バンク毎に 1 回に 1 個のメモリアクセスしかできない バンクは 32 個ある バンクコンフリクトなし １回のアクセスでOK
90. シェアードメモリのバンクコンフリクト 各ワード (4byte) は順番に各バンクに割り振られて いる バンク 0 1 2 3 4 5 6 7 ワード 0 1 2 3 4 5 6 7 ワード 32 33 34 35 36 37 38 39 ワード 64 65 66 67 68 69 70 71 …
91. シェアードメモリのバンクコンフリクト バンクコンフリクトするとコンフリクトの回数 メモリアクセスが必要になり、遅くなる アクセスパターンに気をつけてコンフリクトを 回避すると高速化出来る バンクコンフリクトあり 最大３回アクセス必要
92. シェアードメモリのバンクコンフリクト アクセスパターンに気をつけてコンフリクトを 回避すると高速化出来る バンクコンフリクトなし
93. バンクコンフリクトとパディング バンクの数が 4 だとして説明する 二次元配列に縦にアクセス するなど、 4 個飛びで アクセスすると、 バンクコンフリクトする 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3
94. バンクコンフリクトとパディング 末尾にパディングを入れると、バンクコンフリクト を回避できる __shared__ T ary[N][M+1]; 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3
95. コンスタントメモリの活用 各 SM にコンスタントキャッシュがある コンスタントメモリが使えるとき GPU から値を書き換えない 64KB 以内 Warp 内で同じ値を使うときに高速になる そうでないときはあまり速くない シェアードメモリ同様うまく使えれば高速 π などの定数を格納するのに使える
96. 同時マルチスレッディング 1 つの SM は複数のワープを持つ それらのワープのうち、実行可能なものを選んで　 実行する SM ワープ
97. 同時マルチスレッディング 1 つの SM が複数の Warp を持つ それらの Warp のうち、実行可能なものを選んで　 実行する SM 実行
98. 同時マルチスレッディング 同時マルチスレッディングによって、ある Warp の 命令が終わるのを待っている間に他の Warp を実行 するなどして、命令のレイテンシを隠蔽できる レイテンシをなるべく隠蔽するために、なるべく たくさんの Warp を同時に走らせたい 同時に走らせられる Warp 数は、ハードウェア上の リソースにより制限されている
99. 同時実行できる Warp 数の制限 SM 内で共有されるリソース レジスタファイル シェアードメモリ 例えば GTX1080 ではシェアードメモリは 96KB/SM 1Warp あたり 10KB なら同時実行できるのは最大 9Warp 1 ブロックあたりシェアードメモリの使用量を S 、 レジスタファイル使用量を R 、ブロックあたりの Warp 数を W とすると min(64, 32*W, 96K/S, 256K/R)Warp 同時実行可能
100. Occupancy 実際に同時実行可能な Warp 数 /64( 最大値 ) を Occupancy （占有率）という Occupancy を高めるのが効率を上げる方法の１つ スレッド当たりシェアードメモリの使用量を減らすなど
101. ブロックサイズ・グリッドサイズの調整 ブロックサイズにより Occupancy が変わってくる 通常は 128 〜 256 スレッドが良いといわれている グリッドサイズは SM の数より十分多くなるように 各 SM をたくさんのブロックで占有するため
102. Warp 内分岐を抑える Warp 内で分岐する・しない両方のスレッドがある 場合、両方のパスをマスク付きで実行する 分岐前 分岐しない 分岐する
103. Warp 内分岐を抑える 「分岐しない」パスを実行しているとき、 「分岐する」スレッドの演算器はずっと休んでいる 逆に、「分岐する」パスを実行しているとき、 「分岐しない」スレッドの演算器はずっと休んでいる 演算器が休んでいる分、効率が落ちる
104. Warp 内分岐を抑える 各スレッドでバラバラに分岐すると更に酷いことに
105. Warp 内分岐を抑える すべてのパスを実行するとその分時間がかかる できるだけ Warp 内で分岐する・しないを統一
106. 例題： N-Queen N*N のマス目上に N 個のクイーンを互いに効きが ないように置く場合の数を求める Q ＼ ー ー ー Q Q ／ ／ Q ＼ ｜ ／ Q ー ー Q ー ー ー Q ／ ｜ ＼ ／ ｜ ＼ Q ｜ Q ＼ Q ｜ N=4: 2通り N=6: 8通り
107. 例題： N-Queen 縦の効きを考えると各列に置けるクイーンは 1 個 同様に、各行に置けるのも１個 斜めの効きを考えなければ 組み合わせは n! 通りある 1 列ずつ決めていって、 可能な解を探索できる バックトラックという | ー Q ー ー Q ー | | | Q ダメな例
108. 例題： N-Queen 探索するとき、既に置いた駒の位置を記憶しなくて も、効きだけ持っておけば探索できる 効きも、右、右上、右下 だけ持っておけば良い 効きをビット列で管理 ビット演算で処理できる ／ ／ ／ Q ー ー ー Q ー ー ー ー ー ＼ ／ ＼ ☓ ＼ Q ー ー ー ー
109. 例題： N-Queen CPU( シングルスレッド ) 版 uint64_t solve(int N, int depth = 0, uint32_t left = 0, uint32_t mid = 0, uint32_t right = 0) { 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; }
110. 例題： N-Queen 再帰的に探索する Q × × Q × × Q × Q × Q × × Q × Q × × Q × ×
111. 例題： N-Queen 再帰的に探索する × Q × × Q × × Q × Q × Q × Q × Q Q Q Q
112. 例題： N-Queen GPU 版では並列化のため、数列分展開したものを 並列に解き、最後にその合計を求めることにする Q Q Q Q Q Q 1通り 1通り 0通り
113. 例題： N-Queen CPU 版も同様に並列化出来る コードは省略
114. 例題： N-Queen GPU 版 __global__ void kernel(const int N, const int depth, const uint32_t * const left_ary, const uint32_t * const mid_ary, const uint32_t * const right_ary, uint64_t * const result_ary, const size_t size) { int index = threadIdx.x + blockIdx.x * blockDim.x; if (index < size) { result_ary[index] = solve(N, depth, left_ary[index], mid_ary[index], right_ary[index]); } }
115. 例題： N-Queen 9 10 11 12 13 14 15 16 17 18 0 0 0.01 0.1 1 10 100 CPU(naive) CPU(parallel) GPU(naive) N 実行時間 (s)
116. 例題： N-Queen GPU 版は CPU のシングルスレッド版よりは速いが マルチスレッド版よりは遅い
117. 例題： N-Queen 再帰 + ループで解いていた ループの終了条件で分岐している しかも、分岐したあと再帰でさらに分岐する
118. 例題： N-Queen 再帰的な分岐で演算器が休んでいる 分岐 更に分岐
119. 例題： N-Queen 解決方法 再帰をやめ、自前でスタックを持ってエミュレート スタックはシェアードメモリに置ければベスト アクセスが高速なので 容量制限で置けないときはローカルメモリに置く
120. 例題： N-Queen solve 関数 __host__ __device__ uint64_t solve(int N, int depth = 0, uint32_t left = 0, uint32_t mid = 0, uint32_t right = 0) { 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; }
121. 例題： N-Queen スタックで持つべき状態 depth, left, mid, right, pos の 5 つ 再帰とループと同じことを１つのループで実現する 各ループで、 pos から 1bit とり、スタックに新しい状態を push pos から取れるビットがなければ pop スタックが空になったら終了
122. 例題： N-Queen 例 depth left mid right pos 2 0b1000 0b1010 0b0100 0b0001
123. 例題： N-Queen 例 depth left mid right pos 2 0b1000 0b1010 0b0100 0b0001
124. 例題： N-Queen 例 depth left mid right pos 3 0b0010 0b1011 0b0010 0b0100 2 0b1000 0b1010 0b0100 0b0000
125. 例題： N-Queen 例 depth left mid right pos 3 0b0010 0b1011 0b0010 0b0100 2 0b1000 0b1010 0b0100 0b0000
126. 例題： N-Queen 例 depth left mid right pos 4 0b1100 0b1111 0b0011 0b0000 3 0b0010 0b1011 0b0010 0b0000 2 0b1000 0b1010 0b0100 0b0000
127. 例題： N-Queen 例 depth left mid right pos 4 0b1100 0b1111 0b0011 0b0000 3 0b0010 0b1011 0b0010 0b0000 2 0b1000 0b1010 0b0100 0b0000 ++Count
128. 例題： N-Queen 例 depth left mid right pos 3 0b0010 0b1011 0b0010 0b0000 2 0b1000 0b1010 0b0100 0b0000
129. 例題： N-Queen 例 depth left mid right pos 2 0b1000 0b1010 0b0100 0b0000
130. 例題： N-Queen 例 depth left mid right pos
131. 例題： N-Queen 例 depth left mid right pos スタックが空になったので終了
132. 例題： N-Queen ソースコードの全体は https://github.com/primenumber/nqueen
133. 例題： N-Queen 9 10 11 12 13 14 15 16 17 18 0.0001 0.0010 0.0100 0.1000 1.0000 10.0000 100.0000 CPU(naive) CPU(parallel) GPU(naive) GPU(optimized) N 実行時間 (s)
134. ストリームによるタスク並列 CUDA ではストリームという機能を用いて、タスク （カーネルやメモリ転送）並列の並列化ができる
135. ストリーム ストリームごとにキューがあるイメージ 各キューの前のタスクが終わると次のタスクが実行 ストリーム1 ストリーム2 タスクA タスクB タスクC タスクD タスクE タスクF GPU
136. ストリーム ストリームを作るには cudaStreamCreate 関数 ストリームのキューが空になるまで待つには cudaStreamSynchronize 関数 ストリームを削除するには　　 cudaStreamDestroy 関数
137. メモリ転送とカーネル実行の並列化 カーネルを複数のストリームに分割すると、メモリ 転送とカーネル実行を同時に行える メモリ転送(HtoD) メモリ転送(DtoH) カーネル実行 HtoD DtoH カーネル HtoD DtoH カーネル 高速化
138. メモリ転送とカーネル実行の並列化 メモリ転送を非同期に行うには cudaMemcpyAsync を使う その際、ホスト側のアドレスには cudaMallocHost で確保したアドレスを渡す カーネルをあるストリームで実行するには、三重山 括弧の第 4 引数にストリームを渡す
139. メモリ転送とカーネル実行の並列化 cudaStream_t str[3]; for (int i = 0; i < 3; ++i) { cudaStreamCreate(str + i); cudaMemcpyAsync(a_d + n*i, a_h + n*i, sizeof(int) * n, cudaMemcpyHostToDevice, str[i]); kernel<<<1024, 256, 0, str[i]>>>(a_d + n*i, b_d + n*i, n); cudaMemcpyAsync(b_h + n*i, b_d + n*i, sizeof(int) * n, cudaMemcpyDeviceToHost, str[i]); } for (int i = 0; i < 3; ++i) { cudaStreamSynchronize(str[i]); cudaStreamDestroy(str[i]); }
140. GPGPU で性能を出すコツ CPU にも共通する部分があるテクニック ソフトウェアパイプライニング レジスタブロッキング アトミック命令
141. ソフトウェアパイプライニング for (int i = 0; i < n; ++i) { A(i); B(i); } 各 i に対し A(i), B(i) に依存関係があるとする A(i) と A(i+1), A(i+2),... 、 B(i) と B(i+1), B(i+2),... には依存関係がないとする A(i), B(i) にそれなりのレイテンシがあるとすると
142. A(0) B(0) A(1) B(1) A(2) B(2) A(3) B(3) A(4) B(4)
143. ソフトウェアパイプライニング for (int i = 0; i < n; ++i) { A(i); } for (int i = 0; i < n; ++i) { B(i); } とすることで
144. ソフトウェアパイプライニング A(0) B(0) A(1) B(1) A(2) B(2) A(3) B(3) A(4) B(4)
145. ソフトウェアパイプライニング レイテンシが隠蔽され高速化された！ 命令を実行できるものから順不同で実行してくれる （アウトオブオーダー） CPU ならあまり配慮は　 いらない GPU は Warp 内はインオーダーで実行する 多くは同時マルチスレッディングで隠蔽される Occupancy が上げられないときに効果あり
146. レジスタブロッキング レジスタ上のデータを再利用できるようにレジスタ にデータを置くこと CPU でも当然有効だが、 GPU は CPU に比べて　 レジスタが豊富 ただし、レジスタは SM 内で共有される 大量に使うと Occupancy が低下する
147. 例題：行列同士の積 シェアードメモリを使っているので既に Occupancy は低い ソフトウェアパイプライニングするために、 1 スレッドで 8x8 要素分の答えを同時に計算する 同時にレジスタブロッキングでメモリアクセスを減らす ブロック全体では 64*64 のブロック化になる 4096 回の積和算に対し 128 回の読み込み ループ 8 回分の読み込みを同時に行うことで、コアレス アクセスにしてさらに高速化
148. 例題：行列同士の積 ブロック化の様子 64 64 8 8 64 64 = ×
149. 例題：行列同士の積 ソースコード https://github.com/primenumber/matmul
150. 例題：行列同士の積 10 100 1000 10000 0 0.01 0.1 1 10 100 1000 10000 100000 CPU(optimized) GPU(parallel) GPU(shared memory) GPU(register blocking) 行列のサイズ 実行時間 (ms)
151. 例題：行列同士の積 （行列が大きいときに）更に数倍以上速くなった N=4096 で理論値の 60% ほど ここからさらに倍は速くならない
152. アトミック命令 次の操作を考える Load A to x x = x+1 Store x to A 操作の後 A の値は 1 増える これを二並列で実行したら、当然 2 増えてほしい
153. アトミック命令 マルチスレッドで同じ場所に対して読み書きすると Load A to x スレッド1 Load A to y x = x+1 y = y+1 Store x to A Store y to A スレッド2
154. アトミック命令 マルチスレッドで同じ場所に対して読み書きすると Load A to x スレッド1 Load A to y x = x+1 y = y+1 Store x to A Store y to A スレッド2 Aの値は2増える
155. アトミック命令 マルチスレッドで同じ場所に対して読み書きすると Load A to x スレッド1 Load A to y x = x+1 y = y+1 Store x to A Store y to A スレッド2
156. アトミック命令 マルチスレッドで同じ場所に対して読み書きすると Load A to x スレッド1 Load A to y x = x+1 y = y+1 Store x to A Store y to A スレッド2 Aに対して2回足したはず なのに1しか増えない！
157. アトミック命令 Read-modify-write などの操作を他のスレッドに　 邪魔されずに一度に行う命令 modify には加算、乗算などの処理が入る
158. アトミック命令 アトミック命令を使うと一貫性が保証される Load A to x スレッド1 Load A to y x = x+1 y = y+1 Store x to A Store y to A スレッド2 Atomic Atomic
159. アトミック命令 アトミック命令を使うと並列化の幅が広がる GPU の性能を引き出すのに役立つ 一方、同じ場所に集中してアクセスすると一貫性の 保証にコストがかかり遅くなる原因に 次の例題で扱う
160. 例題：合計値の計算 多数の数の合計値の計算を行う 5 8 3 4 2 9 1 2 6 2 合計 42
161. 例題：合計値の計算 簡単のため、 1 〜　までの合計値を計算する 簡単に式で書けることには目を瞑る アトミック命令で全部足してしまえば計算できる 2M
162. 例題：合計値の計算 デバイス側コード __global__ void sum_vec_naive( unsigned long long * result) { int i = threadIdx.x + blockIdx.x * threadsPerBlock; atomicAdd(result, i + 1); }
163. 例題：合計値の計算 20 21 22 23 24 25 26 27 28 29 30 1 10 100 1000 GPU(naive) CPU(naive) M 実行時間 (ms)
164. 例題：合計値の計算 CPU より遅い！ アトミックになっている場所は並列に実行できないので ブロック内で先に和を取り、結果を足すようにする 5 8 3 2 9 1 6 4 6 9 7 3 2 1 8 5 result
165. 例題：合計値の計算 CPU より遅い！ アトミックになっている場所は並列に実行できないので ブロック内で先に和を取り、結果を足すようにする 5 8 3 2 9 1 6 4 6 9 7 3 2 1 8 5 result tmp tmp
166. 例題：合計値の計算 同様にシェアードメモリに対してアトミック命令で 加算すると、ブロック内で並列化できない reduction を行ってブロック内の和を得る 5 8 3 2 9 1 6 4 13 5 10 10 ＋ ＋ ＋ ＋ 18 20 ＋ ＋ ＋ 38
167. 例題：合計値の計算 int i = threadIdx.x + blockIdx.x * threadsPerBlock; __shared__ unsigned long long tmp[threadsPerBlock]; tmp[threadIdx.x] = i + 1; __syncthreads(); for (uint32_t s = 1; s < threadsPerBlock; s *= 2) { if (threadIdx.x % (2*s) == 0) { tmp[threadIdx.x] += tmp[threadIdx.x + s]; } __syncthreads(); } if (threadIdx.x == 0) { atomicAdd(result, tmp[0]); }
168. 例題：合計値の計算 20 21 22 23 24 25 26 27 28 29 30 0.1 1 10 100 1000 GPU(naive) CPU(naive) GPU(reduction) M 実行時間 (ms)
169. 例題：合計値の計算 CPU よりは速くなった 実はまだ速く出来る まず、シェアードメモリへのアクセスが バンクコンフリクトしているので、それを解消する
170. 例題：合計値の計算 int i = threadIdx.x + blockIdx.x * threadsPerBlock; __shared__ unsigned long long tmp[threadsPerBlock]; tmp[threadIdx.x] = i + 1; __syncthreads(); for (uint32_t s = threadsPerBlock/2; s > 0; s >>= 1) { if (threadIdx.x < s) { tmp[threadIdx.x] += tmp[threadIdx.x + s]; } __syncthreads(); } if (threadIdx.x == 0) { atomicAdd(result, tmp[0]); }
171. 例題：合計値の計算 変更前 5 8 3 2 9 1 6 4 13 5 10 10 ＋ ＋ ＋ ＋ 18 20 ＋ ＋ ＋ 38 ⓪ ② ④ ⑥ ⓪ ④ ⓪
172. 例題：合計値の計算 変更後 5 8 3 2 9 1 6 4 14 9 9 6 23 15 38 ⓪ ② ⓪ ① ⓪ ① ③
173. 例題：合計値の計算 更に最適化していく 半数のスレッドは 1 回も reduction に参加しない ブロックあたりのスレッド数を半分にする 代わりに、 1 スレッドで 2 個の和を計算する ブロックあたりのスレッド数を固定すると、 ループを展開できる ループの条件分岐が不要になる
174. 例題：合計値の計算 int i = threadIdx.x + blockIdx.x * threadsPerBlock; __shared__ unsigned long long tmp[halfOfThreadsPerBlock]; tmp[threadIdx.x] = (i + 1) + (i + halfOfThreadsPerBlock + 1); __syncthreads(); if (threadIdx.x < 64) { tmp[threadIdx.x] += tmp[threadIdx.x + 64]; } 続く
175. 例題：合計値の計算 続き __syncthreads(); if (threadIdx.x < 32) { tmp[threadIdx.x] += tmp[threadIdx.x + 32]; __syncwarp(); tmp[threadIdx.x] += tmp[threadIdx.x + 16]; __syncwarp(); tmp[threadIdx.x] += tmp[threadIdx.x + 8]; __syncwarp(); tmp[threadIdx.x] += tmp[threadIdx.x + 4]; __syncwarp(); tmp[threadIdx.x] += tmp[threadIdx.x + 2]; __syncwarp(); tmp[threadIdx.x] += tmp[threadIdx.x + 1]; } if (threadIdx.x == 0) { atomicAdd(result, tmp[0]); }
176. 例題：合計値の計算 __syncwarp 関数 ワープ内でだけ同期を取る その分 __syncthreads より低コスト
177. 例題：合計値の計算 20 21 22 23 24 25 26 27 28 29 30 0.1 1 10 100 1000 GPU(naive) CPU(naive) GPU(reduction) GPU(bank conflict free) GPU(2 number/thread) GPU(loop unroll) M 実行時間 (ms)
178. 例題：合計値の計算 更に数倍速くなった
179. アトミック命令 異なるブロック同士で同期を取る手段の１つ 他には Co-operative group アトミック命令を使うと mutex が実装できる 基本的な演算でなくても任意の演算がほかのスレッドに 邪魔されずに実行できる 使いすぎると並列性を阻害するので注意
180. プロファイラを使おう プロファイリングをすることで、どこが ボトルネックになっているか見極め手助けになる デバイス側までプロファイリングするときは、 nvprof または NVIDIA Visual Profler が使える プロファイラでは各カーネルの実行時間や、 メモリ転送効率、 Occupancy などの情報を 収集できる プロファイラの使い方等はこの講座の範囲を超える ので省略
181. プロファイラを使おう
182. まとめ 本当に GPU による高速化が必要か考えよう GPU でプログラムを高速化するには 並列化する CPU と同様の高速化テクニック GPU 特有の高速化テクニック を組み合わせる プロファイラを使おう
183. 参考文献 [0] CUDA Toolkit Documentation http://docs.nvidia.com/cuda/index.html [1] CUDA C Programming Guide http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html [2] Better Performance at Lower Occupancy http://www.nvidia.com/content/GTC-2010/pdfs/2238_GTC2010.pdf
184. 参考文献 [3] CUDA By Example 汎用 GPU プログラミング入門 インプレスジャパン ISBN: 978-4-8443-2978-7 [4] CUDA C プロフェッショナルプログラミング インプレス ISBN: 978-4-8443-3891-8 [5] はじめての CUDA プログラミング 工学社 ISBN: 978-4-7775-1477-9