CUDAの高速化の復習2023年版 Reduction編
今回は最近のCUDA Samplesのコードを参考にCUDAでreductionを速くするテクニックのまとめになります。
私はCUDAを2009年のころから研究で使っていました。当時は頑張って勉強していたので自分の研究分野以外のCUDA応用の論文などを読んで高速化テクニックを勉強していました。ただ、2015年に企業に就職して以降、GPUを使うことはあってもCUDAのコードを直接自分で書くことはほとんどなくなってしまいました。このため、最近CUDAのコードを速くするためにどうすればいいのか？みたいな議論のときに、「昔は～」みたいな老害なコメントしかできない状態になってしまっていました。
この状態はさすがにまずいということでCUDAの勉強をし直そうと思い、この記事はそのまとめの第一弾でreductionをテーマで勉強した内容のまとめになります。私のように昔（2009年ごろ）CUDAを勉強したけど最近のCUDAわからんって人向けに記事は書いています。
今回測定のために書いたコードはこちらです。
https://github.com/shu65/cuda_reduction
この記事にはほとんど登場しないですが、昔風に書いたコードも一緒に含めています。
また元にしたNVIDIAのCUDA Samplesはこちらです。
https://github.com/NVIDIA/cuda-samples/tree/v11.8/Samples/2_Concepts_and_Techniques/reduction
計算時間に関してはCUDA Versionは12.0、GPUはV100を使った結果になります。計測はwarmupとしてreductionを実行したあと100回測定して平均計算時間を示しています。reductionの要素数は32M個に固定して測定します。また、reductionの計算では最初の段階はGPUを使いますが、CUDAの各thread blockの結果をさらにreductionする段階では簡単のためにCPUを使うことにします。
目次
Reductionとは？
Reductionとは配列を入力としてとり、配列の全要素の和などを計算する処理になります。C++のコードを見るのが一番わかりやすいと思うので、以下にC++のコードを示します。
int reduce(const int *in, size_t n)
{
int ret = 0;
for (size_t i = 0; i < n; ++i)
{
ret += in[i];
}
return ret;
}「和など」と書きましたが、各要素に対して使える演算は加算以外にも積、max, minなどいくつかあります。今回は説明を簡単にするために、加算で説明します。
そもそもなんでreductionをテーマにしようとしたかというと、私が勉強した当時、CUDAの高速化のテクニックがいろいろ詰まった題材としてreductionがよく紹介されていて、私もreductionを写経して勉強した経験がありました。このため、昔と現在の違いが分かりやすいということで選びました。
ちなみに昔の私が読んでたNVIDIAのreductionの資料はまだ公開されているようです。
https://developer.download.nvidia.com/assets/cuda/files/reduction.pdf
現在のCUDA Samplesで公開されているコードも基本的にはこの高速化テクニックに沿って実装されているようなので、この記事でも同じように沿って説明します。
Reductionを高速化していく
Reduce1: Baseline
まずはできるだけシンプルなreducetionのコードを最近のCUDA Sampleのコードを参考にしながら書いたコードを示します。
__global__ void reduce_gpu_v1_kernel(const int *g_in, size_t n, int *g_out)
{
cg::thread_block cta = cg::this_thread_block();
extern __shared__ int sdata[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
{
sdata[tid] = g_in[i];
}
else
{
sdata[tid] = 0;
}
cg::sync(cta);
for (int s = 1; s < blockDim.x; s *= 2)
{
if (tid % (2 * s) == 0)
{
sdata[tid] += sdata[tid + s];
}
cg::sync(cta);
}
if (tid == 0)
{
g_out[blockIdx.x] = sdata[0];
}
}このコードを見てblock内のスレッドの同期をする関数として__syncthreads() を使うんじゃないの？と思った方、いますよね？ちなみに私は思いました。昔の資料を参考にすると昔のreductionのコードは以下のような感じでした。
__global__ void reduce_gpu_old_v1_kernel(const int *g_in, size_t n, int *g_out)
{
extern __shared__ int sdata[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n)
{
sdata[tid] = g_in[i];
}
else
{
sdata[tid] = 0;
}
__syncthreads();
for (int s = 1; s < blockDim.x; s *= 2)
{
if (tid % (2 * s) == 0)
{
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid == 0)
{
g_out[blockIdx.x] = sdata[0];
}
}ほとんど同じですが、同期の部分で昔のコードでは__syncthreads() を使っているのに対して最近のコードではcooperative_groups:sync() を使うようになっています。
最近のCUDAを追ってない方はCooperative Groupsって何？と思った方もいると思うので簡単に説明します。Cooperative GroupsはCUDA 9から導入されたもので、様々な単位でスレッドの同期などを行うための仕組みになります。
CUDAでは通常スレッドのまとまりとしてblock、gridなどの単位があります。ただ、同期ができるスレッドの単位はこれまでblockくらいで、CUDAのカーネル関数内で他の単位で同期するのは結構面倒でした。
この問題を解決して、様々な単位、例えばblockよりも少ないスレッド数や、grid単位で同期したりできる仕組みがCooperative Groupsです。より詳しく知りたい方はこれらの資料に詳しく書かれていますのでご覧ください。
- https://developer.nvidia.com/blog/cooperative-groups/
- https://www.nvidia.com/content/apac/gtc/ja/pdf/2017/1041.pdf
このreduce1のほうの平均時間は以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
これをベースにして他の改良をしたらどうなるか？を示していきます。
Reduce2: Branch Divergenceの削減
Reduce1のコードの問題点の1つ目として以下のif文で実行されるスレッドが飛び飛びになってしまっているという問題点があります。
for (int s = 1; s < blockDim.x; s *= 2)
{
if (tid % (2 * s) == 0)
{
sdata[tid] += sdata[tid + s];
}
cg::sync(cta);
}この結果、1つのwarpでバラバラの処理が実行されることになり、branch divergenceが発生することになります。これはCUDAのコードの高速化をする際には注意する点の一つとなっています。詳しく知りたい方はこちらをご覧ください。
https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/#branching-and-divergence
これを解決するためにblock内のスレッドの0番から連続したスレッドがif文の中を実行するようにします。コードとしては以下の通りです。
for (uint32_t s = 1; s < blockDim.x; s *= 2)
{
uint32_t index = 2 * s * tid;
if (index < blockDim.x)
{
sdata[index] += sdata[index + s];
}
cg::sync(cta);
}コード全体としては以下のようになります。
https://github.com/shu65/cuda_reduction/blob/main/src/reduction_gpu_old.cu#L90-L119
このwarpのdivergenceの削減は昔から重要な高速化ポイントの一つで、昔の資料でも2.33倍高速化すると書かれています。ではこれを今回の環境で測定すると以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
ご覧の通り、現在でもwap divergenceはちゃんと削減すると効果があることがわかりました。
次はshared memoryのbank conflictの削減のための工夫です。shared memoryは高速アクセスできるのですがbank confictに注意する必要があります。memory bankは連続アドレスに割り当てられ、複数のスレッドが同じbankを使う場合は遅延が発生します。V100の場合、以下の資料にいくつか例でどういうときにbank conflictが起きるかが書かれているので、詳しく知りたい方は参考にしてください。
https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#shared-memory-5-x
注意点として昔のGPUとはmemory bankが16でしたが、最近のGPUはmemory bankが32らしいので注意してください。
bank conflictを削減するためには、飛び飛びのアドレスにアクセスするのではなく、各スレッドが連続したアドレスにアクセスするように改良します。先ほどしめしましたが、元々のコードは以下の通りです。
for (uint32_t s = 1; s < blockDim.x; s *= 2)
{
uint32_t index = 2 * s * tid;
if (index < blockDim.x)
{
sdata[index] += sdata[index + s];
}
cg::sync(cta);
}このコードでは例えば0番と1番のスレッドは5回目のイテレーションで32個隣をみるようになるので、このタイミングでmemory bankが衝突するようになります。
このコードに対してshared memoryのbank confilict削減するために以下のようにforループを工夫します。
for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1)
{
if (tid < s)
{
sdata[tid] += sdata[tid + s];
}
cg::sync(cta);
}こうすることでmemory bankが衝突することはなくなります。本当にbank conflictが削減しているのか？を確認したい場合、Nsight Computeというプロファイラを使うと確認できます。ちなみに昔からあったプロファイラのnvprofとnvvpはAmpereからサポートされなくなったので注意してください。
このコードを測定すると結果は以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
reduce3では1.2倍なので若干速くなっていますが、昔の資料をみると2倍速くなってたらしいのでだいぶ効果がうすれたなかという印象があります。
Reduce4: スレッドの実行効率向上
次にスレッドの実行効率の向上を図ります。
for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1)
{
if (tid < s)
{
sdata[tid] += sdata[tid + s];
}
cg::sync(cta);
}このループはよく見ると最初のイテレーションでblockの半分のスレッドはif文に入らないことが分かります。この結果半分のスレッドは一度もreductionの加算を実行しないことになります。これではせっかくスレッドを立ち上げたのにもったいないことになります。このため、効率をもう少しあげるために、以下のように最初にshared memoryに代入する部分を工夫します。
int tid = threadIdx.x;
int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
int sum_value = 0;
if (i < n)
{
sum_value = g_in[i];
}
if ((i + blockDim.x) < n)
{
sum_value += g_in[i + blockDim.x];
}
sdata[tid] = sum_value;このように最初shared memoryに足す前に各スレッドが2か所データを読み込んで加算してshared memoryに足すようにします。これで少なくとも1回は各スレッドがreductionの加算を実行することになります。またこの結果、カーネル実行時に起動するblock数を半分にできます。この結果計算時間は以下のようになります。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
| reduce4 | 0.250 | 1.814 | 3.876 |
結果からわかる通りやってみると1.8倍速くなっていて、昔も1.7倍の高速化があったらしいので、この改良は今でも効果的なことがわかります。
Reduce6: 完全なloop unroll
さて、次は1つ飛んでReduce6のloop unrollについて説明します。
昔のReduce5はwarp周りの改良なのですが、このwarp周りの事情が昔と今で変わっているのと、同じ部分に別の最適化の話が新しくCUDA Sampleに追加されているので最後にまとめて説明します。
では、reduce6の完全なloop unrollについてです。loop unrollはCUDAをやっていれば最後の手段的に出てくるテクニックで、ループの回数が分かっているならfor文やwhile文を使わずに直書きするというテクニックです。こうすることで、ループを抜けるかの条件判定などをなくすことができます。この結果、速くなるというものです。
今回のreductionでいえば、blockのスレッド数に依存してfor文のループ回数が決まっているので、blockのスレッド数を決め打ちしてあげればfor文なしで書けることになります。この時templateを使えば条件分岐で指定のblockのスレッド数の関数を呼び出すということもでき、ある程度柔軟にスレッド数を指定できるカーネル関数もできます。reductionの加算部分のfor文をunrollするとコードになります(長いので一部省略してます)。
if (kBlockSize >= 512)
{
if (tid < 256)
{
sdata[tid] += sdata[tid + 256];
}
cg::sync(cta);
}
if (kBlockSize >= 256)
{
if (tid < 128)
{
sdata[tid] += sdata[tid + 128];
}
cg::sync(cta);
}
...
if (kBlockSize >= 2)
{
if (tid < 1)
{
sdata[tid] += sdata[tid + 1];
}
cg::sync(cta);
}今回私の書いたコードは簡単のためにカーネル関数の呼び出しもとでblock数を512で決め打ちで書いています。
https://github.com/shu65/cuda_reduction/blob/main/src/reduction_gpu.cu#L284-L285
templateを使ったスレッド数の分岐に関して詳しく知りたい場合はCUDAのSampleの以下のコードをご覧ください。
この変更を加えた実行時間は以下のようになります。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
| reduce4 | 0.250 | 1.814 | 3.876 |
| reduce6 | 0.239 | 1.043 | 4.043 |
この変更も昔は1.41倍速くなっていたらしいのですが、ほとんど効果がなくなっているような印象です。個人的にはloop をunrollするとコードのメンテナンス性が非常に悪くなって嫌いなので、これの効果が小さくなっていることは私としてはちょっとうれしいです。
Reduce7: 1スレッドあたりの仕事を増やす
Reduce7では、昔の資料の最後の最適化で起動するスレッド数を減らしつつ、1スレッドあたりのreductionの加算の回数を増やすということをします。これを実現するためにreductionの最初はシーケンシャルに加算を実行していき、その後、いままでの加算する担当のスレッドを半分ずつ減らすというアルゴリズムにします。まず変更前のコードは以下のとおりです。
int tid = threadIdx.x;
int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
int sum_value = 0;
if (i < n)
{
sum_value = g_in[i];
}
if ((i + blockDim.x) < n)
{
sum_value += g_in[i + blockDim.x];
}
sdata[tid] = sum_value;
cg::sync(cta);この部分、以下のようにwhileループを追加します。
int tid = threadIdx.x;
int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
int grid_size = 2 * blockDim.x * gridDim.x;
int sum_value = 0;
while (i < n)
{
sum_value += g_in[i];
if ((i + kBlockSize) < n)
{
sum_value += g_in[i + kBlockSize];
}
i += grid_size;
}
sdata[tid] = sum_value;
cg::sync(cta);これに加えてカーネル関数の起動時のblock数を減らします。今回は決め打ちで128に固定しています。結果は以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
| reduce4 | 0.250 | 1.814 | 3.876 |
| reduce6 | 0.239 | 1.043 | 4.043 |
| reduce7 | 0.214 | 1.119 | 4.523 |
この最適化に関しては昔も1.42倍と効果が小さかったですが、さらに効果が小さくなっている印象です。
Reduce8 : スレッド数がwarpサイズになった以降の同期の最適化(Reduce5の部分)
昔の資料ではreduce7までの最適化ですが、ここからCUDAの新しい機能による最適化に関して紹介していきます。まずはスレッド数がwarpサイズになったところからの同期に関してです。昔、勉強してた人は「え、reduce5がそれでは？」と思った方もいるかもしれませんが、reduce5のコードそのままだと実は危険なので、改良したのものを紹介します。
昔、warpサイズは32で、連続したwarpサイズ分のスレッドは同時に実行されるので同期が必要ないという話がありました。これに基づいて同期を削除して高速化したのがreduce5です。ただ、CUDA 9からwarpに関するいろいろな機能追加があり、その中でwarpサイズ分のスレッドがすべて同時に実行されるという保証がなくなってしまいました。このためwarpサイズ分以下のスレッド数しかないとはいえ、同期なしだと何等かのエラーが発生する可能性があります。
このため、reduce5に同期を加える形にします。コードとしては以下の通りです(長いので一部省略してます)。
cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta);
if (cta.thread_rank() < 32)
{
if (kBlockSize >= 64)
{
if (tid < 32)
{
sdata[tid] += sdata[tid + 32];
}
tile32.sync();
}
if (kBlockSize >= 32)
{
if (tid < 16)
{
sdata[tid] += sdata[tid + 16];
}
tile32.sync();
}
...
if (kBlockSize >= 2)
{
if (tid < 1)
{
sdata[tid] += sdata[tid + 1];
}
tile32.sync();
}
}こちらのコードではまず、warpサイズの32でスレッド数が十分になったタイミングでtiled_partitionを使って32区切りのthread_block_tile を作り、そのうち0番から31番までのスレッドだけその後の処理をするようにします。さらに同期はこのタイル単位で行うことでwarpサイズ内での同期にしています。thread_block_tile を使わずに__syncwarp() を呼んで同期するという手もありますが、速度的にはそこまで変化しなかったので、ここではthread_block_tile を使うバージョンだけ測定します。
測定結果は以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
| reduce4 | 0.250 | 1.814 | 3.876 |
| reduce6 | 0.239 | 1.043 | 4.043 |
| reduce7 | 0.214 | 1.119 | 4.523 |
| reduce8 | 0.214 | 0.999 | 4.519 |
やってみるとこれに関しては全然速くならないという結果でした。ただ、ループのunrollはreduce6ですでにやってあるし、同期なしにはできなくなってしまっているのでこんなものかもしれないと思っています。
Reduce9 : shfl_down()の使用
最後にshfl_down() の使用による高速化です。先ほど紹介した通りCUDA 9のタイミングでwarp levelの関数がいくつか追加されました。warpレベルの同期もその一つですが、それに加えてデータ交換の関数もいくつか追加されています。これらの関数はCooperative Groupsにも同様の関数があり、今回はそちらを使うようにします。reduce8の段階のコードは以下の通りです。
cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta);
if (cta.thread_rank() < 32)
{
if (kBlockSize >= 64)
{
if (tid < 32)
{
sdata[tid] += sdata[tid + 32];
}
tile32.sync();
}
if (kBlockSize >= 32)
{
if (tid < 16)
{
sdata[tid] += sdata[tid + 16];
}
tile32.sync();
}
...
if (kBlockSize >= 2)
{
if (tid < 1)
{
sdata[tid] += sdata[tid + 1];
}
tile32.sync();
}
}shfl_down() を使う場合は以下のようになります。
cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta);
if (cta.thread_rank() < 32)
{
sum_value = sdata[tid];
if (kBlockSize >= 64)
{
sum_value += sdata[tid + 32];
}
for (int offset = tile32.size() / 2; offset > 0; offset /= 2)
{
sum_value += tile32.shfl_down(sum_value, offset);
}
}このときCUDA Sampleのコードに合わせてループのunrollをやめているので注意してください。またshfl_down()を使うケースではshared memoryを使わなくなっているので同期が不要になります。それ以外は基本処理の流れは同じです。
このコードの実行結果は以下の通りです。
| 平均計算時間(msec.) | 高速化率 | トータルの高速化率 | |
| reduce1 | 0.968 | 1.000 | 1.000 |
| reduce2 | 0.564 | 1.716 | 1.716 |
| reduce3 | 0.453 | 1.246 | 2.137 |
| reduce4 | 0.250 | 1.814 | 3.876 |
| reduce6 | 0.239 | 1.043 | 4.043 |
| reduce7 | 0.214 | 1.119 | 4.523 |
| reduce8 | 0.214 | 0.999 | 4.519 |
| reduce9 | 0.213 | 1.008 | 4.553 |
結果をみると少しだけ速くなっていることがわかります。
終わりに
今回、改めてCUDAの復習ということでreductionについての今時のコードの解説を書きました。やってみると結構昔と結果が変わって驚きました。なので、ちゃんと改めて復習してよかったと思っています。とりあえず、これで少しは老害発言をへらせるんじゃないかと思っています。
CUDA Sampleには今回紹介しなかった他のバージョンのreductionも書かれているのでもし気になる方は見てみてください。
