[CUDA + OpenMP] CPUの複数スレッドで同時カーネル実行

##はじめに
この記事は，CUDAとOpenMPを使用した並列計算に関して私がハマったことについての記録です．

OpenMPでCPU側に複数のスレッドが生成されていて，その各スレッドがそれぞれカーネルを実行するというときにちょっとした作法があるみたいです．

以下

- 動作環境
- CUDA Streamの簡単な説明
- OpenMPと一緒に使う方法の説明
- 簡単なサンプル
 
を載せています．

間違っていたらご指摘ください．

##環境
私は以下のような環境で動かしています．

Windows 10
Intel(R) Core(TM) i7-6700K CPU
Nvidia GeForce GTX 980
CUDA 9.1
Microsoft Visual Studio Community 2017 VisualStudio/15.0.0+26228.29

##CUDA Stream
###CUDA Streamとは
CUDAではメモリの転送やカーネル実行のスケジュールをStream単位で管理しています．
Streamはデフォルトでは1本で，それをみんなが使っています．
このStreamは実は追加で増やすことができます．
###Streamを生成する
Streamを新しく生成して使用するメリットは待ち時間が少なくなることです．
レジが1つから2つに増えたみたいなイメージです．

Streamを増やす方法は簡単でhost関数で次のように記述してやります．

```C++
cudaStream_t stream;
cudaStreamCreate(&stream);
```
このようにするとストリームが1つ増えて2つになります．
また，こんなような書き方もできます．

```C++
cudaStream_t stream[CPU_THREAD];
for (int i = 0; i < CPU_THREAD; i++) cudaStreamCreate(&stream[i]);
```
`#define CPU_THREAD (8)`などと定義しておけば，CPUのスレッド分ストリームを生成してくれることになります．

###Streamを使用する
次にStreamを使用してみます．
Streamを使用する場面は

- メモリ転送
- カーネル実行

の2つが主です．順に見ていきましょう．

まず，メモリ転送ですが，CUDAのメモリ転送は2種類あります．

```C++
cudaMemcpy(&device, &host, sizeof(type) * N, cudaMemCpyHostToDevice);
cudaMemcpyAsync(&device, &host, sizeof(type) * N, cudaMemCpyHostToDevice, sm);
```
上の方の`cudaMemcpy`はCPUと同期して処理するのが特徴で，デフォルトのストリームを使うことが決まっています．
一方，下の方の`cudaMemcpyAsync`はCPUと非同期で処理するのが特徴で，最後に`cudaStream_t`型の変数を引数に取ることでストリームを選択できます．

カーネル実行では次のように書き分けます．

```C++
kernel <<<grid, block>>> ();
kernel <<<grid, block, 0, stream>>> ();
```
上の方はデフォルトのストリームを使用する場合で，下の方はストリームを選択して使用する場合です．
また，もしもカーネル実行でストリームを選択して使用したら最後に

```C++
cudaStreamSynchronize(stream);
```
と記述して同期を取ってあげなければなりません．

###コンパイルオプション
ストリームの機能を使用するには上記のようにコ―ディングしたうえで以下のコンパイルオプションを追加しなければなりません．

```
--default-stream per-thread
```

##OpenMPと一緒に使う
OpenMPを使ってfor分を並列化する場合は次のように書きます．

```C++
#pragma omp parallel for
for (int i = 0; i < N; i++)
{
    // 何かしらの処理
}
```
たったこれだけでCPUのスレッドが手分けしてN回の**何かしらの処理**を並列に実行してくれます．

さて，この**何かしらの処理**の部分でカーネル実行を行うことを考えましょう．

```C++
#pragma omp parallel for
for (int i = 0; i < N; i++)
{
    kernel <<<grid, block>>> ();
}
```
もしも，このような書き方をしてしまうと，CPUが並列計算してくれている恩恵があまりありません．
なぜなら，全てのカーネル実行においてデフォルトのストリームを使用しているので待ちが発生してしまっているからです．

なので，次のように書き換えます．

```C++
#pragma omp parallel for
for (int i = 0; i < N; i++)
{
    int thread_id = omp_get_thread_num();

    kernel <<<grid, block, 0, stream[thread_id]>>> ();

    cudaStreamSynchronize(stream[thread_id]);
}
```
このようにすることで複数のストリームを使用して効率的に処理することができます．

##サンプル
最後にサンプルをあげておきます．
サンプルの末尾にはコメントアウトで出力結果が書いてあります．ストリームを的確に使用することで2倍以上高速化しています．

とくに意味のないコードですが，お役に立てば幸いです．（コマンドラインオプション忘れずに）


```C++
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <chrono>
#include <stdio.h>

#include <omp.h>

#define CPU_THREAD (8)

__global__ void addKernel(int *c, int *a, int *b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    c[i] += a[i] + b[i];
}

void normal_execution()
{
	fprintf(stdout, "\nnormal execution\n");

	const unsigned int size = 1 << 10;
	static int a_host[size];
	static int b_host[size];
	static int c_host[size];
	int *a_device, *b_device, *c_device;

	cudaMalloc(&a_device, sizeof(int) * size);
	cudaMalloc(&b_device, sizeof(int) * size);
	cudaMalloc(&c_device, sizeof(int) * size);

	// initialize
	for (int i = 0; i < size; i++)
	{
		a_host[i] = i;
		b_host[i] = i + 1;
		c_host[i] = 0;
	}

	dim3 grid(size >> 10, 1, 1);
	dim3 block(1024, 1, 1);

	std::chrono::system_clock::time_point  start, end;
	int elapsed;
	start = std::chrono::system_clock::now();

#pragma omp parallel for
	for (int i = 0; i < 1 << 10; i++)
	{
		cudaMemcpy(a_device, a_host, sizeof(int) * size, cudaMemcpyHostToDevice);
		cudaMemcpy(b_device, b_host, sizeof(int) * size, cudaMemcpyHostToDevice);
		cudaMemcpy(c_device, c_host, sizeof(int) * size, cudaMemcpyHostToDevice);
		addKernel << < grid, block, 0 >> > (c_device, a_device, b_device);
	}

	cudaMemcpy(c_host, c_device, sizeof(int) * size, cudaMemcpyDeviceToHost);

	end = std::chrono::system_clock::now();
	elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
	fprintf(stdout, "%d milli sec.\n", elapsed);

	// check
	for (int i = 0; i < 5; i++)
		fprintf(stdout, "%2d + %2d = %2d\n", a_host[i], b_host[i], c_host[i]);
}

void conccurent_execution()
{
	fprintf(stdout, "\nconcurrent execution\n");

	const unsigned int size = 1 << 10;
	static int a_host[size];
	static int b_host[size];
	static int c_host[size];
	int *a_device, *b_device, *c_device;

	cudaMalloc(&a_device, sizeof(int) * size);
	cudaMalloc(&b_device, sizeof(int) * size);
	cudaMalloc(&c_device, sizeof(int) * size);

	// initialize
	for (int i = 0; i < size; i++)
	{
		a_host[i] = i;
		b_host[i] = i + 1;
		c_host[i] = 0;
	}

	dim3 grid(size >> 10, 1, 1);
	dim3 block(1024, 1, 1);

	cudaStream_t stream[CPU_THREAD];
	for (int i = 0; i < CPU_THREAD; i++) cudaStreamCreate(&stream[i]);

	std::chrono::system_clock::time_point  start, end;
	int elapsed;
	start = std::chrono::system_clock::now();

#pragma omp parallel for
	for (int i = 0; i < 1 << 10; i++)
	{
		int thread_id = omp_get_thread_num();
		cudaMemcpyAsync(a_device, a_host, sizeof(int) * size, cudaMemcpyHostToDevice, stream[thread_id]);
		cudaMemcpyAsync(b_device, b_host, sizeof(int) * size, cudaMemcpyHostToDevice, stream[thread_id]);
		cudaMemcpyAsync(c_device, c_host, sizeof(int) * size, cudaMemcpyHostToDevice, stream[thread_id]);
		addKernel << < grid, block, 0, stream[thread_id] >> > (c_device, a_device, b_device);
		cudaStreamSynchronize(stream[thread_id]);
	}

	cudaMemcpy(c_host, c_device, sizeof(int) * size, cudaMemcpyDeviceToHost);

	end = std::chrono::system_clock::now();
	elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
	fprintf(stdout, "%d milli sec.\n", elapsed);

	// check
	for (int i = 0; i < 5; i++)
		fprintf(stdout, "%2d + %2d = %2d\n", a_host[i], b_host[i], c_host[i]);
}

int main()
{
	normal_execution();
	conccurent_execution();

    return 0;
}

/*

normal execution
172 milli sec.
0 +  1 =  1
1 +  2 =  3
2 +  3 =  5
3 +  4 =  7
4 +  5 =  9

concurrent execution
74 milli sec.
0 +  1 =  1
1 +  2 =  3
2 +  3 =  5
3 +  4 =  7
4 +  5 =  9
*/
```