cudaAddKernel

/* double型の配列A・B・Cに対して「C=A+B」を実行している(ベクタ演算) */
/*
Debug(Win32)で
CPU計算時間:4.35151e-005[s] -> 376.513[MFLOPS]
GPU計算時間:4e-005[s] -> 409.6[MFLOPS]
最大絶対誤差:0.0000000000000000
Release(x64)で
CPU計算時間:1.51892e-005[s] -> 1078.66[MFLOPS]
GPU計算時間:3.7824e-005[s] -> 433.164[MFLOPS]
最大絶対誤差:0.0000000000000000
*/

/*
XXX: GT610M の性能の推測

GT610M : CC=2.1, 48コア, 900MHz, 14.4Gbyte/s
GT430  : CC=2.1, 96コア, 1400MHz, 28.8Gbyte/s (手持ち:現在はスロットが足りないので使ってない)
CUDA-Z 0.7.181 では
  int32 = 89.4 Giop/s
  float = 178.3 Gflop/s
  double = 22.3 Gflop/s

GT610M では、22.3 * 48/96 * 900/1400 = 7.16 Gflop/s

さらにaddKernelはメモリ律速だし(メモリ読み出し2回、double加算1回、メモリ書き出し1回）。
手元のGTS450をNsightで実行したら、
memory load 31.76 GB/s
memory store 15.88 GB/s
Achived FLOPS 2.13 Gflop/s
GT610Mではloadのバンド幅がこれより少ないので、確実にメモリ律速。

*/

/* プリプロセッサ(CUDA実行用) */
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

/* プリプロセッサ(その他) */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <windows.h>

/* 実行用関数。返り値はGPU演算の実行可否(cudaSuccessと等しければOK) */
cudaError_t addWithCuda(double *c, double *a, double *b, unsigned int size);

/* カーネル関数(ポインタ相手に演算。ゆえに返り値なし) */
__global__ void addKernel(double *c, double *a, double *b, int size_x)
{
    //blockDim.xはブロックサイズ、blockIdx.xはブロック番号、threadIdx.xはスレッドID
    int x = blockDim.x * blockIdx.x + threadIdx.x;
    int y = blockDim.y * blockIdx.y + threadIdx.y;
    int i = y * size_x + x;

    //加算(注:Cにおいては配列≒ポインタ)
    c[i] = a[i] + b[i];
}

int main(void){
    /* 問題サイズ。冷静に考えるとC文法で#defineしないのは変だが、
       CUDAでは許されている模様 */
    const int arraySize = 16 * 1024 * 1024;

    /* 各配列を定義 */
    double *a = NULL, *b = NULL, *c = NULL, *c_ = NULL;
    /* 時間計測用 */
    LARGE_INTEGER Frequency, time_s, time_g;
    double elapsedTime;
    /* その他変数 */
    int i;
    cudaError_t cudaStatus;

    /* XXX: スタックは小さい(せいぜい2MB)だから、ヒープから確保する */
    a = (double*) malloc(sizeof(double) * arraySize);
    b = (double*) malloc(sizeof(double) * arraySize);
    c = (double*) calloc(arraySize, sizeof(double));
    c_ = (double*) calloc(arraySize, sizeof(double));

    /* 配列初期化。とりあえず[0,1)とした */
    srand((unsigned) time(NULL));
    for(i = 0; i < arraySize; i++){
        a[i] = (double)(rand()) / RAND_MAX;
        b[i] = (double)(rand()) / RAND_MAX;
    }

    /* CPUによる演算 */
    memset(&Frequency, 0x00, sizeof Frequency);
    memset(&time_s,    0x00, sizeof time_s);
    memset(&time_g,    0x00, sizeof time_g);
    QueryPerformanceFrequency(&Frequency);
    QueryPerformanceCounter(&time_s);
    for(i = 0; i < arraySize; i++){
        c[i] = a[i] + b[i];
    }
    QueryPerformanceCounter(&time_g);
    elapsedTime = (double)(time_g.QuadPart - time_s.QuadPart) / Frequency.QuadPart;
    printf("CPU計算時間:%3.9lf[s] -> %g[MFLOPS]\n", elapsedTime, 1.0 * arraySize / elapsedTime / 1000.0 / 1000.0);

    /* GPUを利用して並列演算 */
    cudaStatus = addWithCuda(c_, a, b, arraySize);
    if (cudaStatus != cudaSuccess) {    //失敗時の動作
        fprintf(stderr, "addWithCuda failed!");
        return 1;
    }

    /* 誤差計算(デバッグ用) */
    double max_diff = 0.0;
    for(i = 0; i < arraySize; i++){
        if(max_diff < fabs(c[i] - c_[i]))
            max_diff = fabs(c[i] - c_[i]);
    }

    /* 結果表示 */
    /*for(i = 0; i < arraySize; i++){
        printf("%d %f\t%f\t%f\t%f\n", i, a[i], b[i], c[i], c_[i]);
    }*/
    printf("最大絶対誤差:%1.16lf\n", max_diff);

    /* 終了用関数呼び出し */
    cudaStatus = cudaDeviceReset();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceReset failed!");
        return 1;
    }

    if (a) free(a);
    if (b) free(b);
    if (c) free(c);
    if (c_) free(c_);

    return 0;
}

/* 実行用関数。カーネル関数をここから呼び出す */
cudaError_t addWithCuda(double *c, double *a, double *b, unsigned int size)
{
    //カーネル関数に値を渡すためのポインタ
    double *dev_a = 0;
    double *dev_b = 0;
    double *dev_c = 0;
    //ステータス用変数
    cudaError_t cudaStatus;
    //時間計測用
    cudaEvent_t start, stop;
    float elapsedTime;          //経過時間[ms]

    /* 1次元配列を2次元(size_x * size_y)とみなして、それをさらにCUDA blockに分割: CC=2.1の制限: gridsize>=65536を超えるため */
    const dim3 blockSize(256, 1);
    const int size_x = blockSize.x * 1024;  // 1ブロックは最大1024スレッド
    const int size_y = (size + size_x - 1) / size_x;
    const dim3 gridSize((size_x + blockSize.x - 1)/ blockSize.x, (size_y + blockSize.y - 1)/ blockSize.y);

    printf("size: %d\n", size);
    printf("size_x,y: %d,%d\n", size_x, size_y);
    printf("blockSize: %d,%d\n", blockSize.x, blockSize.y);
    printf("gridSize: %d,%d\n", gridSize.x, gridSize.y);


    /* どのCUDAデバイスを利用するかを選択(今回は番号0のもの)
       ・cudaGetDeviceCount(int *count)で総数を取得できる
       ・cudaGetDevice(int *current_device)で現在の番号を取得
       ・cudaGetDeviceProperties(int *device, cdaDeiceProp*prop)で情報取得 */
    cudaStatus = cudaSetDevice(1);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");
        goto Error;
    }

    /* カーネル関数に渡せるようにメモリを確保 */
    cudaStatus = cudaMalloc((void**)&dev_c, size * sizeof(double));
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
        goto Error;
    }
    cudaStatus = cudaMalloc((void**)&dev_a, size * sizeof(double));
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
        goto Error;
    }
    cudaStatus = cudaMalloc((void**)&dev_b, size * sizeof(double));
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
        goto Error;
    }

    /* 入力用にデータをコピーする */
    cudaStatus = cudaMemcpy(dev_a, a, size * sizeof(double), cudaMemcpyHostToDevice);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
        goto Error;
    }
    cudaStatus = cudaMemcpy(dev_b, b, size * sizeof(double), cudaMemcpyHostToDevice);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
        goto Error;
    }
    cudaThreadSynchronize();

    /* 時間計測の準備 */
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start,0);

    /* 各エレメントに対しカーネル関数を実行(それぞれブロック数、スレッド数*/
    addKernel<<<gridSize, blockSize>>>(dev_c, dev_a, dev_b, size_x);

    /* 時間計測結果を取得 */
    cudaEventRecord(stop,0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime,start,stop);

    printf("GPU計算時間:%3.9lf[s] -> %g[MFLOPS]\n", elapsedTime / 1000, size / elapsedTime / 1000.0);

    /* エラーが起きていないかをチェック */
    cudaStatus = cudaGetLastError();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
        goto Error;
    }

    /* 演算が完了するまで待機 */
    cudaStatus = cudaDeviceSynchronize();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);
        goto Error;
    }

    /* 出力用にデータをコピーする */
    cudaStatus = cudaMemcpy(c, dev_c, size * sizeof(double), cudaMemcpyDeviceToHost);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
        goto Error;
    }

/* エラー発生時時はすぐにメモリを解放して終了 */
Error:
    cudaFree(dev_c);
    cudaFree(dev_a);
    cudaFree(dev_b);

    return cudaStatus;
}