*CUDAで行列演算 その1 加減算 [#z5d8526e]
CUDAでのプログラミングの基礎を学ぶために，数値解析などでよく使う行列の演算(和，差，積，行列式，逆行列など)をC for CUDAで実装してみる．まず最初に行列同士の加減算を行う．C for CUDAでの基本的な処理手順は以下のようになる．
+ホスト側メモリ確保,入力データ準備
+デバイス側メモリ確保
+ホストからデバイスへデータ転送
+CUDAカーネル呼び出し
+結果をデバイスからホストへ転送
+メモリ解放
**ホスト側メモリ確保 [#o1d84552]
行列を1次元配列上に格納する．演算対象の行列をAとB，結果を格納する行列をCとし，行列の大きさをNrow, Ncolとする．このとき，rowメジャーとすると
| A0 A1 A2 A3 |
| A4 A5 A6 A7 |
| A8 A9 A10 A11 |
の順に格納される．ここで，Nrow = 3, Ncol = 4とした．
このとき，i行j列の要素はA[i*Ncol+j]で取得できる．
ホストメモリを格納する変数名はhA,hBのようにhを接頭子とすると，
ホストメモリは次のようにして確保できる．
#code(C){{
unsigned int size= sizeof(float)*nrow*ncol;
float *hA = (float*)malloc(size);
float *hB = (float*)malloc(size);
float *hC = (float*)malloc(size);
unsigned int size_byte = sizeof(float)*nrow*ncol;
float *hA = (float*)malloc(size_byte);
float *hB = (float*)malloc(size_byte);
float *hC = (float*)malloc(size_byte);
}}
メモリの確保はCPU上で行われるので，C/C++言語の他の手法でも問題ない．．
また，CPU側で処理される関数上でSTLのvector等を使って確保した場合は，
#code(C,nooutline,nomenu){{
vector<float> A;
A.resize(size);
CuTest(&A[0]);
}}
CuTestはCudaのホスト側関数である．
デバイスメモリに転送する前に，
上記の行列の要素アクセスを使って配列に行列の入力値を代入しておく．
ここでは実験的にチェックしたいので，乱数を用いて
#code(C){{
void RandomInit(float* data, int size, float max)
{
for(int i = 0; i < size; ++i){
data[i] = max*(rand()/(float)RAND_MAX);
}
}
}}
のような関数で初期化を行う．
**デバイス側メモリ確保 [#m776bcca]
CUDAのプログラミングモデルでは，ホストとデバイスがそれぞれ自身のメモリ
(ホストメモリとデバイスメモリ)を持つことを仮定しています．
CUDAカーネルはデバイスメモリの中だけ操作することが可能であり，
カーネルに渡すデータはデバイスメモリに確保されている必要があります．
デバイスメモリは，linear memory か CUDA arrays のどちらかで確保される．
よく使われるのはlinear memoryだと思われる．
CUDA arraysは内部のデータレイアウトが不透明であり，テクスチャフェッチに最適化されている．
ここでは，linear memoryを使用する．
メモリ確保命令は cudaMalloc である．
#code(C){{
// デバイスメモリの確保
float *dA, *dB, *dC;
cutilSafeCall(cudaMalloc((void**)&dA, size));
cutilSafeCall(cudaMalloc((void**)&dB, size));
cutilSafeCall(cudaMalloc((void**)&dC, size));
}}
cutilSafeCallはSDKが提供するデバッグ用の関数ラッパである．
**ホストからデバイスへデータ転送 [#m7009637]
デバイスメモリを確保したら，ホストメモリからデバイスメモリへデータをコピーする．コピーするには cudaMemcpy 関数を用いる．
#code(C){{
// ホストからデバイスへの転送
cutilSafeCall(cudaMemcpy(dA, hA, size, cudaMemcpyHostToDevice));
cutilSafeCall(cudaMemcpy(dB, hB, size, cudaMemcpyHostToDevice));
}}
cudaMemcpyは第一引数に転送先，第二引数に転送元，第三引数にデータサイズ，最後の引数で転送方向を指定する．cudaMemcpyHostToDeviceでホストからデバイス，
cudaMemcpyDeviceToHostでデバイスからホストへのコピーである．
**CUDAカーネル呼び出し [#l7a8f65a]
CUDAカーネルの呼び出しには，<<< ... >>> を用いる．
カーネル関数名をmatrixAddとすると，
#code(C){{
dim3 block(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid((nrow+block.x-1)/block.x, (ncol+block.y-1)/block.y);
matrixAdd<<< grid, block >>>(dC, dA, dB, nrow, ncol);
}}
のようにして呼び出す．
ここで，BLOCK_SIZE=16とした．
"<<< >>>" の中にはグリッド内のブロック数(grid,2Dまで)，ブロック内のスレッド数(block,3Dまで)とオプションで，ストリームID，シェアードメモリサイズも指定できる．今回の演算では行列の縦横のサイズ(2D)に合わせてスレッドを確保する．
ブロック内のスレッド数を固定としているので，ブロック数をnrow/block.x, ncol/block.yとすればよいである．ただし，nrow < block.x or ncol < block.yの場合を考慮して，上記のようにしてある．
dim3はCUDAのベクトル型変数でグリッドやブロックのサイズを指定するのに用いる．他のベクトル型については，
[[CUDAのベクトル型]]にまとめる．
呼び出される側のカーネル関数は，
#code(C){{
__global__
void matrixAdd(float* C, float* A, float* B, int nrow, int ncol)
{
int row = blockIdx.y*blockDim.y+threadIdx.y;
int col = blockIdx.x*blockDim.x+threadIdx.x;
int idx = row*ncol+col;
if(row < nrow && col < ncol){
C[idx] = A[idx]+B[idx];
}
}
}}
のように定義する．
__global__修飾子については[[CUDA関数修飾子]]を参照．
また，blockIdxやblockDim，threadIdxは[[CUDAビルトイン変数]]である．
SDKでカーネル実行エラーをチェックするには，カーネル呼び出し直後に，
cutilCheckMsg("Kernel execution failed");
を呼び出す．
**結果をデバイスからホストへ転送 [#ga6e2975]
演算結果をデバイスからホストへ転送する．
#code(C){{
// デバイスからホストへ結果を転送
cutilSafeCall(cudaMemcpy(hC, dC, size, cudaMemcpyDeviceToHost));
}}
**メモリ解放 [#bf7a7561]
すべて終了したらホストメモリ，デバイスメモリを解放する．
#code(C){{
// ホストメモリ解放
free(hA);
free(hB);
free(hC);
// デバイスメモリ解放
cutilSafeCall(cudaFree(dA));
cutilSafeCall(cudaFree(dB));
cutilSafeCall(cudaFree(dC));
}}
デバイスメモリの解放にはcudaFreeを用いる．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
