CUDAプログラミング

次に、CUDA SDKのmatrixMultという行列積のプログラムに沿って、CUDAのプログラミングを見て行こう。プログラミングガイドに例題が載っている行列積と基本的には同じプログラムであるが、CUDA SDKでは、通常実行用のGPUバイナリを生成する以外に、デバグ機能を組み込んだバイナリを生成したり、GPUのない環境でCPUで実行するエミュレーションバイナリを生成したりする機能があり、SDKの例題は、直接CUDAを呼び出すのではなく、cutil.hで定義された呼び出し形式で書かれているという点が異なっている。

このcutil.hは、例えば後述のCUDA_SAFE_CALLは、通常実行モードバイナリでは単純に引数として与えられたCUDA関数の呼び出しに置き換えるが、_DEBUGが定義されていると、引数の関数を呼び出し、その戻り値のエラーを検査してエラーの場合はエラーメッセージを出力するという一連の文に置き換えるというようなことを行う。

行列積の計算を行うmatrixMul.cuの最初の部分は、次のようになっている。

// Program main

int
main(int argc, char** argv)
{
    runTest(argc, argv);

    CUT_EXIT(argc, argv);
}

このmainのプログラムはrunTestという関数を呼び出してCUDAの実行を行い、最後にCUT_EXIT関数を呼び出して後始末をするだけのプログラムである。そして、runTestの最初の部分は、次のようになっている。

//! Run a simple test for CUDA

void
runTest(int argc, char** argv)
{
    CUT_DEVICE_INIT(argc, argv);

    // set seed for rand()
    srand(2006);

    // allocate host memory for matrices A and B
    unsigned int size_A = WA * HA;
    unsigned int mem_size_A = sizeof(float) * size_A;
    float* h_A = (float*) malloc(mem_size_A);
    unsigned int size_B = WB * HB;
    unsigned int mem_size_B = sizeof(float) * size_B;
    float* h_B = (float*) malloc(mem_size_B);

    // initialize host memory
    randomInit(h_A, size_A);
    randomInit(h_B, size_B);

まず、CUT_DEVICE_INIT関数を呼んでGPUの初期化を行う。WA、HAは行列Aの幅(列数)と高さ(行数)、WB、HBは行列Bの幅と高さであり、これらの行列のサイズを計算して、mallocでホストのメモリ領域を確保している。そして、randomInit関数は、これらの行列A、Bにランダムに生成した数字を入れて初期化する関数である。

次に、GPU側のデバイスメモリ(グローバルメモリ)を確保する。CUDA_SAFE_CALL関数は、GPU側で実行するCUDA関数の実行を指示する関数で、実態はcudaMalloc関数が実行され、最初の呼び出しでは、行列Aを格納するためのd_Aでポイントされるデバイスメモリ領域が確保される。続いて、同様にして行列Bを格納する領域を確保する。次に、CUDA_SAFE_CALL関数を使ってcudaMemcpy関数を実行させ、ホストメモリのh_Aでポイントされる領域から、GPUのデバイスメモリのd_Aでポイントされる領域にデータをコピーする。最後のパラメタのcudaMemcpyHostToDeviceは、転送方向がホストからデバイスであることを指定している。また、行列Bも同様にデバイスメモリに転送する。

    // allocate device memory
    float* d_A;
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_A, mem_size_A));
    float* d_B;
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_B, mem_size_B));

    // copy host memory to device
    CUDA_SAFE_CALL(cudaMemcpy(d_A, h_A, mem_size_A,
                              cudaMemcpyHostToDevice) );
    CUDA_SAFE_CALL(cudaMemcpy(d_B, h_B, mem_size_B,
                              cudaMemcpyHostToDevice) );

そして、行列Aと行列Bの積である行列Cを格納する領域をデバイスメモリに確保し、更に、ホストメモリにも行列Cを受け取る領域を確保する。

    // allocate device memory for result
    unsigned int size_C = WC * HC;
    unsigned int mem_size_C = sizeof(float) * size_C;
    float* d_C;
    CUDA_SAFE_CALL(cudaMalloc((void**) &d_C, mem_size_C));

    // allocate host memory for the result
    float* h_C = (float*) malloc(mem_size_C);

次に、実行時間の測定のためにcutCreateTimer関数でGPU側にタイマーを定義し、cutStartTimer関数でタイマーをスタートする。

    // create and start timer
    unsigned int timer = 0;
    CUT_SAFE_CALL(cutCreateTimer(&timer));
    CUT_SAFE_CALL(cutStartTimer(timer));

ここからがGPUコンピューティングの本番で、dim3 threadsは2次元のスレッドブロックを定義している。この定義のBLOCK_SIZEは16となっており、256スレッドからなるブロックを定義している。そして次のdim3 gridでグリッドに含まれる2次元のブロックを定義している。ここで、threads.xはブロックのX方向のスレッド数、threads.yはY方向のスレッド数である。従って、全体では行列Cの各要素に対するスレッドが作られることになる。後述の性能測定では、最大2048元の行列積を計算しており、何と、その場合には200万スレッドが実行されている。

そして、matrixMul<<< grid,threads >>>が、計算の本体であるmatrixMul関数を、前に定義したグリッド数、スレッド数で実行するというという構文である。なお、matrixMul関数の内容については後述する。