![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
CUDAで行列演算 その2 乗算†CUDAで行列の乗算を行います．このページで示すコードは CUDA Programming Guideの 3.2.2 Shared Memoryとほとんど同じですのでそちらも参照してください． 行列の定義†行列の積は和や差と異なり，それぞれの行列のサイズが異なります． A(N×M) * B(M×L) = C(N×L) CUDAで行列演算：加減算ではすべての行列のサイズが等しかったので楽でしたが， 行列のサイズが異なると管理が面倒になるので，新しくMatrix構造体を定義します． struct Matrix
{
int width;
int height;
float *elements;
};
構造体の定義はC++風で問題ないようです(いちいちtypedefつかわなくてよい)． widthとheightが行列の列数と行数で，elementsが行列の要素を格納する配列です． 上記の行列Aのサイズだと，width=M，height=Nです． ホスト側の行列の定義は以下のようにしました． ただし，行列のサイズは適当です(面倒が少ない用にブロックサイズの倍数にはしてますが)． 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で行列演算：加減算のものと同じで，テスト用に乱数を行列の要素に代入する関数です． デバイス側のメモリ確保とデータ転送も前回と同様に 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)); カーネル†カーネルの実行は前回と同様です． 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のそれぞれの要素を計算します． カーネル関数の実装は， __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文を用いて計算しています． ホストへの結果の転送と後片付け†結果をホストに転送し，確保したメモリを解放します． // デバイスからホストへ結果を転送 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; 計算時間†CUDAで実装したコードの計算時間を計算して，CPU実装と比較してみました．実行環境は以下です．
時間の計測にはSDKの関数を用いました． 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倍として計算してみました． 計算時間は以下の表です．
時間の単位は ms(ミリ秒)です．一回計測しただけなので正確性は低いです． GPUでの計算部分のみ(カーネル呼出し部のみ)を計測していて，デバイスメモリの確保・解放，デバイスメモリへの転送などは含めていませんが， それでも行列のサイズが十分大きいときにはGPUでの実装が有用であると考えられます． |
