#繰り返し処理をCUDAで書く（配列同士の足し算）
繰り返し処理を並列化するとき、その目的によって並列化の手段が違う。役に立つかどうかは置いといて、ここでは基本的なことをまとめます。
##目的
配列同士の足し算ベクトル和を考えます。例えばそれぞれ、大きさが $N$ 個の配列 $A$ と配列 $B$ の各要素の足し算の結果を、配列 $C$ に入れたいという問題を考えたとき、C言語では以下のように書けます。
for(int i=0;i<N;i++) C[i] = A[i] + B[i];
これをCUDAで並列化するにはどうするかという一例を解説します。
##CUDAプログラミングにおける基礎知識
###スレッド管理について
GPUは多くのスレッドを持っています。CUDAプログラミングではそのスレッドをうまく利用して大規模な計算を高速で処理します。すなわち、並列処理を行う場合、この「スレッド」を意識してプログラムを組む必要があります。
CUDAプログラミングではGPUがもつ大量のスレッドのうち必要な数をプログラマーが指定して並列処理を行います。
その際、プログラマーはCUDAで決められたスレッドの管理方式に従って、その数を指定します。スレッドは３つの階層で管理され、それぞれ最下層「スレッド」、中層「ブロック」、最上層「グリッド」という名前が付けられています。
- スレッド...階層の最小単位
- ブロック...複数のスレッドの塊
- グリッド...ブロック全体の集合
例えば、全部で $512×2048 = 1048576$ 個のスレッドを使いたいとします。このとき、ブロックサイズ（ブロック一つ当たりのスレッド数）を $512$ とすると、必要なグリッドサイズは、 $(512 \times 2048)/512 = 2048 $となります。したがって、グリッドサイズを $2048$ 、ブロックサイズを $512$ と指定すれば、所望のスレッドを呼び出すことができます。
グリッドの中でブロックは、「x方向、y方向、z方向」と３次元的に配置されています。そしてブロックの中でもスレッドは同じようにx,y,zと３次元的に配置されています。
CUDA7.5で提供されているSampleの一つに、deviceQueyというソリューションファイルがありますが(参考：CUDA 7.5 をWindows10にインストールし、Sampleコードを実行するまで)、それを実行させると、お使いのGPUで、実際にスレッドをどのように管理できるかがわかります。
ちなみに私の環境(GeForce GTX 750 Ti)でdeviceQueyを実行すると、以下のような出力がでます。
CUDA Device Query (Runtime API) version (CUDART static linking)
Detected 1 CUDA Capable device(s)
Device 0: "GeForce GTX 750 Ti"
CUDA Driver Version / Runtime Version 7.5 / 7.5
CUDA Capability Major/Minor version number: 5.0
Total amount of global memory: 2048 MBytes (2147483648 bytes)
( 5) Multiprocessors, (128) CUDA Cores/MP: 640 CUDA Cores
GPU Max Clock rate: 1085 MHz (1.08 GHz)
Memory Clock rate: 2700 Mhz
Memory Bus Width: 128-bit
L2 Cache Size: 2097152 bytes
Maximum Texture Dimension Size (x,y,z) 1D=(65536), 2D=(65536, 65536), 3D=(4096, 4096, 4096)
Maximum Layered 1D Texture Size, (num) layers 1D=(16384), 2048 layers
Maximum Layered 2D Texture Size, (num) layers 2D=(16384, 16384), 2048 layers
Total amount of constant memory: 65536 bytes
Total amount of shared memory per block: 49152 bytes
Total number of registers available per block: 65536
Warp size: 32
Maximum number of threads per multiprocessor: 2048
Maximum number of threads per block: 1024
Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
Max dimension size of a grid size (x,y,z): (2147483647, 65535, 65535)
Maximum memory pitch: 2147483647 bytes
Texture alignment: 512 bytes
Concurrent copy and kernel execution: Yes with 1 copy engine(s)
Run time limit on kernels: Yes
Integrated GPU sharing Host Memory: No
Support host page-locked memory mapping: Yes
Alignment requirement for Surfaces: Yes
Device has ECC support: Disabled
CUDA Device Driver Mode (TCC or WDDM): WDDM (Windows Display Driver Model)
Device supports Unified Addressing (UVA): Yes
Device PCI Domain ID / Bus ID / location ID: 0 / 1 / 0
Compute Mode:
< Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >
deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 7.5, CUDA Runtime Version = 7.5, NumDevs = 1, Device0 = GeForce GTX 750 Ti
Result = PASS
ここで、以下の三行に注目すると、CUDA7.5でGeForce GTX 750 Tiのスレッドは以下のように管理できることがわかります。
Maximum number of threads per block: 1024
Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
Max dimension size of a grid size (x,y,z): (2147483647, 65535, 65535)
- ブロック一つあたりで扱える最大スレッド数 $1024$
- ブロックの中で扱える次元数 $1024×1024×64$
- グリッドの中で扱える次元数 $214743647×65535×65535$
この環境でスレッドすべてを使うことを考えると、 $214743647×65535×65535$ のブロックを用意して、ブロックには最大でスレッドを $1024$ 個ずつもたせることができるので、このGPUを使ってCUDAで同時に処理できる最大のスレッド数は $214743647×65535×65535×1024$ 個（ $\approx$ $10^{20}$ 個 ）となります。
###ホスト、デバイス
CUDAプログラミングではCPUのことを「ホスト」、GPUのことを「デバイス」と呼び、区別します。
ホストで作られた命令をデバイスに渡して並列処理を行い、その結果をデバイスからホストへ移してホストによってその結果を出力するのが、CUDAプログラミングの基本的な流れです。
###dim3型 変数
ホスト側でデバイスのスレッドを指定するとき、dim3型という3次元ベクトルの変数を使って指定します。
###ビルドイン変数
デバイス側にスレッドが用意されると、各スレッドには、自分の位置を示すIDが降られます。このID番号は各スレッドに用意された「ビルドイン変数」というものに格納され、自由に呼び出すことができます。
この変数は宣言せずに利用できますが、内容を書き換えることはできません。
###カーネル関数
CUDAで書かれるコードは、大きくわけるとホストで処理される部分と、デバイスで処理される部分の二つに分類できます。
デバイスで行う並列処理を記述した関数のことを「カーネル関数」といい、プログラム上ではそれとわかるよう宣言をします。
##コード
では、具体的なCUDAのコードを見ていきます。
ここで行う処理は、それぞれ要素数が $N$ 個の配列 $A$ と $B$ を用意し、 $A$ にスカラー値 $k$ を掛けて、それと $B$ を足したものを 同じ要素数の配列 $C$ に代入するというものです。
#include <stdio.h>
//カーネル関数vec_sumの宣言
__global__
void vec_sum(float k, float *a, float *b, float *c)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
c[i] = k*a[i] + b[i];
}
int main(void)
{
//N = 512×2048
int N = 1<<20;
//a,b,cはホスト用、d_a,d_b,d_cはデバイス用のポインタ
float *a, *b, *c *d_a, *d_b, *d_c;
//ホスト側の配列を用意
a = (float*)malloc(N*sizeof(float));
b = (float*)malloc(N*sizeof(float));
c = (float*)malloc(N*sizeof(float));
//デバイス側の配列を用意
cudaMalloc(&d_a, N*sizeof(float));
cudaMalloc(&d_b, N*sizeof(float));
cudaMalloc(&d_c, N*sizeof(float));
//a,bの配列にそれぞれ1，2を代入し、cを初期化
for (int i = 0; i < N; i++) {
a[i] = 1.0f;
b[i] = 2.0f;
c[i] = 0.0f;
}
//ホスト側の配列の内容をデバイス側にコピー
cudaMemcpy(d_a, a, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_b, b, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_c, c, N*sizeof(float), cudaMemcpyHostToDevice);
//スレッドの設定
int blocksize = 512;
//ブロックあたりのスレッド数（blocksize)を512、
//ブロックの総数（gridsize）をN/512用意する
//したがって総スレッド数は blocksize × gridsize = N 個
dim3 block (blocksize, 1, 1);
dim3 grid (N / block.x, 1, 1);
// カーネル関数の呼び出し
vec_sum<<<grid, block>>>(2.0f, d_a, d_b,d_c);
//計算結果をホストへコピー
cudaMemcpy(c, d_c, N*sizeof(float), cudaMemcpyDeviceToHost);
float maxError = 0.0f;
//計算結果の確認
for (int i = 0; i < N; i++) maxError = max(maxError, abs(c[i]-4.0f));
printf("Max error: %f", maxError);
//メモリの開放
free(a);
free(b);
free(c);
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
return 0;
}
- カーネル関数部分
__global__
void vec_sum(float k, float *a, float *b, float *c)
{
int i = blockIdx.x*blockDim.x + threadIdx.x;
c[i] = k*a[i] + b[i];
}
コードの冒頭で宣言されているこの関数が、デバイス側で実行されるカーネル関数です。
__global__ は関数に対する修飾子で、デバイス側で呼び出したいカーネル関数をこの下に宣言します。ここでは関数vec_sumがカーネル関数となります。
カーネル関数の中に登場するblockIdx.x,blockDim.x,thereadIdx.xは「ビルドイン変数」で、これらは宣言せずに利用できます。
このビルドイン変数には処理を実行するスレッドのIDの値が書き込まれています。ビルドイン変数は読み込み専用で、内容を書き換えることはできません。
それぞれの変数の範囲は、指定したスレッド数によって定まります。
今、グリッドサイズ $2048$ 、ブロックサイズ $512$ と指定しています。そのとき、blockIdx.x,blockDim.x,thereadIdx.xそれぞれの値は以下のようになります。
- blockIdx.x $= 0$ ~ $2047$
- blockDim.x $= 512$
- threadIdx.x $= 0$ ~ $511$
したがって、カーネル関数におけるi = blockIdx.x*blockDim.x + threadIdx.xは、スレッドに応じて $0$ ~ $1048575$ の値をとります。
よって、このカーネル関数が呼び出されたスレッドでは、それに対応するビルドイン変数によって i の値が決定し、その値に応じた配列要素の演算を行うということが分かります。
ビルドイン変数はこのほかにも
-
gridDim.x,gridDim.y,gridDim.z -
blockIdy.y,blockDim.y,thereadIdy.y -
blockIdz.z,blockDim.z,thereadIdz.z
などが利用できます。これらを駆使して多次元的に並列処理を行うこともできますが、今回は単純な例なので1次元（ｘ軸）だけを使って演算処理をします。
- スレッドの設定
int blocksize = 512;
dim3 block (blocksize, 1, 1);
dim3 grid (N / block.x, 1, 1);
カーネル関数を呼び出すときに指定するスレッド数のためにベクトル変数dim3を宣言します。dim3は上記のように
dim3 hoge(x,y,z);
と宣言されます。ベクトルの各要素x,y,zにはそれぞれの次元の最大値を指定します。また、それぞれの要素はhoge.x,hoge.y,hoge.zで呼び出すことができて、以下のように内容を変更することもできます。
hoge.x = 20;
hoge.y = 2;
hoge.z = 1;
今回の例では各要素の値は以下のようになっています。
- block.x $= 512$
- block.y $= 1$
- block.z $= 1$
- grid.x $= N/512 = 2048$
- grid.y $= 1$
- grid.z $= 1$
- カーネル関数の呼び出し
ホスト側で用意した配列を、デバイス側で用意した配列にコピーしたあと、カーネル関数を呼び出します。そのときグリッドとブロックの数を指定してデバイス側の変数を引数としします。
vec_sum<<<grid, block>>>(2.0f, d_a, d_b,d_c);
Cとの違いは3重のカッコ記号<<<g,b>>>ですね。このカッコによって、カーネル関数を呼び出すスレッドを指定しています。この例ではグリッド内のブロック数を2048個、ブロック内のスレッド数を512個、それぞれ一次元軸上に指定しています。
実は、1次元だけしか使わない場合、dim3で宣言した変数を使わなくても、以下のように直接書くこともできます。
vec_sum<<<2048, 512>>>(2.0f, d_a, d_b,d_c);
##実行結果
さて、先に示したコードを実行してみます。このままだと面白くないので、すこし手を加えて時間計測の結果を表示するようにして、その結果を表示します。すると以下のようになりました。
cpu time:0.003000 sec
gpu time:0.000000 sec
Max error: 0.000000
cpu timeが、ホスト側で、要素が$512 \times 2048$個からなる ベクトル和の計算処理にかかった時間、gpu time がデバイス側での計算処理にかかった時間です。ホスト側で行った計算は以下になります。
for(int i=0;i<N;i++) C[i] = 2.0f*A[i] + B[i];
100万回の処理ではcpuでもほぼ一瞬で終わってしまいますが、gpuに至っては少なくとも計測不能なほど速いということがこの結果からわかります。
##その他、疑問など
上の結果のような100万回の処理ではgpuのありがたみがどうも見えてこないので、Nの数を1000万や１億に増やして実行したのですが、gpu側で計算結果にエラーが発生してしまい、うまくいきませんでした。
なにが原因なのかを知るには、もうすこし勉強する必要がありそうです。
次の更新では、総和処理の並列化について書きたいと思います。
##参考サイト・文献
[1] 著：John Cheng,Max Grossman,Ty McKercher,訳：株式会社クイープ,監訳：森野慎也, CUDA C プロフェッショナル プログラミング, 株式会社impres, 2015
[2] 青木尊之，額田彰，はじめてのCUDAプログラミング，工学社，2009
[3] Mark Harris , An Easy Introduction to CUDA C and C++, https://devblogs.nvidia.com/parallelforall/easy-introduction-cuda-c-and-c/ ,2012 (2016/07/6 アクセス)