iPX社員によるブログ

iPX社員が"社の動向"から"自身の知見や趣味"、"セミナーなどのおすすめ情報"に至るまで幅広い話題を投下していくブログ。社の雰囲気を感じ取っていただけたら幸いです。

CUDA Parallel Reduction

初めに

ご無沙汰しております。iPXの小川です。
弊社では定期的に社内でのCUDA勉強会を実施しております。

新入社員も増え、CUDA教育を行う機会も増えましたが、CUDAアーキテクチャの説明が文や図解だけだと難しく感じている今日このごろです。

そこで本日は私が理解を果たした過程で、一番寄与した Parallel Reduction を題材にサンプルコードを載せつつ寄り道をしながら解説してきたいと思います。

問題点の考察

CUDAの説明はNVIDIA様のオープンな資料・サンプルが充実し、「CUDA C プロフェッショナル プログラミング」などの良書も存在します。
ではなにゆえ壁に直面するのか、その問題点を実体験を元に考察してみます。

CUDA学習の中でよく直面したケースは、下記の黄金パターンでした。

  • 概論の習得(理解したつもり)
  • サンプルコードの読み解き
  • !理解不能!

何故なのか。

一つの原因は、資料中の説明がある側面から見てなされているからだ、と考えます。
CUDAは一言では並列コンピューティングプラットフォームですが、汎用的なプログラミングの為、デバイスの抽象化やメモリ管理の抽象化などあらゆる抽象化がなされています。
その全てをいきなり説明することも理解することも難しいので、段階によって側面に分解され説明がなされるのは順当な理由です。

よって自分が今学んでいる箇所はデバイスから見た側面なのか、プログラミングモデルから見た側面なのか、メモリモデルから見た側面なのか、説明の視点を理解することが迷子にならない最良の方法です。

CUDAの全てではありませんが、よく使う部分を抜粋して図解しようと思います。
記事で何をやっているかわからなくなった場合は図解を眺め直すことをおすすめします。

f:id:ipx-writer:20170831001430p:plain

f:id:ipx-writer:20170831001448p:plain

f:id:ipx-writer:20170831001503p:plain

※実際にはデバイスによってSMの数などは異なります。

それでは次項よりサンプルコードを載せて解説していきます。


CUDAプログラミングモデル

この記事をご覧になられている方は既にCUDAプログラミングモデルを習得済みかもしれません。
復習を兼ねてベクトル和のCUDAカーネルを実装します。
なお開発環境は Visual Studio 2015 ですが、Eclipse や Nsight でも同様に動作すると思います。

// CUDA
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
// Thrust
#include <thrust/execution_policy.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/copy.h>
// C++
#include <iostream>

template <class T>
__global__ void plus_kernel(const T* g_x, const T* g_y, T* g_o) {
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	g_o[i] = g_x[i] + g_y[i];
}

int main(void) {
	using namespace thrust;

	const int blocks = 1;
	const int threads = 32;
	const int N = blocks * threads;

	counting_iterator<int> ci(1);
	device_vector<int> d_x(ci, ci + N);
	device_vector<int> d_y(ci, ci + N);
	
	device_vector<int> d_o(blocks*threads);
	plus_kernel<int> <<< blocks, threads >>> (
		raw_pointer_cast(d_x.data()),
		raw_pointer_cast(d_y.data()),
		raw_pointer_cast(d_o.data())
		);

	copy(host, d_o.cbegin(), d_o.cend(), std::ostream_iterator<int>(std::cout, " "));

	return 0;
}

単純に配列の各要素を足した結果を返却しています。
ここでグリッド内に存在するブロックを1、ブロック内に存在するスレッドを32と定義しました。
論理コンポーネントのサイズはプログラマによって変更可能なのでした。


ここで Thrust についてご存じない方がいらっしゃるかもしれません。
今回は本質であるCUDAカーネルに注視したいので、Thrust にGPUメモリ確保・解放を担わせています。
簡単に概要を説明しておきます。

Thrust

STLライクなCUDAラッパーライブラリです。
CUDA Toolkit インストール時に一緒にインストールされているので、特別な環境構築は不要です。

キーワード* 意味*
execution_policy.h Thrustで実装されている関数の実行プロファイルをhost、deviceで明示的に指定出来ます。copy関数がSTLのcopyと衝突するので、hostを明示的に指定するのに使用しています。
device_vector バイスメモリの確保・解放を管理してくれます。dataメンバ関数で配列の先頭アドレスを取り出すことが出来ます。
counting_iterator 指定した定数からの連番を返すイテレータです。
raw_pointer_cast device_vectorから取り出したアドレスをCUDAカーネルに渡せるよう生ポインタへ変換します。

ご興味のある方は CUDA Toolkit Documentation を参照してください。

docs.nvidia.com

Reduction

私は日本語で畳み込みと読んでいますが、配列の次元を削減するような処理を指します。
これは配列の各要素に対してある計算を行い、最終的に一つの値にする様が配列を畳んでいるように見えるからです。
例として1次元配列ではよく合計やmax/minなどの計算を行いますが、これらは正に畳み込み処理です。

CPU実装

先ずはCPUで実装してみます、for文を使う場合下記になるでしょう。
STLの場合は accumulate関数が畳み込みの用途で使えます。
accumulateは所謂 map,reduce,filterパターンでいうことろのreduceです。
C++では transform,accumulate,copy_ifになります。

#include <iostream>
#include <numeric> // accumulateに必要

template <class T, int N>
T reduction_cpu(const T* x) {
	T sum{};
	for (int i = 0; i < N; ++i) {
		sum += x[i];
	}
	return sum;
}

int main(void) {
	const int N = 5;
	int x[N] = { 1,2,3,4,5 };
	std::cout << "reduction_cpu:" << reduction_cpu<int, N>(x) << std::endl;
	std::cout << "reduction_cpu_stl:" << std::accumulate(x, x + N, 0) << std::endl;
	return 0;
}

CUDA実装

それではGPUで実装してみます。
一つのアイデアはCUDAカーネルへ渡した引数に足し込むことです。
CUDAカーネルに渡された引数はグローバルメモリに配置される為、全てのスレッドから参照可能です。
しかしこれは失敗に終わります、実験してみます。

template <class T>
__global__ void reduction1_kernel(const T* g_x, T* g_o) {
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	g_o[0] += g_x[i];

int main(void) {
	using namespace thrust;

	const int blocks = 1;
	const int threads = 5;
	const int N = blocks * threads;

	counting_iterator<int> ci(1);
	device_vector<int> d_x(ci, ci + N);
	device_vector<int> d_o(1);

	reduction1_kernel<int> << <blocks, threads >> > (
		raw_pointer_cast(d_x.data()),
		raw_pointer_cast(d_o.data())
		);

	copy(host, d_o.cbegin(), d_o.cend(), std::ostream_iterator<int>(std::cout, " "));
	return 0;
}

1+2+3+4+5=15 を期待しましたが、結果は5でした。
これは各スレッドが並列に動作しているからだと暫定理解しましょう(後にウォープを解説します)。
よって畳み込むためにはスレッド間で値を交換する必要があります。

共有メモリ

冒頭の図解によるとブロック内において、各スレッドは共有メモリを参照することが可能でした。
ブロック内における各スレッドの足並みを揃えるには __syncthreads関数を用います。
各スレッドが足並みを揃える様がまるでバリアを張ったようである為、バリア同期と呼ばれます。

template <class T, int N>
__global__ void reduction2_kernel(const T* g_x, T* g_o) {
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	int tid = threadIdx.x;

	__shared__ T s_x[N];
	s_x[tid] = (i < N) ? g_x[i] : T{};
	__syncthreads();

	for (int s = 1; s < blockDim.x; s *= 2) {
		if ((tid % (2 * s)) == 0) {
			s_x[tid] += s_x[tid + s];
		}
		__syncthreads();
	}
	if (tid == 0)
		g_o[blockIdx.x] = s_x[0];
}

急に難しくなった!詐欺や!と思われるのも仕方がありません。
私も初めてみたときはそう感じました、複雑に見えますがfor文以外はとても単純です。
ひとつづつ分解します。

下記はブロック内におけるスレッド数分の共有メモリ領域を確保する処理です。
グローバルメモリ(g_x)から共有メモリ(s_x)へ値を転送します。
そのままだと先走るスレッドが存在するかもしれないので、バリア同期で足並みをそろえます。
これで共有メモリの初期化が完了です。

__shraed__ T s_x[N];
s_x[tid] = (i < N) ? g_x[i] : T{};
__syncthreads();


下記は共有メモリの0番目に合計された値が格納されているので、グローバルメモリ(g_o)へ合計値を返却する処理です。
図解より共有メモリはブロック内の各スレッドが共有できるメモリです。
この共有メモリ(s_x[0])に合計値がforループによって格納されることになります。
そしてスレッドID(tid)が0番目の場合は該当するブロックID(blockIdx.x)が示す領域へ返却します。
tid == 0 の条件が無くとも動作はしますが、無駄な処理となります。
仮に条件を抜いてしまうと、ブロックに存在する全スレッドが共有メモリ(s_x[0])へのアクセスを行うので、順番待ちが発生してしまいます。

if (tid == 0)
	g_o[blockIdx.x] = s_x[0];

最後にforループ部分です、少し複雑なので図解します。

for (int s = 1; s < blockDim.x; s *= 2) {
	if ((tid % (2 * s)) == 0) {
		s_x[tid] += s_x[tid + s];
	}
	__syncthreads();
}

f:id:ipx-writer:20170831025420p:plain

一度理解してしまえばとても理にかなっている事が分かります。
2,4,8,.. と飛ばして共有メモリを足し込むことによって、最終的に共有メモリ(s_x[0])には、ブロック内における全スレッドの合計値が格納されるという仕組みです。

これで共有メモリを使用した Parallel Reduction を行うことが出来ました。
しかしCUDAでは共有メモリよりも高速なレジスタが存在します。
次項ではレジスタとウォープの関係について説明します。

レジスタ

レジスタはCUDAの中でも最も高速なメモリ領域です。
このレジスタを積極的に活用する事が基本的な戦略となります。
併せて、CC(コンピュートケーパビリティ)が3.0以上の場合はシャッフル関数が使用可能です。
このシャッフル関数はレジスタ経由でスレッド間の値を交換することが出来ます。
ただし、シャッフル関数はウォープ内(32スレッドをまとめた論理コンポーネント)でのみ動作が可能です。

これはCUDAはSIMT(Single Instruction Multiple Thread)アーキテクチャを採用しており、32スレッドは並列で動作しますが、協調的に動作するスレッドによるデータ並列コードも記述出来る仕組みに起因しています。

重要なのはスレッドはレジスタが使える、ウォープはレジスタの値を交換出来るという点です。
前項で共有メモリを使用して記述した Parallel Reduction をシャッフル関数を使い高速化してみます。

※compute, sm の設定は忘れずに行ってください。

f:id:ipx-writer:20170831033848p:plain

template <class T>
__device__ T warpReduction(T val) {
	for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
		val += __shfl_down(val, offset);
	}
	return val;
}

template <class T>
__device__ T blockReduction(T val) {
	static __shared__ T s[32];
	int laneId = threadIdx.x & (32 - 1);
	int wid = threadIdx.x / 32;

	val = warpReduction<T>(val);
	if (laneId == 0)
		s[wid] = val;
	__syncthreads();

	val = (threadIdx.x < blockDim.x / 32) ? s[laneId] : 0;
	if (wid == 0)
		val = warpReduction(val);

	return val;
}

template <class T>
__global__ void reduction3_kernel(const T *g_x, T *g_o) {
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	int sum = blockReduction<T>(g_x[i]);
	if (threadIdx.x == 0)
		g_o[blockIdx.x] = sum;
}

int main(void) {
	using namespace thrust;
	const int blocks = 2;
	const int threads = 32;
	const int N = blocks * threads;

	counting_iterator<int> ci(1);
	device_vector<int> d_x(ci, ci + N);
	device_vector<int> d_o(blocks);

	reduction3_kernel<int><<<blocks, threads >>> (
		raw_pointer_cast(d_x.data()),
		raw_pointer_cast(d_o.data())
		);

	copy(host, d_o.cbegin(), d_o.cend(), std::ostream_iterator<int>(std::cout, " "));

	return 0;
}

今までの総まとめといった形ですが、__shfl_down について少し解説しておこうと思います。

offsetで指定しているのはスレッドの後半部分となりますです。
ウォープ内のスレッドは32スレッドなので、ループを展開するとoffsetが 16, 8, 4, 2, 1 と変化します。
レーンID(ウォープ内のインデックスの意味で0~31の値)、offset、__shfl_downで取得されるスレッドについてまとめます。

ループ一回目

レーンID offset __shfl_down(val,offset)
0 16 16レーンの__shfl_downに指定したval
1 16 17レーンの__shfl_downに指定したval
2 16 18レーンの__shfl_downに指定したval
3 16 19レーンの__shfl_downに指定したval
4 16 20レーンの__shfl_downに指定したval
5 16 21レーンの__shfl_downに指定したval
6 16 22レーンの__shfl_downに指定したval
7 16 23レーンの__shfl_downに指定したval
8 16 24レーンの__shfl_downに指定したval
9 16 25レーンの__shfl_downに指定したval
10 16 26レーンの__shfl_downに指定したval
11 16 27レーンの__shfl_downに指定したval
12 16 28レーンの__shfl_downに指定したval
13 16 29レーンの__shfl_downに指定したval
14 16 30レーンの__shfl_downに指定したval
15 16 31レーンの__shfl_downに指定したval

ループ二回目

レーンID offset __shfl_down(val,offset)
0 8 8レーンの__shfl_downに指定したval
1 8 9レーンの__shfl_downに指定したval
2 8 10レーンの__shfl_downに指定したval
3 8 11レーンの__shfl_downに指定したval
4 8 12レーンの__shfl_downに指定したval
5 8 13レーンの__shfl_downに指定したval
6 8 14レーンの__shfl_downに指定したval
7 8 15レーンの__shfl_downに指定したval

ループ三回目

レーンID offset __shfl_down(val,offset)
0 4 4レーンの__shfl_downに指定したval
1 4 5レーンの__shfl_downに指定したval
2 4 6レーンの__shfl_downに指定したval
3 4 7レーンの__shfl_downに指定したval

ループ四回目

レーンID offset __shfl_down(val,offset)
0 2 2レーンの__shfl_downに指定したval
1 2 3レーンの__shfl_downに指定したval

ループ五回目

レーンID offset __shfl_down(val,offset)
0 1 1レーンの__shfl_downに指定したval

終わりに

弊社ではCUDAだけではなく、Deep Learning や自動車制御開発、MBDコンサルティング等の事業を行っております。
次回のオートモーティブワールド2018では、それらをフュージョンした弊社ならではの方法論を提示できたらと考えております。

最後までお付き合いいただき、有難う御座いました。

参考文献