*CUDAで行列演算 その2 乗算 [#h119404b]
CUDAで行列の乗算を行います．このページで示すコードは CUDA Programming Guideの
3.2.2 Shared Memoryとほとんど同じですのでそちらも参照してください．
**行列の定義 [#keeb1364]
行列の積は和や差と異なり，それぞれの行列のサイズが異なります．
A(N×M) * B(M×L) = C(N×L)
[[CUDAで行列演算：加減算]]ではすべての行列のサイズが等しかったので楽でしたが，
行列のサイズが異なると管理が面倒になるので，新しくMatrix構造体を定義します．
#code(C){{
struct Matrix
{
int width;
int height;
float *elements;
};
}}
構造体の定義はC++風で問題ないようです(いちいちtypedefつかわなくてよい)．
widthとheightが行列の列数と行数で，elementsが行列の要素を格納する配列です．
上記の行列Aのサイズだと，width=M，height=Nです．
ホスト側の行列の定義は以下のようにしました．
ただし，行列のサイズは適当です(面倒が少ない用にブロックサイズの倍数にはしてますが)．
#code(C){{
Matrix hA, hB, hC;
hA.height = hC.height = 3*BLOCK_SIZE;
hA.width = hB.height = 4*BLOCK_SIZE;
hB.width = hC.width = 5*BLOCK_SIZE;
hA.elements = new float[hA.width*hA.height];
hB.elements = new float[hB.width*hB.height];
hC.elements = new float[hC.width*hC.height];
RandomInit(hA.elements, hA.width*hA.height);
RandomInit(hB.elements, hB.width*hB.height);
}}
RandomInitは[[CUDAで行列演算：加減算]]のものと同じで，テスト用に乱数を行列の要素に代入する関数です．
デバイス側のメモリ確保とデータ転送も前回と同様に
#code(C){{
Matrix dA, dB, dC;
dA.width = hA.width; dA.height = hA.height;
dB.width = hB.width; dB.height = hB.height;
dC.width = hC.width; dC.height = hC.height;
int size;
// デバイスメモリの確保とホストからの転送
size = dA.width*dA.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dA.elements, size));
cutilSafeCall(cudaMemcpy(dA.elements, hA.elements, size, cudaMemcpyHostToDevice));
size = dB.width*dB.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dB.elements, size));
cutilSafeCall(cudaMemcpy(dB.elements, hB.elements, size, cudaMemcpyHostToDevice));
size = dC.width*dC.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dC.elements, size));
}}
**カーネル [#g9ea828d]
カーネルの実行は前回と同様です．
#code(C){{
dim3 block(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid((dC.width+block.x-1)/block.x, (dC.height+block.y-1)/block.y);
matrixMul<<< grid, block >>>(dA, dB, dC);
// カーネル実行エラーのチェック
cutilCheckMsg("Kernel execution failed");
}}
カーネル名はmatrixMulとしています．
カーネルは行列Cのサイズと同じ数実行され，各カーネルスレッドは行列Cのそれぞれの要素を計算します．
カーネル関数の実装は，
#code(C){{
__global__
void matrixMul(Matrix A, Matrix B, Matrix C)
{
int row = blockIdx.y*blockDim.y+threadIdx.y;
int col = blockIdx.x*blockDim.x+threadIdx.x;
if(row < C.height && col < C.width){
float x = 0.0f;
for(int k = 0; k < A.width; ++k){
x += A.elements[row*A.width+k]*B.elements[k*B.width+col];
}
C.elements[row*C.width+col] = x;
}
}
}}
行列Cの各要素をAのrow行とBのcol列からfor文を用いて計算しています．
**ホストへの結果の転送と後片付け [#z394441a]
結果をホストに転送し，確保したメモリを解放します．
#code(C){{
// デバイスからホストへ結果を転送
size = dC.width*dC.height*sizeof(float);
cutilSafeCall(cudaMemcpy(hC.elements, dC.elements, size, cudaMemcpyDeviceToHost));
// デバイスメモリ解放
cutilSafeCall(cudaFree(dA.elements));
cutilSafeCall(cudaFree(dB.elements));
cutilSafeCall(cudaFree(dC.elements));
// ホストメモリ解放
delete [] hA.elements;
delete [] hB.elements;
delete [] hC.elements;
}}
**計算時間 [#od18e551]
CUDAで実装したコードの計算時間を計算して，CPU実装と比較してみました．実行環境は以下です．
-GPU : GeForce GTX285
-CPU : Core 2 Duo 3.16GHz
CPUの方が少し古くGPU有利ですが参考データ程度にということで．
時間の計測にはSDKの関数を用いました．
#code(C){{
unsigned int timer = 0;
cutCreateTimer(&timer);
cutStartTimer(timer);
// 処理
cutilSafeCall(cudaThreadSynchronize());
cutStopTimer(timer);
printf("Processing time: %f (ms) \n", cutGetTimerValue(timer));
cutDeleteTimer(timer);
}
cudaThreadSynchronise()はデバイス側の処理が全て終わるまで待つ関数です．
CPUでの処理用に以下のような関数を用意しました．
#code(C){{
void matrixMulByCPU(const Matrix A, const Matrix B, Matrix C)
{
for(int i = 0; i < C.height; ++i){
for(int j = 0; j < C.width; ++j){
float cval = 0.0f;
for(int k = 0; k < A.width; ++k){
cval += A.elements[i*A.width+k]*B.elements[k*B.width+j];
}
C.elements[i*C.width+j] = cval;
}
}
}
}}
行列の行と列のサイズを同じとし，その大きさはブロックサイズを16×16としてその1倍,10倍,100倍として計算してみました．
計算時間は以下の表です．
| |16|160|1600|
|GPU|0.085|0.42|398|
|CPU|0.010|12.58|21245|
時間の単位は ms(ミリ秒)です．一回計測しただけなので正確性は低いです．
GPUでの計算部分のみ(カーネル呼出し部のみ)を計測していて，デバイスメモリの確保・解放，デバイスメモリへの転送などは含めていませんが，
それでも行列のサイズが十分大きいときにはGPUでの実装が有用であると考えられます．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
