CUDA技術を利用したGPUコンピューティングの実際（後編） ―― FFTを利用した光波の伝播（フレネル回折）をGPUで高速計算

#include <stdio.h>
#include <conio.h>
#include <cutil.h> //CUDAのユーティリティ関連ヘッダ・ファイル

//マクロ宣言
#define BLOCK 16
#define WIDTH 512

//プロトタイプ宣言
void Host(float *a, float *b, float *c);
__global__ void Kernel1(float *a, float *b, float *c);
__global__ void Kernel2(float *a, float *b, float *c);
__global__ void Kernel3(float *a, float *b, float *c);
void SetTimer(unsigned int *t);
float EndTimer(unsigned int *t);

//行列データと計算結果格納用配列
float h_a[WIDTH*WIDTH];
float h_b[WIDTH*WIDTH];
float h_c[WIDTH*WIDTH];

//メイン関数
int main()
{
    int i;
    unsigned int timer;

    //GPUの初期化　←(1)
    CUT_DEVICE_INIT();

    //GPU上にメモリを確保　←(2)
    float *d_a, *d_b, *d_c;
    cudaMalloc((void**) &d_a, sizeof(float)*WIDTH*WIDTH);
    cudaMalloc((void**) &d_b, sizeof(float)*WIDTH*WIDTH);
    cudaMalloc((void**) &d_c, sizeof(float)*WIDTH*WIDTH);
    cudaMemset(d_c, 0, sizeof(float)*WIDTH*WIDTH);
    //d_cを0で初期化する

    //行列データ作成　←(3)
    for(i=0; i<WIDTH*WIDTH; i++)
    {
    h_a[i]=i;
    h_b[i]=i;
    }
    SetTimer(&timer);

    //変数をGPU上のメモリへコピー　←(4)
    cudaMemcpy(d_a, h_a, sizeof(float)*WIDTH*WIDTH,
         cudaMemcpyHostToDevice);
     cudaMemcpy(d_b, h_b, sizeof(float)*WIDTH*WIDTH,
         cudaMemcpyHostToDevice);

    //ブロックとグリッドの定義とカーネルの起動　←(5)
    dim3 grid(WIDTH/BLOCK, WIDTH/BLOCK, 1);
    dim3 threads(BLOCK, BLOCK, 1);
    Kernel1<<< grid, threads >>>(d_a, d_b, d_c);
    ↑(カーネル関数名)

    //計算結果を取得　←(6)
    cudaMemcpy(h_c, d_c,
    sizeof(float)*WIDTH*WIDTH,cudaMemcpyDeviceToHost);
    printf("計算時間=%f(ms)",EndTimer(&timer));
    printf(" 計算結果=%f\n",h_c[WIDTH*WIDTH-1]);

     //GPU上のメモリを開放　←(7)
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);

    //ホスト側での計算（比較用）
    SetTimer(&timer);
    Host(h_a, h_b, h_c);
    printf("ホスト計算時間=%f(ms)",EndTimer(&timer));
    printf("計算結果=%f\n",h_c[WIDTH*WIDTH-1]);

    getch();
}