CUDAで行列演算 その2 乗算

CUDAで行列の乗算を行います.このページで示すコードは CUDA Programming Guideの 3.2.2 Shared Memoryとほとんど同じですのでそちらも参照してください.

行列の定義

行列の積は和や差と異なり,それぞれの行列のサイズが異なります.

A(N×M) * B(M×L) = C(N×L)

CUDAで行列演算:加減算ではすべての行列のサイズが等しかったので楽でしたが, 行列のサイズが異なると管理が面倒になるので,新しくMatrix構造体を定義します.

  1
  2
  3
  4
  5
  6
struct Matrix
{
    int width;
    int height;
    float *elements;
};

構造体の定義はC++風で問題ないようです(いちいちtypedefつかわなくてよい). widthとheightが行列の列数と行数で,elementsが行列の要素を格納する配列です. 上記の行列Aのサイズだと,width=M,height=Nです.

ホスト側の行列の定義は以下のようにしました. ただし,行列のサイズは適当です(面倒が少ない用にブロックサイズの倍数にはしてますが).

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
    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で行列演算:加減算のものと同じで,テスト用に乱数を行列の要素に代入する関数です.

デバイス側のメモリ確保とデータ転送も前回と同様に

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
    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));

カーネル

カーネルの実行は前回と同様です.

  1
  2
  3
  4
  5
  6
  7
    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のそれぞれの要素を計算します.

カーネル関数の実装は,

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
__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文を用いて計算しています.

ホストへの結果の転送と後片付け

結果をホストに転送し,確保したメモリを解放します.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
    // デバイスからホストへ結果を転送
    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実装と比較してみました.実行環境は以下です.

  • GPU : GeForce GTX285
  • CPU : Core 2 Duo 3.16GHz CPUの方が少し古くGPU有利ですが参考データ程度にということで.

時間の計測にはSDKの関数を用いました.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 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倍として計算してみました. 計算時間は以下の表です.

161601600
GPU0.0850.42398
CPU0.01012.5821245

時間の単位は ms(ミリ秒)です.一回計測しただけなので正確性は低いです. GPUでの計算部分のみ(カーネル呼出し部のみ)を計測していて,デバイスメモリの確保・解放,デバイスメモリへの転送などは含めていませんが, それでも行列のサイズが十分大きいときにはGPUでの実装が有用であると考えられます.


トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2011-10-27 (木) 15:09:13