CUDAプログラミング第一歩

# はじめに
本記事はCUDA初学者の書く備忘録のため、新しい情報は記載していません。

## 参考
- http://www.kumikomi.net/archives/2008/10/22gpu2.php?page=1
- https://www.nvidia.co.jp/docs/IO/51174/NVIDIA_CUDA_Programming_Guide_1.1_JPN.pdf
- https://www.nvidia.co.jp/docs/IO/59373/VolumeI.pdf 

## 実行環境
- Ubuntu 16.04
- CUDA 10.1
- GeForce GTX 960

## 処理時間比較
CPU(シングルコア)とGPUでの処理速度の違いと比較するため、`time`コマンドを使用し、CPU、GPUそれぞれの処理にかかる時間を測定した。実行したコードは後の項目で示す。

```bash
#CPU
計算結果=18448244448962797568.000000

real    511m55.772s
user    511m37.328s

#GPU
knakajima@itage:~/work/cuda_test$ time ./a.out 
計算結果=18448244448962797568.000000

real    0m19.579s
user    0m13.976s
```

# GPUサンプルコード
GPUを使用する行列の掛け算を実行するプログラムを以下に示す。
このプログラムでは、GPUのグローバルメモリ上で演算を行っている。

```c++:mul_cuda.cu
#include <stdio.h>

const int BLOCK = 16;
const int WIDTH = 8192;

void Host(double *a, double *b, double *c);
__global__ void Kernel1(double *a, double *b, double *c);
__global__ void Kernel2(double *a, double *b, double *c);
__global__ void Kernel3(double *a, double *b, double *c);

double h_a[WIDTH * WIDTH];
double h_b[WIDTH * WIDTH];
double h_c[WIDTH * WIDTH];

__global__ void Kernel1(double *A, double *B, double *C) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    double tmp = 0.0;
    for (int k = 0; k < WIDTH; k++) {
        int row = k + y * WIDTH;
        int col = x + k * WIDTH;
        tmp += A[row] * B[col];
    }
    C[x + y * WIDTH] = tmp;
}

int main() {
    int i;

    // GPU上にメモリを確保
    double *d_a, *d_b, *d_c;
    cudaMalloc((void **)&d_a, sizeof(double) * WIDTH * WIDTH);
    cudaMalloc((void **)&d_b, sizeof(double) * WIDTH * WIDTH);
    cudaMalloc((void **)&d_c, sizeof(double) * WIDTH * WIDTH);
    cudaMemset(d_c, 0, sizeof(double) * WIDTH * WIDTH);

    for (i = 0; i < WIDTH * WIDTH; i++) {
        h_a[i] = i;
        h_b[i] = i;
    }

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

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

    //計算結果を取得
    cudaMemcpy(h_c, d_c, sizeof(double) * WIDTH * WIDTH,
               cudaMemcpyDeviceToHost);
    printf("計算結果=%f\n", h_c[WIDTH * WIDTH - 1]);

    // GPU上のメモリを開放
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
}
```
## 専用組み込み関数について
`cudaMalloc()`でGPU上のメモリを確保する。
`cudaMemset()`では指定したメモリ領域を引数で指定したバイト値で埋める。
`cudaMemcpy()`はGPU、CPU間でメモリをコピーする。
`cudaFree()`でGPUメモリを開放する。

## 関数型修飾子
関数名の頭につける`__global__`などの修飾子を関数型修飾子という。

|修飾子|説明|
|:-:|:-:|
|global|カーネル関数であることを宣言 |
|device|カーネル関数内で使用する関数であることを宣言|
|host|ホストで実行する関数であることを宣言|
なお、`host`と`device`を両方適応することもできる。

##カーネル関数
GPUで行う処理は、カーネル関数に記述する。
カーネル関数の呼び出しは以下のように記述する。
`カーネル関数名<<<グリッド数(dim3) , スレッド数(dim3)>>>()`
[グリッドとスレッドなど、GPUのアーキテクチャに関してはこちらを参照のこと](https://qiita.com/Keech/items/381495ed90e012d69f1a#cuda%E3%81%AE%E3%82%A2%E3%83%BC%E3%82%AD%E3%83%86%E3%82%AF%E3%83%81%E3%83%A3)

# CPUサンプルコード
演算速度を比較するために、CPUで同様の処理を実行するプログラムを用意した。

```c++:mul_cpu.cpp
#include <stdio.h>

const int WIDTH = 8192;

double h_a[WIDTH * WIDTH];
double h_b[WIDTH * WIDTH];
double h_c[WIDTH * WIDTH];

void Host(double *A, double *B, double *C) {
    //ホスト側での行列の乗算
    int i, j, k;
    int row, col;
    double tmp;
    for (i = 0; i < WIDTH; i++) {
        for (j = 0; j < WIDTH; j++) {
            tmp = 0.0;
            for (k = 0; k < WIDTH; k++) {
                row = k + i * WIDTH;
                col = j + k * WIDTH;
                tmp += A[row] * B[col];
            }
            C[j + i * WIDTH] = tmp;
        }
    }
}

int main() {
    int i;
    for (i = 0; i < WIDTH * WIDTH; i++) {
        h_a[i] = i;
        h_b[i] = i;
    }

    Host(h_a, h_b, h_c);
    printf("計算結果=%f\n", h_c[WIDTH * WIDTH - 1]);
}
```