何がしたいのか
int sum(int *array, size_t n) {
int i, total = 0;
for (i=0; i<n; i++) total += a[i];
return total;
}
に相当する処理をGPUで書きたい、という話です。主にCUDAのサンプルの、6_Advanced/reductionに沿って書いています。
またこれや、その日本語版らしきやつも参考にしました。これらは、だいたい、サンプルと同じですが、途中から変わってきます。
(サンプルの方がバージョンが少し新しく、PDFの方が古いように見受けられます。)
がんばって調べて理解しながら書いたつもりですが。まだまだ甘い部分があると思いますので、強い人からのマサカリを歓迎します。
GPUでの計算
GPUで処理を行う場合は、CPUではありえないような、膨大なスレッド数を使ってプログラムを書くことになります。
例えば、2つのベクトルの加算、つまりCPUでいうところの
void add(int *a, int *b, size_t n) {
int i;
for (i=0; i<n; i++) a[i] += b[i];
}
は、GPUで処理する場合、CUDAでは
__global__ void add(int *a, int *b, size_t n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) a[i] += b[i];
}
になります。
これは、配列の要素数だけ(あるいは、場合によってはそれより少し多く)スレッドを立ち上げて、各スレッドは1要素だけを足し算する、という処理をしています。
配列の総和で考慮すべきこと
CUDAのサンプルで、配列の総和は 6_Advanced のフォルダに入っています。たかが総和が、Advancedとあるように、意外と難しいわけです。
ベクトル和のときと同じように、要素数だけスレッドを立ち上げるとしましょう。1つの変数に全てのスレッドが同時に足し算をしにいくことは、GPUでもできません。
(atomicAddはありますが。総和をそれでやったら、多分、とっても遅いです)
また、スレッドはブロックごとに分けられ、1ブロックあたり(GPUによってはもっと少ないが、最近のGPUだと)1024スレッドまでしか扱えません。スレッド間の同期は基本的にはブロック単位となります。
その都合上、一度のカーネル呼び出しではなく、各ブロックごとに集計する処理を何度か呼び出すことで、総和を行っています。
ここでは、サンプルの流れにそって、とりあえずやってみたreduce0から、いろんな観点からの効率化を図り、reduce6まで見てみます。
コピーライト表記
本記事では、CUDAのサンプルコードを引用しています。サンプルコードのライセンスについては、以下を参照下さい。
// This software contains source code provided by NVIDIA Corporation.
// http://docs.nvidia.com/cuda/eula/index.html#nvidia-cuda-samples-end-user-license-agreement
This article contains CUDA's sample source code or derivative of the code provided by NVIDIA Corporation.
reduce0: とりあえずやってみる
難しいことを考えず、まずやってみましょう。
コメントを読めば、だいたい何をやってるか分かると思います。
__syncthreads()で、ブロック内のスレッド間の同期をとれます。
__syncthreads()はif文の中には入れないでください。ループ内で使う場合は、ブロック内の全スレッドが、同じ回数だけループが回るように気をつけてください。
同期を取るということは、全てのスレッドが__syncthreads()に辿り着くのを待つということです。辿り着かないスレッドがあれば、いつまで経っても同期できません。
__global__ void reduce0(int *g_idata, int *g_odata, unsigned int n)
{
// shared memoryは、うまく使えば速度が速い。(ここでは、うまく使ってない。後で直す)
// 同じブロック内のスレッド間で共有される。
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
sdata[tid] = (i < n) ? g_idata[i] : 0;
__syncthreads();
// 例えばtidが0..7の場合、
// sdata[0] += sdata[1]; sdata[2] += sdata[3]; sdata[4] += sdata[5]; sdata[6] += sdata[7];
// sdata[0] += sdata[2]; sdata[4] += sdata[6];
// sdata[0] += sdata[4];
// のように、sdata[0]に集められる。
for (unsigned int s=1; s<blockDim.x; s*=2) {
if (tid % (2*s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
reduce1: divergent branchを減らそう
もう少し、CUDAでのスレッドの実行について話をしましょう。
CUDAでは「ワープ」という単位でスレッドが実行され、1ワープ = 32スレッドです。
ワープ内の各スレッドは、プログラムカウンタを共有しており、全スレッドが同じタイミングで動きます。では、条件分布があった場合にはどうなるでしょうか。
例えばif文があったとしてます。ひとつでも、if文の条件に合致するスレッドがあれば、そのスレッドがif文の中を実行しなければならず、その間、他のスレッドは待つことになります。また、elseがあった場合も、やはり同様です。
このような状況を、divergent branchと呼んでいます。
ワープ内は集団行動です。32人のうち、1人でもトイレに行きたい人がいれば、みんなが足止めを喰らいます。ならば、トイレに頻繁に行く人の多いワープと、トイレにあまり行かない人の多いワープに分けた方が高速化に寄与しますね。
ということで、こんな感じにできました。
// reduce0 では、例えばs=1のとき、
// tid = 0: sdata[0] += sdata[1];
// tid = 1: 待機
// tid = 2: sdata[2] += sdata[3];
// tid = 3: 待機
// ...
// のようになっていたが、今回は、
// tid = 0: sdata[0] += sdata[1];
// tid = 1: sdata[2] += sdata[3];
// ...
// tid = blockDim.x / 2 - 1: sdata[blockDim.x - 2] += sdata[blockDim.x - 1];
// tid = blockDim.x / 2: 待機
// tid = blockDim.x / 2 + 1: 待機
// ...
// のようになる。
__global__ void reduce1(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
sdata[tid] = (i < n) ? g_idata[i] : 0;
__syncthreads();
for (unsigned int s=1; s<blockDim.x; s*=2) {
int index = 2 * s * tid;
if (index < blockDim.x) {
sdata[index] += sdata[index + s];
}
__syncthreads();
}
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
reduce2: バンクコンフリクトを解消しよう
sharedメモリは、バンクコンフリクトを起こさないようにアクセスしたら、非常に高速です。
sharedメモリにアクセスする際、内部的には「バンク」に対して、欲しいメモリ番地をリクエストしています。メモリ番地ごとに、どのバンクにリクエストを出すかが決まっていて、ワープ内の複数スレッドが、同時に、同じバンクに対して、異なるアドレスでアクセスしたとき、バンクの順番待ちが起こります。これをバンクコンフリクトと呼んでいます。
ちなみに、同じバンクに同時にアクセスしていても、同じアドレスに対するアクセスであれば、バンクコンフリクトは起こりません。
(※同じアドレス/異なるアドレス、というのは厳密な表現ではなく、NVIDIAのC Programming Guideでは、32-bit wordという表現が使われています。つまり、アドレスが違っていても、同じ32-bit wordに収まる範囲であれば、バンクコンフリクトは起こらないと考えられます。)
さて。バンクコンフリクトを起こさないためには、一番簡単なのは、シーケンシャルにアクセスすることです。(バンクが32-bit wordなので、64-bitのdouble型などは、それでもバンクコンフリクトするような気がしますが、今回は32-bit整数を使うので、考えないことにします)
__global__ void reduce2(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
sdata[tid] = (i < n) ? g_idata[i] : 0;
__syncthreads();
for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
reduce3: 何も足し算していないスレッドに、1回くらいは足し算してもらう
なんとなんと。今までのコード、半分のスレッドは、一度も足し算しないまま終わってたんです。
半分ずつに集約していく形なので、仕事がないスレッドが多くなるのは仕方ないですが、さすがに、一度くらいは足し算してもらいましょう。
これに伴い、indexのとり方が変わります。また、呼び出し元でも、必要なブロック数が変わるので注意してください。
__global__ void reduce3(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
int mysum = (i < n) ? g_idata[i] : 0;
if (i + blockDim.x < n) mysum += g_idata[i + blockDim.x];
sdata[tid] = mysum;
__syncthreads();
for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
reduce4: [Kepler以降]Shuffle operationを使って高速に
Shuffleとは、同じワープ内の別のスレッドのローカル変数にアクセスする仕組みで、Keplerから導入されました。
ここには、Shuffleを使う利点として、sharedメモリを介したものだと複数回の操作が必要になるものを一度の操作で行えるために高速化(メモリ帯域幅の有効利用と、レイテンシの低減)が望めること、(今回はあまり気にしていないですが)容量が少なく貴重なsharedメモリを使わなくていいことが挙げられています。また、同じワープ内にしか作用しないため、__syncthreads()する必要がないことも挙げられています。(これについては、sharedメモリを使った場合も、ワープごとに依存関係がないように書けば、同様だと思いますが)
そのため、これ以降はKepler以降のカードじゃないと動きません。
また、コンパイル時に、 -arch=sm_30 オプションなどでKepler以降のカードをターゲットにする必要があります。
(サンプルのMakefileを見ればわかりますが、サンプルでは -gencode で、複数のアーキテクチャを指定しているようです)
なお、コンパイルオプションを指定せずに、Shuffle operationを使った場合、コンパイルエラーは
reduce.cu(172): error: identifier "__shfl_down" is undefined
のようになります。一見、ヘッダのインクルード忘れを疑いますが、そうではありません。
サンプルでは、shuffleがなければループアンローリングでもやっとけ、みたいなことが書いてあって、条件コンパイルで、shuffleがない場合はループアンローリングしたバージョンが動くようになっています。
私のコードでは、そこは面倒なのでやっていません。
(#error で処理するようにしたのですが、コンパイルオプションを付けてもエラーが出ました。
原因は恐らく、コードがホスト用のコードとデバイス用のコードに分けてコンパイルされる際、ホスト側では __CUDA_arch__ が定義されておらず、その状態でプリプロセッサ命令の解釈が行われたからです。回避のため、条件コンパイルで defined(__CUDA_arch__) も見るようにしましたが、もっといい方法あるはず・・・・)
というわけで、コードです。mysumに部分和を入れて、後で __shfl_down によって、ワープ内の別スレッドのmysumをもらってきています。
__global__ void reduce4(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
int mysum = (i < n) ? g_idata[i] : 0;
if (i + blockDim.x < n) mysum += g_idata[i + blockDim.x];
sdata[tid] = mysum;
__syncthreads();
for (unsigned int s=blockDim.x/2; s>32; s>>=1) {
if (tid < s) {
sdata[tid] = mysum = mysum + sdata[tid + s];
}
__syncthreads();
}
if (tid < 32) {
if(blockDim.x >= 64) mysum += sdata[tid + 32];
for (int offset = 32/2; offset>0; offset>>=1) {
mysum += __shfl_down(mysum, offset);
}
}
if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}
reduce5: 最初のループをアンローリングしよう
C++のテンプレート機能を使って、スレッド数をコンパイル時に決めたら、
最初のループを効率よくアンローリングできる、という魂胆のようです。
BLOCKDIM をコンパイル時に決めないといけないので、呼び出す側のコードがとてもアレになります。
理由はよくわからないですが。私の環境では、かなり効果は微妙でした。
また、サンプルコードでは、ブロックあたりのスレッド数が1024にならない前提のコードでしたが、最近のアーキテクチャでは1024スレッドまで対応しています。
その場合、コードの修正が必要となりますので、ご注意ください。(以下コードでは、1024スレッドでも動くようにしています。たぶん今はないけど、2048スレッド以降はごめんなさい。)
template <unsigned int BLOCKDIM>
__global__ void reduce5(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * (BLOCKDIM * 2) + threadIdx.x;
int mysum = (i < n) ? g_idata[i] : 0;
if (i + BLOCKDIM < n) mysum += g_idata[i + BLOCKDIM];
sdata[tid] = mysum;
__syncthreads();
if (BLOCKDIM >= 1024 && tid < 512) sdata[tid] = mysum = mysum + sdata[tid + 512];
__syncthreads();
if (BLOCKDIM >= 512 && tid < 256) sdata[tid] = mysum = mysum + sdata[tid + 256];
__syncthreads();
if (BLOCKDIM >= 256 && tid < 128) sdata[tid] = mysum = mysum + sdata[tid + 128];
__syncthreads();
if (BLOCKDIM >= 128 && tid < 64) sdata[tid] = mysum = mysum + sdata[tid + 64];
__syncthreads();
if (tid < 32) {
if (BLOCKDIM >= 64) mysum += sdata[tid + 32];
for (int offset = 32/2; offset>0; offset>>=1) {
mysum += __shfl_down(mysum, offset);
}
}
if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}
reduce6: ブロック数を制限してみる
ここで突然、 Brent's theorem という言葉が出てきます。
これは、並列計算における計算時間を見積もる定理で、 $p$ 個のプロセッサがあった場合の計算時間 $T_p$ を、 $1$ 個のプロセッサでの計算時間 $T_1$ と、プロセッサが無尽蔵にあった場合の計算時間 $T_\infty$ で表すと次のようになります。
$T_p = T_\infty + \frac{T_1 - T_\infty}{p}$
で。さらに、計算コストという概念も考えましょう。計算コストとは、プロセッサごとの計算時間の和で、各プロセッサが同じ計算時間を計算したとみなせるなら、プロセッサ数と計算時間の積になります。計算コストはプロセッサ数と計算時間の積で表せるものとして、総和計算の計算コストを見積もりましょう。
まず、1プロセッサで総和を取ろうと思ったら、実行時間は明らかに $T_1 = O(N)$ です。また、プロセッサが無尽蔵にあった場合は、配列の要素数だけプロセッサを使って、要素数を半分に減らす処理を繰り返すことになるので $T_\infty = O(\log N)$ でいけます。そうすると、
$T_p = T_\infty + \frac{T_1 - T_\infty}{p} = O(\log N + \frac{N - \log N}{p})$
となり、計算コストは
$p T_p = O(p \log N + N)$
となります。なので、$p = N$のときは、計算コストは $O(N \log N)$ となり、1プロセッサでの計算コスト $O(N)$ よりも悪いんですね。
ところが、もしも $p = N / \log N$とすれば、計算コストは $O(N)$となり、1プロセッサの場合の計算量と一致します。
これをするには、各スレッドで、今までよりも多くの要素を足し算し、ブロック数を抑えることが重要です。
……という話をした後、ズコーなんですが。reduce6では、ブロック数の最大を64と置いて実行します。
Brent's theoremって言葉も、ちゃんとNVIDIAの資料にも、サンプルコードのコメント欄にも書いてあるんですが。一体なんだったのか。
あと、サンプルコードでは、nが2の累乗かどうかを静的に渡して、計算量を減らすって工夫があったんですが。それすると、結果が合わないような気がします。
私が写経ミスったのか、普通にバグってるのか知りたいんですが、サンプルコードのやつを適当にオプションつけて動かしても、答えが合わないことが結構あるんですよね。なんなんだろ。
template <unsigned int BLOCKDIM, bool N_IS_POW2>
__global__ void reduce6(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * (BLOCKDIM * 2) + threadIdx.x;
unsigned int gridsize = BLOCKDIM * 2 * gridDim.x;
int mysum = 0;
while (i < n) {
mysum += g_idata[i];
if (/*N_IS_POW2 ||*/ i + BLOCKDIM < n) mysum += g_idata[i + BLOCKDIM];
i += gridsize;
}
sdata[tid] = mysum;
__syncthreads();
if (BLOCKDIM >= 1024 && tid < 512) sdata[tid] = mysum = mysum + sdata[tid + 512];
__syncthreads();
if (BLOCKDIM >= 512 && tid < 256) sdata[tid] = mysum = mysum + sdata[tid + 256];
__syncthreads();
if (BLOCKDIM >= 256 && tid < 128) sdata[tid] = mysum = mysum + sdata[tid + 128];
__syncthreads();
if (BLOCKDIM >= 128 && tid < 64) sdata[tid] = mysum = mysum + sdata[tid + 64];
__syncthreads();
if (tid < 32) {
if (BLOCKDIM >= 64) mysum += sdata[tid + 32];
for (int offset = 32/2; offset>0; offset>>=1) {
mysum += __shfl_down(mysum, offset);
}
}
if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}
実行結果
複数回動かしたり平均取ったりはしていないです。
時間の単位はmsです。
N: 268436690
CPU: 630.4120
threads 32 64 128 256 512 1024
reduce0: 20.4983 20.9556 22.0257 23.6568 27.0491 31.7910
reduce1: 14.4885 11.6330 9.9901 11.8105 13.5425 16.4100
reduce2: 14.6104 11.7311 9.7764 10.9015 12.3254 14.1772
reduce3: 7.4528 6.1384 5.8372 5.7876 6.2441 7.2453
reduce4: 5.8222 4.9664 4.5815 4.7934 4.5182 5.2623
reduce5: 6.8359 5.0480 4.9372 4.6059 4.5440 4.8092
reduce6: 21.0829 10.3854 5.8892 5.4824 5.2538 5.0353
reduce5の効果がすごく疑問です。やらない方がよくね?
reduce6は、計算時間よりも計算コストを下げる、ということに焦点を置いているので、時間が遅くなるのは仕方ないのかもしれません。
アーキテクチャもNの数も違うのでそんなもんかもしれないですが、NVIDIAの資料での結果とは傾向が少し違うようには思います。
パフォーマンスの評価
以下の2つの観点から、パフォーマンスを追求することになります。
分からない数字は、グラボ屋さんのWebページか、NVIDIAのWebページか、CUDAのサンプルについている 1_Utilities/deviceQuery をコンパイルして、deviceQueryで確認できます。
計算速度: GFLOP/s
2 * CUDA Core数 * GPUクロック数 (GHz) … ただし、float型の場合
ここで出てくる2は、1命令で2個のfloatを扱える、という意味らしいですが、後述の理由から、今回はGFLOP/sは見ません。
メモリ帯域幅: GB/s
メモリのバス幅(bit) / 8 * メモリのクロック数 (GHz) * 2
ここで出てくる2は、積んでるメモリがDDR (Double-Data-Rate)という意味です。
Pascalアーキテクチャ(GeForce GTX 10xx)だとGDDR5X積んでてDDRじゃなくQDRで動く、という話も出てくるんですが、
NVIDIAの出してるスペックとdeviceQueryとを突き合わせて計算した感じ、2で正しいような気がします。
どちらも理論上の最大値なので、これに達することはないはずですが。どれだけ近づいたかを指標にできます。
また、総和計算では、計算の数は大したことないので、メモリがボトルネックになるはず、ということで、計算速度は見ないでメモリ帯域幅でパフォーマンスを測ります。
今回、メモリのバス幅256ビット、クロック数5505 MHzのGTX 1080を使用したので、おおよそ350 GB/sが理論値ということになります。
どのようにして実測値を求めるのかは、正直よくわからなかったのですが、
理論上の最高速が、入力した配列の各要素に1回ずつ読み込みアクセスすること、と考えて、
アクセスする容量(GB) = 要素数 * sizeof(int) / (1024*1024*1024)
実行時間(s) = 実行時間(ms) / 1000
メモリ帯域幅(GB/s) = アクセスする容量(GB) / 実行時間(s)
とすれば、だいたいNVIDIAの資料の値と近い数字は出ることが分かりました。
(NVIDIAの資料の値の方が、もう少し、いい値が出ている。多分、俺が何かを見落としているんでしょう)
試した中で、最も計算時間の短かった、reduce4の512スレッドでの時間 4.5182 ms で、メモリ帯域幅を求めてみると。
アクセスする容量 = 268436690 * sizeof(int) / (1024*1024*1024) = ほぼ1GB
実行時間 = 4.5182 ms = 0.0045182 s
メモリ帯域幅 = 約 221 GB/s
350 GB/sの理論値に対し、6割強は出ました。
もう少しがんばりたい気はしますが、とりあえずここまで。
コード全体と、雑に書きすぎてちょっと恥ずかしいhost側コードはこちらに置いています。
https://gist.github.com/gyu-don/7aa8c013e966579862323a764b28f794
最後に
マサカリ歓迎です。
みなさま、よきCUDAライフを。
参考資料
http://developer.download.nvidia.com/compute/cuda/1.1-Beta/x86_website/projects/reduction/doc/reduction.pdf
http://gpu-computing.gsic.titech.ac.jp/Japanese/Lecture/2010-06-28/reduction.pdf
https://stackoverflow.com/questions/15055877/how-to-get-memory-bandwidth-from-memory-clock-memory-speed
https://en.wikipedia.org/wiki/Pascal_(microarchitecture)#Performance
http://docs.nvidia.com/cuda/cuda-c-best-practices-guide/index.html#branching-and-divergence
https://devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/
https://everything2.com/title/Brent%2527s+theorem