• トップ
  • G-DEPについて
  • ご購入ガイド
  • サポート
  • お問い合わせ

G-DEPトップ  >  第9回 プログラミング演習 ~基礎編②~

第9回 プログラミング演習 ~基礎編②~

第9回 「プログラミング演習 ~基礎編②~」

<< 第8回   |   目次に戻る   |   第10回 >>

前回から時間が空いてしまいましたが、今回はもう少し発展させて行列演算を使った並列処理問題を考えてみます。

計算の内容としては正方行列A、B、Cがあり、C = A × B(乗算)という単純なものです。

例えば、2行3列のCの要素は、以下のように計算します。

A21 × B13 + A22 × B23 + A23 × B33 + … + A2j × Bi3

 

GPUの場合、行列Cの1要素に付き1スレッドを使い、それに関連するAとBの数値をメモリから読み込んで並列に計算します。

つまりは、上の図で言うとi × j 個のスレッドで計算することになります。これが、CPUでの計算になると、一つの要素を順番に計算していくことになります。

明らかにGPUで並列計算した方が、速そうですよね。

 

それでは実際にCPUとGPUでどれだけ違うのか、問題を通して見ていきたいと思います。

(問題は第4回で説明したものを基に解説を加えて再構成しています。)

 

CPUによる計算


まずはCPUバージョンを見ていきましょう。

上でも述べましたが、GPU計算と異なるのは単体のCPUで全ての要素について順番に計算を行っているということです。

// CPU Matrix

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>

#define N 1024 // 正方行列のサイズを指定(N×N)

int main( int argc, char** argv)
{
    int col, row, scan;

    int* matA;
    int* matB;
    int* matC;

    // 行列A~Cのメモリ確保
    matA = (int*)malloc(sizeof(int) * N * N);
    matB = (int*)malloc(sizeof(int) * N * N);
    matC = (int*)malloc(sizeof(int) * N * N);

    // 行列A、Bの各要素に乱数を入れる
    // 行列Cは初期値0

    for (col = 0; col < N; col++){
        for (row = 0; row < N; row++){
            matA[col * N + row] = rand() % (N * N);
            matB[col * N + row] = rand() % (N * N);
            matC[col * N + row] = 0;
        }
    }

    // 行列の計算C = A + B
    for (col = 0; col < N; col++) {
        for (row = 0; row < N; row++) {
            for (scan = 0; scan < N; scan++) {
                matC[col * N + row] += matA[col * N + scan] * matB[scan * N + row];
            }
        }
    }

    // 各行列のメモリ解放
    free(matA);
    free(matB);
    free(matC);

    return 0;
}

 
1.ヘッダファイルの読み込み
 
2.「#define N 1024 // 正方行列のサイズを指定(N×N)
 今回取り扱う行列は正方行列です。その一辺のサイズを指定するもので、ここでは1024×1024の行列ということになります。
 
3.main()関数
 
4.「int col, row, scan」~「int* matC」
 変数の型を指定します。colはカラム(列)、row(行)、scanは行列の要素にアクセスして和を求めるための変数として使用しています。
 
5.行列A~Cのメモリ確保
 行列の各要素データを保管するためのメモリ領域を確保します。int型をN×N個用意しています。
 
6.行列変数の初期値を代入
 行列A、Bにはrand()によって乱数を代入し、行列Cには0を代入します。
 
【行列の各要素へのアクセスについて】
ある行列のi行目j列目の要素へアクセスしたい場合、以下の指定で読み込むことが出来ます。
 行列名[ i * N + j ]
 
7.行列の計算
 Cのある要素についてその行と列に対応する行列A、Bの要素をscanという変数によって読み込み計算をします。一つの要素について終わったら今度は同じ行で列を変え、1行終わったら次の行で同じ処理を繰り返します。
 
8.各行列のメモリ解放
 計算が終わったらメモリを解放します。

 

ここまでが一連の流れになります。

まずは行列演算のプログラムについて見ましたが、これでは単に計算するだけなので結果の表示もどれだけ計算時間かかったのかも分からず、面白みに欠けます。

そこで今度はタイマーを付けて、実際どれくらいの時間がかかるのか調べてみます。赤文字で示している所が追加した文です。

【タイマー付】

// CPU Matrix

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <time.h>

#define N 1024 // 正方行列のサイズを指定(N×N)

int main( int argc, char** argv)
{
    int col, row, scan;

    int* matA;
    int* matB;
    int* matC;

    // 時間計測の関数セット
    time_t timeStart, timeStop;


    // 行列A~Cのメモリ確保
    matA = (int*)malloc(sizeof(int) * N * N);
    matB = (int*)malloc(sizeof(int) * N * N);
    matC = (int*)malloc(sizeof(int) * N * N);
    
    // 時間計測開始
    time(&timeStart);


    // 行列A、Bの各要素に乱数を入れる
    // 行列Cは初期値0
    for (col = 0; col < N; col++){
        for (row = 0; row < N; row++){
            matA[col * N + row] = rand() % (N * N);
            matB[col * N + row] = rand() % (N * N);
            matC[col * N + row] = 0;
        }
    }

    // 行列の計算C = A + B
    for (col = 0; col < N; col++) {
        for (row = 0; row < N; row++) {
            for (scan = 0; scan < N; scan++) {
                matC[col * N + row] += matA[col * N + scan] * matB[scan * N + row];
            }
        }
    }

    // 時間計測終了
    time(&timeStop);


    // 計算時間結果の表示
    printf("It takes %dsec\n", timeStop - timeStart);


    // 各行列のメモリ解放
    free(matA);
    free(matB);
    free(matC);

    return 0;
}

これらは他の計算プログラムでも使うことができますので、覚えておくと便利かもしれません。

ポイントは、メインの計算が始まる(初期化含め)前にtimeStartを入れ、計算が終了する( for が } で閉じた後)にtimeStopを入れることです。

計算時間はPCの性能にもよりますが、筆者のPCでは10秒程度でした。場合によっては1分近くかかる可能性もあります。

正方行列の一辺の要素数Nを変えると、全体の要素数が変わる分、計算時間も変わってきます。ただし、値を大きく変えてしまうと計算が中々終わらなかったりしてしまいますので、ご注意を。こちらは下の演習で触れていきます。

 

GPU(CUDA)による計算


では、今度はGPUバージョンを見ていきましょう。

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <cutil_inline.h>

#define N 1024 // 正方行列のサイズを指定(N×N)
#define BLOCK 16 // ブロックのサイズを指定

__global__ void
matrixMul(int* inMatA, int* inMatB, int* inMatC);

int main(int argc, char** argv){
    
    // 行列のサイズをバイト単位で算出
    int matrixSize = sizeof(unsigned int) * N * N;
    
    // ホスト側の行列変数設定
    int* hMatA;
    int* hMatB;
    int* hMatC;

    // 行列変数のメモリ確保
    hMatA = (int*)malloc(matrixSize);
    hMatB = (int*)malloc(matrixSize);
    
    // 初期値設定
    int col, row;
    for (col = 0; col < N; col++){
        for (row = 0; row < N; row++){
            hMatA[col * N + row] = rand() % (N * N);
            hMatB[col * N + row] = rand() % (N * N);
        }    
    }

    // デバイス側の行列変数設定
    int* dMatA;
    int* dMatB;
    int* dMatC;
    
    // デバイスメモリ領域の確保
    cudaMalloc((void**)&dMatA, matrixSize);
    cudaMalloc((void**)&dMatB, matrixSize);
    cudaMalloc((void**)&dMatC, matrixSize);

    // ホストからデバイスへの変数の受け渡し
    cudaMemcpy(dMatA, hMatA, matrixSize, cudaMemcpyHostToDevice);
    cudaMemcpy(dMatB, hMatB, matrixSize, cudaMemcpyHostToDevice);

    // ブロックサイズとグリッドサイズの設定
    dim3 block(BLOCK, BLOCK);
    dim3 grid( N / BLOCK, N / BLOCK);

    // カーネルの起動
    matrixMul<<<grid, block>>>(dMatA, dMatB, dMatC);
    cudaThreadSynchronize();

    // 結果の領域確保とデバイス側からのメモリ転送
    hMatC = (int*)malloc(matrixSize);
    cudaMemcpy(hMatC, dMatC, matrixSize, cudaMemcpyDeviceToHost);

    // ホスト・デバイスメモリの解放
    free(hMatA);
    free(hMatB);
    free(hMatC);
    cudaFree(dMatA);
    cudaFree(dMatB);
    cudaFree(dMatC);

    // 終了処理
    cudaThreadExit();
    cutilExit(argc, argv);
}
 
__global__ void
matrixMul(int* inMatA, int* inMatB, int* inMatC)
 {
  int col = blockIdx.x * blockDim.x + threadIdx.x;
  int row = blockIdx.y * blockDim.y + threadIdx.y;
  int scan;
  int target = 0;

  // 行列の演算を行う
  for (scan = 0; scan < N; scan++) {
      target += inMatA[col * N + scan] * inMatB[scan * N + row];
      __syncthreads();
  }

  inMatC[col * N + row] = target;
}
 
1.ヘッダファイルの読み込み
 
2.#define N 1024 // 正方行列のサイズを指定(N×N)
 
3.#define BLOCK 16 // ブロックのサイズを指定
 1ブロック中のスレッド数を16としています。
 
4.「__global__ void matrixMul(int* inMatA, int* inMatB, int* inMatC);」
 行列の演算を行う関数の定義を行います。
 
5.main()関数
 
6.host側の行列変数の型を指定
 
7.host側の行列変数のためのメモリ領域確保
 
 hMatA、hMatBのデータ保管のためのメモリ領域を確保します。
 
8.初期値設定
 hMatA、hMatBに初期値を代入します。
 
9.device側の行列変数の型を指定
 
10.device側の行列変数のためのメモリ領域確保
 
11.host側→device側への変数の値を転送
 device側で計算させるためにhost側の変数の値を送ります。 ※ cudaMemcpy( 転送先、転送元、サイズ、転送方向)
 
12.ブロックサイズとグリッドサイズの設定
 dim3型変数の定義と、ブロックおよびグリッドのサイズを決めています。今回の場合、1ブロック中に16×16=256スレッド、1グリッド中に(1024/16)×(1024/16)=4096ブロックを使用しています。全部で256×4096=1048576スレッドが動き、つまりこれは
1024×1024の行列Cの要素一つずつに1スレッドを割り当てていることになります。
 
13.カーネルの起動
 やっとここでdevice側で計算させるための関数を起動します。今まで設定した変数の値、ブロックサイズ、グリッドサイズの情報が指定されて計算されます。
 
cudaThreadSynchronize(); について】
 現在のスレッドが何番目のものであるのかを取得し、それを元に計算し、結果を行列Cに入れます。このとき各スレッドが並列に動作しているため、スレッドごとに同期を取っています。今回はある要素での計算結果が他の要素の値へ影響を与えることがないためなくても必要ではありませんが、今後必ず使っていくものなので取り入れました。
 
14.devicet側→host側への変数の値を転送
 カーネルで計算された結果をhost側の変数に返します。
 
15.メモリの解放
 計算が終了したため、メインメモリ・デバイスメモリ両方の領域を解放します。
 
16.終了処理

 

今度はタイマーを加えてみましょう。CUDAの場合、構文が異なります。

【タイマー付】

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <cutil_inline.h>

#define N 1024 // 正方行列のサイズを指定(N×N)
#define BLOCK 16 // ブロックのサイズを指定

__global__ void
matrixMul(int* inMatA, int* inMatB, int* inMatC);

int main(int argc, char** argv){
    
    // 行列のサイズをバイト単位で算出
    int matrixSize = sizeof(unsigned int) * N * N;
    
    // ホスト側の行列変数設定
    int* hMatA;
    int* hMatB;
    int* hMatC;

    // 行列変数のメモリ確保
    hMatA = (int*)malloc(matrixSize);
    hMatB = (int*)malloc(matrixSize);
    
    // 初期値設定
    int col, row;
    for (col = 0; col < N; col++){
        for (row = 0; row < N; row++){
            hMatA[col * N + row] = rand() % (N * N);
            hMatB[col * N + row] = rand() % (N * N);
        }    
    }

    // デバイス側の行列変数設定
    int* dMatA;
    int* dMatB;
    int* dMatC;
    
    // デバイスメモリ領域の確保
    cudaMalloc((void**)&dMatA, matrixSize);
    cudaMalloc((void**)&dMatB, matrixSize);
    cudaMalloc((void**)&dMatC, matrixSize);

    // タイマーを作成して計測開始
    unsigned int timer = 0;
    cutCreateTimer( &timer);
    cutStartTimer( timer);


    // ホストからデバイスへの変数の受け渡し
    cudaMemcpy(dMatA, hMatA, matrixSize, cudaMemcpyHostToDevice);
    cudaMemcpy(dMatB, hMatB, matrixSize, cudaMemcpyHostToDevice);

    // ブロックサイズとグリッドサイズの設定
    dim3 block(BLOCK, BLOCK);
    dim3 grid( N / BLOCK, N / BLOCK);

    // カーネルの起動
    matrixMul<<<grid, block>>>(dMatA, dMatB, dMatC);
    cudaThreadSynchronize();

    // 結果の領域確保とデバイス側からのメモリ転送
    hMatC = (int*)malloc(matrixSize);
    cudaMemcpy(hMatC, dMatC, matrixSize, cudaMemcpyDeviceToHost);

    // 時間計測終了
    CUT_SAFE_CALL( cutStopTimer( timer));


    // 計算時間結果の表示(ミリ秒であることに注意
    printf("It takes %f (msec)\n", cutGetTimerValue( timer));
    cutDeleteTimer( timer);


    // ホスト・デバイスメモリの解放
    free(hMatA);
    free(hMatB);
    free(hMatC);
    cudaFree(dMatA);
    cudaFree(dMatB);
    cudaFree(dMatC);

    // 終了処理
    cudaThreadExit();
    cutilExit(argc, argv);
}
 
__global__ void
matrixMul(int* inMatA, int* inMatB, int* inMatC){
  int col = blockIdx.x * blockDim.x + threadIdx.x;
  int row = blockIdx.y * blockDim.y + threadIdx.y;
  int scan;
  int target = 0;

  // 行列の演算を行う
  for (scan = 0; scan < N; scan++) {
      target += inMatA[col * N + scan] * inMatB[scan * N + row];
      __syncthreads();
  }

  inMatC[col * N + row] = target;
}

「数値をメインメモリ→デバイスメモリに転送、カーネル起動、計算、結果をデバイスメモリ→メインメモリに返す」までを計算時間としています。

今回使用したCUDAタイマーは単位がミリ秒(1/1000秒)である点に注意してください。

こちらだと、だいたい106ミリ秒で終わりました。その速さ、なんとCPUの約100倍です。すごいですね!

 

それでは、今回の演習です。

演習


行列サイズNを変えると、どれだけ計算時間が変わるでしょうか?

⇒結果はコチラ

 

 

 いかがでしたか?簡単な問題ではありましたが、同じ計算内容でも並列計算においてはCPUとGPUでの処理速度の速さが歴然としていることが分かって頂けたでしょうか。

次回以降はもっとGPUコンピューティングの面白みを知ってもらうため、もう少し踏み込んだスレッドの使い方やイメージ描写などにも触れていきたいと思います。

 

<< 第8回   |   目次に戻る   |   第10回 >>


 

(執筆 G-DEP Associate Research Engineer 東京大学大学院工学系研究科 岡安優)