N x Mの配列に大量の値を加算する。

最近CUDAを勉強中です。今回はCUDAの高速化について簡単にテストしてみました。

条件

• 出力はWxHの配列
• 各要素に対して任意の値(今回は正弦波なので振幅、周波数、位相を入力)をN回加算する。

変数

• width: 配列のXサイズ
• height: 配列のYサイズ
• n: 加算する正弦波の数
• max_thread_size: 1ブロック内で同時に稼働できるスレッド数。今回は1024。
• d_mag: 振幅
• d_freq: 周波数
• d_phi: 位相
• d_tmap: 出力配列

①脳筋コーディング

widthブロック内にぞれぞれheightスレッド立てて、カーネルではn個の正弦波をループで加算する。

__global__ void cuSigGen0
(
	int width,
	int height,
	int n,
	float *d_mag,
	float *d_freq,
	float *d_phi,
	float *d_tmap
)
{
	int idx = blockDim.x * blockIdx.x + threadIdx.x;

	if (idx >= width*height) return;

	d_tmap[idx] = 0.0f;

	for (int i = 0; i < n; i++)
	{
		d_tmap[idx] += (d_mag[i] / (float)n) * sinf(2 * kPi*(fc + d_freq[i] * 1e2f)*(float)idx / (float)width + 2 * kPi*d_phi[i]);
	}
}

②-1 2次元ブロックでatomicAddを使う

CUDAではブロックを2次元に持つことができるため、width x heightのブロックを作成して、各ブロックにmax_thread_sizeスレッド立てる。カーネルではn個の正弦波をmax_thread_size個ずつ求めて配列に加算していく。このとき同じ配列要素に複数のスレッドが同時に書き込みを行う可能性がある。
CUDAでは複数のスレッドが同じメモリを書き込みに行かないようにatomic命令が用意されている。一般的には最後の手段であり、あまり使わないほうが良いとされる。

__global__ void cuSigGen1
(
	int width,
	int height,
	int n,
	int max_thread_size,
	float *d_mag,
	float *d_freq,
	float *d_phi,
	float *d_tmap
)
{
	int bidx = blockIdx.x;
	int bidy = blockIdx.y;
	int tidx = threadIdx.x;

	if (bidy >= height || bidx >= width) return;

	int bid = width * bidy + bidx;

	int done = 0;

	for (; done < n; done += max_thread_size)
	{
		if (tidx + done >= n) return;

		float tmp = (d_mag[tidx + done] / (float)n) * sinf(2 * kPi*(fc + d_freq[tidx + done] * 1e2f)*(float)bid / (float)width + 2 * kPi*d_phi[tidx + done]);
		atomicAdd(&d_tmap[bid], tmp);
	}
}

②-2 2次元ブロックでシェアードメモリを使ってreductionを行う

基本的な作戦は先と同じだが、atomicAddを避けるためにカーネル内にシェアードメモリを作成し、その中で加算を行ってから出力に書き込む。こうすることで同じ配列要素に複数スレッドが書き込むことがなくなる。

__device__ __forceinline__ void smSum
(
	float *sm,
	int tidx,
	int max_thread_size
)
{
	__syncthreads();

	for (int k = 1; k < max_thread_size; k *= 2)
	{
		sm[tidx] += (tidx % (2 * k) == 0) ? sm[tidx + k] : 0.0f;

		__syncthreads();
	}
}

__global__ void cuSigGen2
(
	int width,
	int height,
	int n,
	int max_thread_size,
	float *d_mag,
	float *d_freq,
	float *d_phi,
	float *d_tmap
)
{
	extern __shared__ float sm[];

	int bidx = blockIdx.x;
	int bidy = blockIdx.y;
	int tidx = threadIdx.x;

	if (bidy >= height || bidx >= width) return;

	int bid = width * bidy + bidx;

	for (int done = 0; done < n; done += max_thread_size)
	{
		//float tmp = (d_mag[tidx + done] / (float)n) * sinf(2 * kPi*(fc + d_freq[tidx + done] * 1e2f)*bid / (float)width + 2 * kPi*d_phi[tidx + done]);
		
		sm[tidx] = (tidx + done < n) ? (d_mag[tidx + done] / (float)n) * sinf(2 * kPi*(fc + d_freq[tidx + done] * 1e2f)*(float)bid / (float)width + 2 * kPi*d_phi[tidx + done]) : 0.0f;

		smSum(sm, tidx, max_thread_size);

		if (tidx == 0) d_tmap[bid] += sm[0];
	}
}

③ 1次元ブロックのカーネルをホスト側でループする。

CUDAの場合スレッドやブロックの立てすぎでも速度が低下するときがある。先の手法のブロックを1次元に減らし、ホスト側でカーネルをループする。

__global__ void cuSigGen3
(
	int width,
	int height,
	int n,
	int max_thread_size,
	int offset,
	float *d_mag,
	float *d_freq,
	float *d_phi,
	float *d_tmap
)
{
	extern __shared__ float sm[];

	int bidx = blockIdx.x;
	int tidx = threadIdx.x;

	int bid = bidx;

	int done = 0;

	for (; done < n; done += max_thread_size)
	{
		//float tmp = (d_mag[tidx + done] / (float)n) * sinf(2 * kPi*(fc + d_freq[tidx + done] * 1e2f)*(bid + offset) / (float)width + 2 * kPi*d_phi[tidx + done]);

		sm[tidx] = (tidx + done < n) ? (d_mag[tidx + done] / (float)n) * sinf(2 * kPi*(fc + d_freq[tidx + done] * 1e2f)*(float)(bid + offset) / (float)width + 2 * kPi*d_phi[tidx + done]) : 0.0f;

		smSum(sm, tidx, max_thread_size);

		if (tidx == 0) d_tmap[bid + offset] += sm[0];
	}
}

int step = 1 << 15;
grid = dim3(step, 1, 1);
block = dim3(max_thread_size, 1, 1);

for (int offset = 0; offset < width*height; offset += step)
	{
		cuSigGen3 << < grid, block, max_thread_size * sizeof(float) >> >
			(
				width,
				height,
				n,
				max_thread_size,
				offset,
				d_mag,
				d_freq,
				d_phi,
				d_tmap
				);

		cudaDeviceSynchronize();
	}

時間測定

使ったGPUはRTX2070で諸元は以下の通り。

• width = 1024
• height = 512
• n = 100～1000000
• 周波数などはランダム(シードは0固定)

結果(単位はms)

	100	1000	10000	100000	1000000
①	2.457216	19.241217	158.047623	1398.913086	13832.567383
②-1	56.279041	512.576782	3783.871826	36508.042969	-
②-2	113.077759	85.194466	714.860535	7208.444336	79598.742188
③	95.010529	89.894913	727.293640	7305.295410	86867.156250

脳筋コード速すぎる。。。試行錯誤はなんだったんだ。。。
あとatomicAddはみんなが言う通りあまり使わないほうがよさそうですね。
なお、立てるスレッド数などでも高速化が狙えるが、今回はそういった最適化はしていない。
脳筋コードではnに比例して増加するのでどうにかいいやり方はないものかな。。。あっ。