8. CUDA

8.1 CUDAとは

NVIDIA社のグラフィックスボード(GPU)が搭載されたコンピュータでは、 その高い演算能力を汎用的な科学技術計算に用いることができます。 そのためのプログラミング言語をCUDAと呼びます。
CPUがフロントエンドとして動作し、計算の主要部はGPUが行います。

8.2 CUDAプログラミング

CUDA言語
CUDAのソースコードの拡張子は".cu"、コンパイラ・リンカのコマンドは"nvcc"です。
CUDA言語はCにいくつかの拡張機能を持たせたものですが、言語仕様は基本的にC++です。
ヘッダーファイルは特に必要ありません。
以下の条件を満たすC99スタイルのCのソースコードはそのままnvccでコンパイルすることができます。
(1)引数を持たない関数の宣言は"(void)"ではなく"()"とする。
(2)mallocの返り値は(void *)にキャストする。
C++のソースコードはそのままnvccでコンパイルすることができます。
CUDAは同じソースコードをWindowsとLinuxの両方でコンパイルすることができます。

Compute Capability
グラフィックスボードのハードウェアバージョンを"Compute Capability"と呼びます。
これが上がるごとに計算機能の上限が上がります(必ずしも計算速度が上がる訳ではありません)。
"Kepler"のバージョンは3.X、"Maxwell"のバージョンは5.X、 "Pascal"のバージョンは6.Xです。
グラフィックスボードの"Compute Capability"については このページ を参考にしてください。
機能の違いについては[12]のAppendix Gを参考にしてください。

CUDAプログラミングの指針
CUDAプログラミングでは以下の点が最も重要です。[13]
(1)並列計算できるアルゴリズムを採用する。
(2)CPU/GPU間のデータ転送を最小限にする。
(3)カーネルコードではメモリアクセスをスレッド順とする。(coalescing of memory)
これらを満たさないときは速度は数分の一以下になりGPUを使う意味がなくなります。

単精度と倍精度
GPUコードは倍精度を使用すると極端に性能が下がります。 GPUコードでは計算精度が許す限りなるべく単精度を使ってください。

メモリー転送速度
CUDA付属のサンプルプログラムbandwidthTest.exeを実行するとホスト・デバイス間とデバイス内のメモリー転送速度(バンド幅)を確認することができます。
表8-1はNVIDIA GeForce GTX 1070の結果です。
これからホスト・デバイス間(PCI-Express 3.0)はデバイス内に比べて大幅に遅いことがわかります。
通常はプログラミングしやすいmemory=pageableを使用します。
ホスト・デバイス間のバンド幅がネックになるアプリケーションではGPUの性能は発揮できませんので、 アルゴリズムを見なおす必要があります。

表8-1 メモリー転送速度(bandwidthTestの結果, NVIDIA GeForce GTX 1070)
No.転送方向メモリー転送速度
--memory=pinned--memory=pageable
1host -> device6.4 GB/sec3.8 GB/sec
2device -> host6.4 GB/sec3.7 GB/sec
3device -> device191 GB/sec191 GB/sec

8.3 unified memory

CUDA6からunified memoryが追加されました(managed memoryとも呼ばれることもあります)。 これはhostとdeviceの両方からアクセスできるメモリーであり、 これによりメモリー関係のプログラミングが大幅に簡素化されます。
しかし現状では(2016年現在)まだ未完成の技術であり、特にPascal世代のGPUでは性能が大幅に低下します。
従ってここでは従来の方法(host-device memoryと呼ぶ)を既定値とし、 unified memoryに容易に移行できるようなプログラムを考えます。
なお、unified memoryはCompute Capability 3.0以上の64ビット版でサポートされます。

8.4 自作メモリー関数

メモリー関係の関数をGPU/CPUの別とhost-device memory/unified memory の別の違いを引数で吸収した自作メモリー関数を表8-2とリスト8-1に示します。
これを用いると同じソースコードで異なるアーキテクチャーに対応することができます。
なお、cuda_malloc関数は配列を確保すると同時に0に初期化しています。

表8-2 自作メモリー関数とその実体
自作メモリー関数GPUCPU
host-device memoryunified memory
cuda_malloccudaMalloccudaMallocManagedmalloc
cuda_freecudaFreecudaFreefree
cuda_memsetcudaMemsetcudaMemsetmemset
cuda_memcpycudaMemcpycudaMemcpymemcpy


リスト8-1 自作メモリー関数 (cuda_memory.cu)
     1	/*
     2	cuda_memory.cu
     3	
     4	CUDA memory utilities
     5	*/
     6	
     7	#include <stdlib.h>
     8	#include <string.h>
     9	
    10	// malloc and clear
    11	void cuda_malloc(int gpu, int um, void **ptr, size_t size)
    12	{
    13		if (gpu) {
    14			if (um) {
    15				cudaMallocManaged(ptr, size);
    16			}
    17			else {
    18				cudaMalloc(ptr, size);
    19			}
    20			cudaMemset(*ptr, 0, size);
    21		}
    22		else {
    23			*ptr = malloc(size);
    24			memset(*ptr, 0, size);
    25		}
    26	}
    27	
    28	// free
    29	void cuda_free(int gpu, void *ptr)
    30	{
    31		if (gpu) {
    32			cudaFree(ptr);
    33		}
    34		else {
    35			free(ptr);
    36		}
    37	}
    38	
    39	// memset
    40	void cuda_memset(int gpu, void *ptr, int c, size_t size)
    41	{
    42		if (gpu) {
    43			cudaMemset(ptr, c, size);
    44		}
    45		else {
    46			memset(ptr, c, size);
    47		}
    48	}
    49	
    50	// memcpy
    51	void cuda_memcpy(int gpu, void *dst, const void *src, size_t size, cudaMemcpyKind kind)
    52	{
    53		if (gpu) {
    54			cudaMemcpy(dst, src, size, kind);
    55		}
    56		else {
    57			memcpy(dst, src, size);
    58		}
    59	}

8.5 CUDAプログラミング例

リスト8-2にベクトルの和をCUDAで並列計算するプログラムを示します。
ベクトルの和(C=A+B)を並列計算するには、ループをグリッドとブロックに分割します。


リスト8-2 CUDAプログラム (vadd_cuda.cu)
     1	/*
     2	add two vectors (CUDA)
     3	
     4	Compile,Link:
     5	> nvcc -O3 -o vadd_cuda vadd_cuda.cu cuda_memory.cu
     6	
     7	Usage:
     8	> vadd_cuda [-gpu|-cpu] [-hdm|-um] [-device <device>] <num> <loop>
     9	*/
    10	
    11	// GPU/CPU
    12	__host__ __device__
    13	static void vadd_calc(float a, float b, float *c)
    14	{
    15		*c = a + b;
    16	}
    17	
    18	// GPU
    19	__global__
    20	static void vadd_gpu(int n, const float *a, const float *b, float *c)
    21	{
    22		int tid = threadIdx.x + (blockIdx.x * blockDim.x);
    23		if (tid < n) {
    24			vadd_calc(a[tid], b[tid], &c[tid]);
    25		}
    26	}
    27	
    28	// CPU
    29	static void vadd_cpu(int n, const float *a, const float *b, float *c)
    30	{
    31		for (int i = 0; i < n; i++) {
    32			vadd_calc(a[i], b[i], &c[i]);
    33		}
    34	}
    35	
    36	// GPU/CPU
    37	static void vadd(int gpu, int n, const float *a, const float *b, float *c)
    38	{
    39		if (gpu) {
    40			int block = 256;
    41			int grid = (n + block - 1) / block;
    42			vadd_gpu<<<grid, block>>>(n, a, b, c);
    43		}
    44		else {
    45			vadd_cpu(n, a, b, c);
    46		}
    47	}
    48	
    49	#include <stdlib.h>
    50	#include <stdio.h>
    51	#include <string.h>
    52	#include <time.h>
    53	
    54	extern void cuda_malloc(int, int, void **, size_t);
    55	extern void cuda_free(int, void *);
    56	extern void cuda_memcpy(int, void *, const void *, size_t, cudaMemcpyKind);
    57	
    58	int main(int argc, char **argv)
    59	{
    60		int    gpu = 1;
    61		int    um = 0;
    62		int    device = 0;
    63		int    n = 1000;
    64		int    nloop = 1000;
    65		float  *a, *b, *c;
    66		clock_t t0 = 0, t1 = 0;
    67	
    68		// arguments
    69		while (--argc) {
    70			argv++;
    71			if      (!strcmp(*argv, "-gpu")) {
    72				gpu = 1;
    73			}
    74			else if (!strcmp(*argv, "-cpu")) {
    75				gpu = 0;
    76			}
    77			else if (!strcmp(*argv, "-hdm")) {
    78				um = 0;
    79			}
    80			else if (!strcmp(*argv, "-um")) {
    81				um = 1;
    82			}
    83			else if (!strcmp(*argv, "-device")) {
    84				if (--argc) {
    85					device = atoi(*++argv);
    86				}
    87			}
    88			else if (argc == 2) {
    89				n = atoi(*argv);
    90			}
    91			else if (argc == 1) {
    92				nloop = atoi(*argv);
    93			}
    94		}
    95	
    96		// GPU info and set device
    97		if (gpu) {
    98			int ndevice;
    99			cudaDeviceProp prop;
   100			cudaGetDeviceCount(&ndevice);
   101			if (device >= ndevice) device = ndevice - 1;
   102			cudaGetDeviceProperties(&prop, device);
   103			printf("GPU-%d : %s, C.C.%d.%d, U.M.%s\n", device, prop.name, prop.major, prop.minor, (um ? "ON" : "OFF"));
   104			cudaSetDevice(device);
   105		}
   106	
   107		// alloc device or managed memory
   108		size_t size = n * sizeof(float);
   109		cuda_malloc(gpu, um, (void **)&a, size);
   110		cuda_malloc(gpu, um, (void **)&b, size);
   111		cuda_malloc(gpu, um, (void **)&c, size);
   112	
   113		// alloc host memory
   114		float *h_a = (float *)malloc(size);
   115		float *h_b = (float *)malloc(size);
   116	
   117		// setup problem
   118		for (int i = 0; i < n; i++) {
   119			h_a[i] = i;
   120			h_b[i] = i + 1;
   121		}
   122	
   123		// copy host to device
   124		cuda_memcpy(gpu, a, h_a, size, cudaMemcpyHostToDevice);
   125		cuda_memcpy(gpu, b, h_b, size, cudaMemcpyHostToDevice);
   126	
   127		// timer
   128		t0 = clock();
   129	
   130		// calculation
   131		for (int loop = 0; loop < nloop; loop++) {
   132			vadd(gpu, n, a, b, c);
   133		}
   134		if (gpu) cudaDeviceSynchronize();
   135	
   136		// timer
   137		t1 = clock();
   138	
   139		// copy device to host
   140		float *h_c = (float *)malloc(size);
   141		cuda_memcpy(gpu, h_c, c, size, cudaMemcpyDeviceToHost);
   142	
   143		// sum
   144		double sum = 0;
   145		for (int i = 0; i < n; i++) {
   146			sum += h_c[i];
   147		}
   148	
   149		// output
   150		double exact = (double)n * n;
   151		double sec = (double)(t1 - t0) / CLOCKS_PER_SEC;
   152		printf("nvector=%d nloop=%d %e(%e) %s[sec]=%.3f\n",
   153			n, nloop, sum, exact, (gpu ? "GPU" : "CPU"), sec);
   154	
   155		// free
   156		free(h_a);
   157		free(h_b);
   158		free(h_c);
   159		cuda_free(gpu, a);
   160		cuda_free(gpu, b);
   161		cuda_free(gpu, c);
   162	
   163		return 0;
   164	}

ソースコードの説明
12行目: 関数vadd_calcはGPUとCPUの両方から呼ばれますのでキーワード"__host__ __device__"を付ける必要があります。
19行目: カーネルコード(GPUコード)にはキーワード"__global__"を付ける必要があります。
22行目: スレッド番号を計算する定型文です。 ブロックまたはグリッドが2次元または3次元のときはレファレンスマニュアルを参考にしてください。
23行目: 配列の大きさはブロックの大きさの倍数とは限りませんのでこの条件判定が必要です。
24行目: 配列"a","b","c"はスレッド番号に関して連続的にアクセスされていることに注意してください(coalescing of memory)。 このようになっていないときは計算時間が大幅に増えます。
42行目: GPUでは<<< >>>の中(execution configuration)にグリッドとブロックを指定し、 カーネル関数を呼びます。 カーネル関数の引数の中の配列はdevice memoryです。
45行目: CPUの引数の中の配列はhost memoryです。
109-111行目: device memoryを確保します。unified memoryのときはCPUから直接参照することができます。
114-125行目: device memoryの内容を初期化します。CPU側で配列を作成しGPUに転送しています。 device memoryに一定値を代入するときはCPUを経由せずcuda_memset関数で直接代入することもできます。
134行目: unified memoryではカーネルコードが終了した後、 CPUからunified memoryを参照する前にcudaDeviceSynchronize関数でCPU/GPU間の同期をとる必要があります。 host-device memoryでは同期は不要です。 なお、cudaMemcpy関数を呼び出すと暗黙で同期がとられます。
141行目: GPUからCPUに計算結果を転送します。
159-161行目: device memoryを解放します。

CPU版とGPU版の開発
CUDAプログラムでは変数の指すメモリーがCPUとGPUのどちらにあるかを常に念頭に置く必要があり、 ここでバグが発生しやすくなります。
そこで開発にあたっては先ずメモリー関係をリスト8-1のメモリー関数で記述したCPU版を作成し、 計算ロジックが正しく動作していることをCPUで十分確認してから、 メモリー関数の引数をGPUに変えてGPUで動作確認します。
このようにすればGPU版の開発時間を最小にすることができます。
ただしhost-device memoryではCPUからGPUメモリーを参照することができないために、 cuda_memcpy関数を用いてGPUからCPUにメモリーをコピーする必要が発生することがあります。

変数の命名ルール
host-device memoryでは同じ意味をもつ変数がCPUとGPUの両方にあるときは、 CPUの変数の頭に"h_"(host memoryの意味)、 GPUの変数の頭に"d_"(device memoryの意味)を付ける方法がよく用いられます。
このようにすればその変数がCPUにあるかGPUにあるかが一目でわかります。

CPU/GPUコードの共通化
本プログラムの特徴はGPUではvadd_gpu関数を呼び出し、CPUではvadd_cpu関数を呼び出し、 それぞれから同じvadd_calc関数を呼び出していることです。
計算ロジックはvadd_calc関数にすべて集約され、 vadd_gpu, vadd_cpu関数は変数の受け渡しとスレッド番号の計算という機械的な処理をしているだけです。 実際vadd_gpu関数とvadd_cpu関数の違いは22-23行目と31行目だけです。
このようにすれば計算ロジックとGPUメモリーの操作を切り離すことができ、 CPUで十分デバッグすることができ、GPUコードの追加作業を最小限にすることができます。
単にベクトルの和を計算するのには冗長なように見えますが、 計算処理が複雑になると開発効率に差が出るようになります。

8.6 コンパイル・リンク・実行

コンパイル・リンク方法
コンパイル・リンク方法はVC++とgccで同じです。
> nvcc -O3 -o vadd_cuda -arch=sm_30 vadd_cuda.cu cuda_memory.cu
VC++では多数の"warning C4819"が出ることがあります。 そのときは以下のようにコンパイルオプションを追加してください。
> nvcc -O3 -Xcompiler "/wd4819" -o vadd_cuda -arch=sm_30 vadd_cuda.cu cuda_memory.cu
なお、Windows環境では"C:\Windows\Temp"フォルダに書き込みをする権限が必要です。 コンパイラーのエラーが出たら確認してください。

プログラムの実行方法
プログラムの実行方法は以下の通りです。
> vadd_cuda [-gpu|-cpu] [-hdm|-um] [-device デバイス番号] 配列の大きさ 繰り返し回数
例えば以下のようになります。
> vadd_cuda 10000000 1000 (GPUで計算するとき)
> vadd_cuda -um 10000000 1000 (GPUのunified memoryで計算するとき)
> vadd_cuda -cpu 10000000 1000 (CPUで計算するとき)
> vadd_cuda -device 1 10000000 1000 (デバイス番号1のGPUで計算するとき)
繰り返し回数は計算時間の測定誤差を小さくするためです。
デバイス番号については下記のページの「2.4 グラフィックスボードのデバイス番号」を参考にしてください。
http://www.e-em.co.jp/gpu/gpu.htm

8.7 CUDAの計算時間

表8-3にWindowsでの計算時間を示します。
配列の大きさ(=N)と繰り返し回数(=L)の積は一定(=10^10)です。従って全体の演算量は同じです。
No.1-3ではGPUはCPUより大幅に速いことがわかります。
No.4ではGPUはカーネル起動回数が増えるために遅くなります。 リスト8-2よりカーネル起動回数は繰り返し回数と同じです。 繰り返し回数が1,000,000回のとき2秒程度余分に時間がかかっていることから、 カーネル起動のオーバーヘッドはおよそ2μsecと評価することができます。 (通常のアプリケーションではカーネルをこのように多数回呼ぶことは少ないです)
なお、Pascal世代のGPU(C.C.6.X)ではunified memoryはhost-device memoryに比べて数十倍遅くなります。

表8-3 CUDAの計算時間(Windows, ()内はCPU1コアとの速度比)
No.配列の大きさN繰り返し回数LCPU1コアGPU
host-device memory
110,000,0001,00011.03秒 (1.0)0.58秒 (19.0)
21,000,00010,0007.12秒 (1.0)0.60秒 (11.9)
3100,000100,0002.83秒 (1.0)0.38秒 (7.4)
410,0001,000,0001.65秒 (1.0)2.38秒 (0.7)