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での実装が有用であると考えられます. |