NVIDIA社のグラフィックスボード(GPU)が搭載されたコンピュータでは、
その高い演算能力を汎用的な科学技術計算に用いることができます。
そのためのプログラミング言語をCUDAと呼びます。
CPUがフロントエンドとして動作し、計算の主要部はGPUが行います。
CUDA言語
CUDAのソースコードの拡張子は".cu"、コンパイラ・リンカのコマンドは"nvcc"です。
CUDA言語はCにいくつかの拡張機能を持たせたものですが、言語仕様は基本的にC++です。
ヘッダーファイルは特に必要ありません。
以下の条件を満たすC99スタイルのCのソースコードはそのままnvccでコンパイルすることができます。
(1)引数を持たない関数の宣言は"(void)"ではなく"()"とする。
(2)mallocの返り値は(void *)にキャストする。
C++のソースコードはそのままnvccでコンパイルすることができます。
CUDAは同じソースコードをWindowsとLinuxの両方でコンパイルすることができます。
Compute Capability
グラフィックスボードのハードウェアバージョンを"Compute Capability"（C.C.と略します）と呼びます。
これが上がるごとに計算機能の上限が上がります(必ずしも計算速度が上がる訳ではありません)。
グラフィックスボードの"Compute Capability"については[25]を参考にしてください。
機能の違いについては[12]の Appendix "Compute Capabilities" を参考にしてください。
その他のハードウェアの違いについては[26]を参考にしてください。
CUDAプログラミングの指針
CUDAプログラミングでは以下の点が最も重要です。[13]
(1)並列計算できるアルゴリズムを採用する。
(2)CPU/GPU間のデータ転送を最小限にする。
(3)カーネルコードではメモリアクセスをスレッド順とする。(coalescing of memory)
これらを満たさないときは速度は数分の一以下になりGPUを使う意味がなくなります。
単精度と倍精度
一部の機種を除いて倍精度を使用すると極端に性能が下がります。
計算精度が許す限りなるべく単精度を使ってください。
メモリー転送速度
CUDA付属のサンプルプログラム[24]の bandwidthTest.exe
を実行するとホスト・デバイス間とデバイス内のメモリー転送速度(バンド幅)を確認することができます。
表8-1はNVIDIA GeForce RTX 2060の結果です。
これからホスト・デバイス間(PCI-Express 3.0)はデバイス内に比べて大幅に遅いことがわかります。
通常はプログラミングしやすいmemory=pageableを使用します。
ホスト・デバイス間のバンド幅がネックになるアプリケーションではGPUの性能は発揮できないので、
アルゴリズムを見なおす必要があります。
| No. | 転送方向 | メモリー転送速度 | |
|---|---|---|---|
| --memory=pinned | --memory=pageable | ||
| 1 | host -> device | 6.7 GB/sec | 6.5 GB/sec |
| 2 | device -> host | 6.7 GB/sec | 6.3 GB/sec |
| 3 | device -> device | 231.1 GB/sec | 235.9 GB/sec |
CUDA6から unified memory が追加されました(managed memoryとも呼ばれることもあります)。
これはhostとdeviceの両方からアクセスできるメモリーであり、
これによりメモリー関係のプログラミングが大幅に簡素化されます。
ここでは従来の方法(host-device memoryと呼ぶ)を既定値とし、
unified memory に容易に移行できるようなプログラムを考えます。
メモリー関係の関数をGPU/CPUの別とhost-device memory/unified memory
の別の違いを引数で吸収した自作メモリー関数を表8-2とリスト8-1に示します。
これを用いると同じソースコードで異なるアーキテクチャーに対応することができます。
なお、cuda_malloc関数は配列を確保すると同時に0に初期化しています。
| 自作メモリー関数 | GPU | CPU | |
|---|---|---|---|
| host-device memory | unified memory | ||
| cuda_malloc | cudaMalloc | cudaMallocManaged | malloc |
| cuda_free | cudaFree | cudaFree | free |
| cuda_memset | cudaMemset | cudaMemset | memset |
| cuda_memcpy | cudaMemcpy | cudaMemcpy | memcpy |
リスト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-2にベクトル和をCUDAで並列計算するプログラムを示します。
ベクトル和(C=A+B)を並列計算するには、ループをグリッドとブロックに分割します。
リスト8-2 CUDAプログラム (cuda_vadd.cu)
1 /*
2 cuda_vadd.cu (CUDA)
3
4 Compile,Link:
5 > nvcc -O2 -o cuda_vadd cuda_vadd.cu cuda_memory.cu
6
7 Usage:
8 > cuda_vadd [-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コードの追加作業を最小限にすることができます。
単にベクトルの和を計算するのには冗長なように見えますが、
計算処理が複雑になると開発効率に差が出るようになります。
コンパイル・リンク方法
コンパイル・リンク方法はVC++とgccで同じです。
> nvcc -O2 -o cuda_vadd cuda_vadd.cu cuda_memory.cu
VC++では多数の"warning C4819"が出ることがあります。
そのときは以下のようにコンパイルオプションを追加してください。
> nvcc -O2 -Xcompiler "/wd4819" -o cuda_vadd cuda_vadd.cu cuda_memory.cu
なお、Windows環境では"C:\Windows\Temp"フォルダに書き込みをする権限が必要です。
コンパイラーのエラーが出たら確認してください。
プログラムの実行方法
プログラムの実行方法は以下の通りです。
> cuda_vadd [-gpu|-cpu] [-hdm|-um] [-device デバイス番号] 配列の大きさ 繰り返し回数
例えば以下のようになります。
> cuda_vadd 10000000 1000 (GPUで計算するとき)
> cuda_vadd -um 10000000 1000 (GPUのunified memoryで計算するとき)
> cuda_vadd -cpu 10000000 1000 (CPUで計算するとき)
> cuda_vadd -device 1 10000000 1000 (デバイス番号1のGPUで計算するとき)
繰り返し回数は計算時間の測定誤差を小さくするためです。
表8-3にベクトル和の計算時間を示します。
配列の大きさ(=N)と繰り返し回数(=L)の積は一定(=1010)です。
従って全体の演算量は同じです。
No.1-2ではGPUはCPUより大幅に速いことがわかります。
また、host-device memory と unified memory の計算時間は同じです。
No.3-4ではGPUはカーネル起動回数が増えるために遅くなります。
リスト8-2よりカーネル起動回数は繰り返し回数と同じです。
No.4では繰り返し回数が1,000,000回のとき7秒程度余分に時間がかかっていることから、
カーネル起動のオーバーヘッドはおよそ7μsecと評価することができます。
(通常のアプリケーションではカーネルをこのように多数回呼ぶことは少ないです)
GPUの性能を生かすには総スレッド数（=ブロック数×スレッド数）が1,000,000以上が望ましいです。
なお、No.3-4のCPU1コアでは使用メモリーがCPUのキャッシュに入るために、
SIMDの章で述べたように自動ベクトル化が行われ計算時間が短縮されています。
なお、Pascal世代のGPU(C.C.6.X)ではunified memoryはhost-device memoryに比べて数十倍遅くなります。
| No. | 配列の大きさN | 繰り返し回数L | CPU1コア | GPU host-device memory | GPU unified memory |
|---|---|---|---|---|---|
| 1 | 10,000,000 | 1,000 | 5.46秒 (1.0) | 0.49秒 (11.1) | 0.49秒 (11.1) |
| 2 | 1,000,000 | 10,000 | 5.98秒 (1.0) | 0.52秒 (11.5) | 0.52秒 (11.5) |
| 3 | 100,000 | 100,000 | 1.14秒 (1.0) | 0.78秒 (1.5) | 0.79秒 (1.4) |
| 4 | 10,000 | 1,000,000 | 0.91秒 (1.0) | 7.61秒 (0.1) | 7.74秒 (0.1) |
リスト8-2の40行目ではブロックサイズ(1ブロックのスレッド数)を256としました。
ブロックサイズを変えたときの計算時間を表8-4に示します。
ここで計算条件は表8-3のNo.2と同じです。
これからブロックサイズを変えても計算時間があまり変わらないことがわかります。
多くのケースで256が最適なので、256を推奨します。
なお、ブロックサイズは2のべき乗が普通です。上限は1024です。
（[12]Appendix, "Compute Capabilities" の "Maximum number of threads per block"参照）
| ブロックサイズ | 計算時間 |
|---|---|
| 32 | 0.62秒 |
| 64 | 0.52秒 |
| 128 | 0.52秒 |
| 256 | 0.52秒 |
| 512 | 0.52秒 |
| 1024 | 0.53秒 |
GPU版プログラムを実行するとき複数のグラフィックスボードが実装されているとそれらにデバイス番号が割り当てられます。
既定状態では性能の高い順にデバイス番号0,1,...が割り当てられ、通常はそのままでかまいませんが、
何らかの理由でデバイス番号を変更する必要があるときは以下の方法のいずれかを行ってください。
(1) NVIDIA コントロールパネル
画面の空いた部分を右クリックして[NVIDIAコントロールパネル]→[3D設定]→[3D設定の管理]の
[グローバル設定]タブの[CUDA - GPU]の[設定]で使用するグラフィックスボードのみをONにしてください。
上から順にデバイス番号0,1,...が割り当てられます。
(2) デバイスマネージャー
Windowsの[設定]→[デバイス]→[デバイスマネージャー]を起動し、
[ディスプレイアダプター]のうち使用しないものを右クリックして[無効]としてください。
(3) 環境変数CUDA_VISIBLE_DEVICES
Windowsの[設定]→[システム]→[バージョン情報]→[システム情報]→[システムの詳細設定]を起動し、
[詳細設定]タブの[環境変数]をクリックし、[システム環境変数]の[新規]をクリックし、
[変数名]に"CUDA_VISIBLE_DEVICES"、[変数値]に以下のような適当な数値を入力してください。