CUDAをはじめてみる（2）

まず、このへんから開発ツールを落としてくる。

http://www.nvidia.co.jp/object/cuda_get_jp.html

それから、libとincludeにパスを通す。
サンプルを参考にして、前にやった数値積分の問題をCUDAで計算する。

http://d.hatena.ne.jp/inamori/20080312/p1

フルソースは下に隠す。が、重要な部分はそんなに多くない。
まず、配列を用意して、初期化する。

    // 配列データ準備

    float   *h_data = (float *)malloc(nData * sizeof(float));

    for(unsigned int i = 0; i < nData; i++) {

        h_data[i] = i * (M_PI / 2 / nData);

    }

データをGPUに転送したあと、GPUを走らせる。

    integral_sin<<< grid, threads >>>(d_data, size_x, size_y, nRep, dx);

<<< grid, threads >>>の部分は固定と思われる。関数名と引数は適当でよい。そして、それにあわせてGPU用の関数を作る。

__global__ void integral_sin(float *data,

                int width, int height, int nRep, float dx) {

    unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;

    unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;

    if( (xIndex < width) && (yIndex < height)) {

        unsigned int index = yIndex * width + xIndex;

        float   x = data[index];

このようにして、先ほどの転送した配列からデータを取り出す。
そしてこの配列の個々の要素を並列に処理する。

        data[index] = result;

ここでは、もう一度同じ配列に戻している。そして、GPUからデータを転送する。
コードが書けたら、VC2008のコマンドプロンプトを立ち上げ、次のように打つ。

>nvcc sin.cu -lcutil32

a.exeという実行ファイルができるので、それを実行する。

データの大きさ（size_xとsize_y）を適当に変えて、CPUオンリーのときと処理時間を比較してみた。CPUはXeon 3.20GHz、グラフィックボードはNVIDIA GeForce 9800 GX2。

       256   512   1024

CPU  589ms 2357ms 9430ms

GPU 23.8ms 30.5ms 61.2ms

オーバーヘッドを除くと、GPUはCPUのだいたい230倍の速度になっている。

// sin.cu

#include 

#include 

#include 

#define _USE_MATH_DEFINES

#include 

#include 
#define BLOCK_DIM 16
__global__ void integral_sin(float *data,

                int width, int height, int nRep, float dx) {

    unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;

    unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;

    if( (xIndex < width) && (yIndex < height)) {

        unsigned int index = yIndex * width + xIndex;

        float   x = data[index];

        float   result = 0;

        for(int i = 0; i < nRep; i++) {

            result += sin(x);

            x += dx;

        }

        data[index] = result;

    }

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

    runTest( argc, argv);

    cutilExit(argc, argv);

}
void runTest( int argc, char** argv) {

    const unsigned int  size_x = 1024;

    const unsigned int  size_y = 1024;

    const unsigned int  nData = size_x * size_y;

    const unsigned int  mem_size = sizeof(float) * nData;

    const int           nRep = 256;

    const double        dx = M_PI / 2 / nData / nRep;

    

    unsigned int timer;

    cutCreateTimer(&timer);

    

    if( cutCheckCmdLineFlag(argc, (const char**)argv, "device") )

        cutilDeviceInit(argc, argv);

    else

        cudaSetDevice( cutGetMaxGflopsDeviceId() );

    

    // 配列データ準備

    float   *h_data = (float *)malloc(nData * sizeof(float));

    for(unsigned int i = 0; i < nData; i++) {

        h_data[i] = i * (M_PI / 2 / nData);

    }

    

    cutStartTimer(timer);

    

    // データをGPUに転送

    float   *d_data;

    cutilSafeCall(cudaMalloc( (void**)&d_data, mem_size));

    cutilSafeCall(cudaMemcpy(d_data, h_data, mem_size, cudaMemcpyHostToDevice));

    

    // setup execution parameters

    dim3 grid(size_x / BLOCK_DIM, size_y / BLOCK_DIM, 1);

    dim3 threads(BLOCK_DIM, BLOCK_DIM, 1);

    

    // GPUで計算実行

    integral_sin<<< grid, threads >>>(d_data, size_x, size_y, nRep, dx);

    cudaThreadSynchronize();

    

    cutilCheckMsg("Kernel execution failed");

    

    // データをGPUから転送

    cutilSafeCall(cudaMemcpy(h_data, d_data, mem_size, cudaMemcpyDeviceToHost));

    float optimizedTime = cutGetTimerValue(timer);

    cutStopTimer(timer);

    

    float   integralGPU = 0.0f;

    for(int i = 0; i < nData; i++) {

        integralGPU += h_data[i];

    }

    printf("time: %0.3f ms\n\n", optimizedTime);

    printf("GPU : %f\n", integralGPU * (M_PI / 2 / nData / nRep));

    

    // cleanup memory

    free(h_data);

    cutilSafeCall(cudaFree(d_data));

    cutilCheckError( cutDeleteTimer(timer));

    

    cudaThreadExit();

}

#include 

#define _USE_MATH_DEFINES

#include 

#include 
int main() {

    const unsigned int  size_x = 256;

    const unsigned int  size_y = 256;

    const unsigned int  nData = size_x * size_y;

    const unsigned int  mem_size = sizeof(float) * nData;

    const int           nRep = 256;

    const double        dx = M_PI / 2 / nData / nRep;

    

    float   *h_data = (float *)malloc(nData * sizeof(float));

    for(unsigned int i = 0; i < nData; i++) {

        h_data[i] = i * (M_PI / 2 / nData);

    }

    

    LARGE_INTEGER nFreq, nStart, nEnd;

    QueryPerformanceFrequency(&nFreq);

    QueryPerformanceCounter(&nStart);

    

    for(int i = 0; i < nData; i++) {

        float   x = h_data[i];

        float   result = 0;

        for(int j = 0; j < nRep; j++) {

            result += sin(x);

            x += dx;

        }

        h_data[i] = result;

    }

    

    QueryPerformanceCounter(&nEnd);

    

    float   integralGPU = 0.0f;

    for(int i = 0; i < nData; i++) {

        integralGPU += h_data[i];

    }

    

    printf("time: %0.3f ms\n\n",

            (nEnd.QuadPart - nStart.QuadPart) * 1e3 / nFreq.QuadPart);

    printf("CPU : %f\n", integralGPU * (M_PI / 2 / nData / nRep));

    

    free(h_data);

}