CUDA第一個程式
- CUDA 目前有兩種不同的 API：Runtime API 和 Driver API，兩種 API 各有其適用的範圍。
CUDA 的初始化
- 要使用 runtime API 的時候，需要 include cuda_runtime.h。所以，在程式的最前面，加上:
#include <stdio.h>
#include <cuda_runtime.h>
bool InitCUDA()
{
int count;
/* 取得支援 CUDA 的裝置的數目，如果系統上沒有支援 CUDA 的裝置，則它會傳回 1，
而 device 0 會是一個模擬的裝置，但不支援 CUDA 1.0 以上的功能。
*/
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
/* 要確定系統上是否有支援 CUDA 的裝置，需要對每個 device 呼叫 cudaGetDeviceProperties 函式，
取得裝置的各項資料，並判斷裝置支援的 CUDA 版本（prop.major 和 prop.minor 分別代表裝置支援
的版本號碼，例如 1.0 則 prop.major 為 1 而 prop.minor 為 0）
透過 cudaGetDeviceProperties 函式可以取得許多資料，除了裝置支援的 CUDA 版本之外，還有裝置的名稱、
記憶體的大小、最大的 thread 數目、執行單元的時脈等等。
*/
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
/* 在找到支援 CUDA 1.0 以上的裝置之後，就可以呼叫 cudaSetDevice 函式，把它設為目前要使用的裝置。 */
cudaSetDevice(i);
return true;
}
int main() {
// 在 main 函式中我們直接呼叫剛才的 InitCUDA 函式，並顯示適當的訊息
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
return 0;
}
nvcc 是 CUDA 的 compile 工具，它會將 .cu 檔拆解出在 GPU 上執行的部份，及在 host 上執行的部份，並呼叫適當的程式進行 compile 動作。在 GPU 執行的部份會透過 NVIDIA 提供的 compiler 編譯成中介碼，而 host 執行的部份則會透過系統上的 C++ compiler 編譯（在 Windows 上使用 Visual C++ 而在 Linux 上使用 gcc）。
- nvcc -o first_cuda first_cuda.cu
編譯後的程式，執行時如果系統上有支援 CUDA 的裝置，應該會顯示 CUDA initialized. 的訊息，否則會顯示相關的錯誤訊息。
利用CUDA進行計算
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#define DATA_SIZE 1048576
int data[DATA_SIZE];
bool InitCUDA(){
/* 初始化CUDA裝置 */
int count;
/* 取得支援 CUDA 的裝置的數目，如果系統上沒有支援 CUDA 的裝置，則它會傳回 1，
而 device 0 會是一個模擬的裝置，但不支援 CUDA 1.0 以上的功能。
*/
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
/* 要確定系統上是否有支援 CUDA 的裝置，需要對每個 device 呼叫 cudaGetDeviceProperties 函式，
取得裝置的各項資料，並判斷裝置支援的 CUDA 版本（prop.major 和 prop.minor 分別代表裝置支援
的版本號碼，例如 1.0 則 prop.major 為 1 而 prop.minor 為 0）
透過 cudaGetDeviceProperties 函式可以取得許多資料，除了裝置支援的 CUDA 版本之外，還有裝置的名稱、
記憶體的大小、最大的 thread 數目、執行單元的時脈等等。
*/
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
/* 在找到支援 CUDA 1.0 以上的裝置之後，就可以呼叫 cudaSetDevice 函式，把它設為目前要使用的裝置。 */
cudaSetDevice(i);
return true;
}
void GenerateNumbers(int *number, int size){
for(int i = 0; i < size; i++) {
number[i] = rand() % 10;
}
}
__global__ static void sumOfSquares(int *num, int* result, clock_t* time) {
/* 在函式前面加上 __global__ 表示這個函式是要在顯示晶片上執行的。
在顯示晶片上執行的程式有一些限制，例如它不能有傳回值。
clock 函式，可以取得目前的 timestamp，很適合用來判斷一段程式執行所花費
的時間（單位為 GPU 執行單元的時脈）
*/
int sum = 0;
int i;
clock_t start = clock();
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
*time = clock() - start;
}
int main() {
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
GenerateNumbers(data, DATA_SIZE);
int* gpudata, *result;
clock_t* time;
/* 呼叫 cudaMalloc 取得一塊顯示記憶體（result 則是用來存取計算結果，在稍後會用到），
並透過 cudaMemcpy 將產生的亂數複製到顯示記憶體中。
*/
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMalloc((void**) &time, sizeof(clock_t));
/* cudaMalloc 和 cudaMemcpy 的用法和一般的 malloc 及 memcpy 類似，不過 cudaMemcpy 則多出一個參數，
指示複製記憶體的方向。 在這裡因為是從主記憶體複製到顯示記憶體，所以使用 cudaMemcpyHostToDevice。
*/
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
/* CUDA 中，要執行一個函式，使用以下的語法:
函式名稱<<<block 數目, thread 數目, shared memory 大小>>>(參數...);
這個程式只使用一個 thread，所以 block 數目、thread 數目都是 1。我們也沒有使用到任何 shared memory，
所以設為 0。
*/
sumOfSquares<<<1, 1, 0>>>(gpudata, result, time);
int sum;
clock_t time_used;
cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("cuda sum: %d, time %d\n", sum, time_used);
sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("CPU sum: %d\n", sum);
return 0;
}
增加thread加速程式
我們做了一個計算一大堆數字的平方和的程式。但這個程式的執行效率並不理想。實際上來說，如果只是要做計算平方和的動作，用 CPU 做會比用 GPU 快得多。這是因為平方和的計算並不需要太多運算能力，所以幾乎都是被記憶體頻寬所限制。因此，光是把資料複製到顯示記憶體上的這個動作，所需要的時間，可能已經和直接在 CPU 上進行計算差不多了。
不過，如果進行平方和的計算，只是一個更複雜的計算過程的一部份的話，那麼當然在 GPU 上計算還是有它的好處的。而且，如果資料已經在顯示記憶體上（例如在 GPU 上透過某種演算法產生），那麼，使用 GPU 進行這樣的運算，還是會比較快的。
剛才也提到了，由於這個計算的主要瓶頸是記憶體頻寬，所以，理論上顯示卡的記憶體頻寬是相當大的。這裡我們就來看看，倒底我們的第一個程式，能利用到多少記憶體頻寬。
我們的第一個程式，並沒有利用到任何平行化的功能。整個程式只有一個 thread。在 CUDA 中，一般的資料複製到的顯示記憶體的部份，稱為 global memory。這些記憶體是沒有 cache 的，而且，存取 global memory 所需要的時間（即 latency）是非常長的，通常是數百個 cycles。由於我們的程式只有一個 thread，所以每次它讀取 global memory 的內容，就要等到實際讀取到資料、累加到 sum 之後，才能進行下一步。這就是為什麼它的表現會這麼的差。
由於 global memory 並沒有 cache，所以要避開巨大的 latency 的方法，就是要利用大量的 threads。假設現在有大量的 threads 在同時執行，那麼當一個 thread 讀取記憶體，開始等待結果的時候，GPU 就可以立刻切換到下一個 thread，並讀取下一個記憶體位置。因此，理想上當 thread 的數目夠多的時候，就可以完全把 global memory 的巨大 latency 隱藏起來了。
要怎麼把計算平方和的程式平行化呢？最簡單的方法，似乎就是把數字分成若干組，把各組數字分別計算平方和後，最後再把每組的和加總起來就可以了。一開始，我們可以把最後加總的動作，由 CPU 來進行。
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#define DATA_SIZE 1048576
#define THREAD_NUM 256
int data[DATA_SIZE];
bool InitCUDA(){
/* 初始化CUDA裝置 */
int count;
/* 取得支援 CUDA 的裝置的數目，如果系統上沒有支援 CUDA 的裝置，則它會傳回 1，
而 device 0 會是一個模擬的裝置，但不支援 CUDA 1.0 以上的功能。
*/
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
/* 要確定系統上是否有支援 CUDA 的裝置，需要對每個 device 呼叫 cudaGetDeviceProperties 函式，
取得裝置的各項資料，並判斷裝置支援的 CUDA 版本（prop.major 和 prop.minor 分別代表裝置支援
的版本號碼，例如 1.0 則 prop.major 為 1 而 prop.minor 為 0）
透過 cudaGetDeviceProperties 函式可以取得許多資料，除了裝置支援的 CUDA 版本之外，還有裝置的名稱、
記憶體的大小、最大的 thread 數目、執行單元的時脈等等。
*/
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
/* 在找到支援 CUDA 1.0 以上的裝置之後，就可以呼叫 cudaSetDevice 函式，把它設為目前要使用的裝置。 */
cudaSetDevice(i);
return true;
}
void GenerateNumbers(int *number, int size){
for(int i = 0; i < size; i++) { number[i] = rand() % 10; } } __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { /* 在函式前面加上 __global__ 表示這個函式是要在顯示晶片上執行的。 在顯示晶片上執行的程式有一些限制，例如它不能有傳回值。
clock 函式，可以取得目前的 timestamp，很適合用來判斷一段程式執行所花費
的時間（單位為 GPU 執行單元的時脈）
程式裡的 threadIdx 是 CUDA 的一個內建的變數，表示目前的 thread 是第幾個 thread（由 0 開始計算）。以我們的例子來說，會有 256 個 threads，所以同時會有 256 個 sumOfSquares 函式在執行，但每一個的 threadIdx.x 則分別會是 0 ~ 255。利用這個變數，我們就可以讓每一份函式執行時，對整個資料不同的部份計算平方和。
*/
const int tid = threadIdx.x;
const int size = DATA_SIZE / THREAD_NUM;
int sum = 0;
int i;
clock_t start;
/* 另外，我們也讓計算時間的動作，只在 thread 0（即 threadIdx.x = 0 的時候）進行。 */
if(tid == 0)
start = clock();
for(i = tid * size; i < (tid + 1) * size; i++) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0)
*time = clock() - start;
}
int main() {
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
GenerateNumbers(data, DATA_SIZE);
int* gpudata, *result;
clock_t* time;
/* 呼叫 cudaMalloc 取得一塊顯示記憶體（result 則是用來存取計算結果，在稍後會用到），
並透過 cudaMemcpy 將產生的亂數複製到顯示記憶體中。
*/
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
/* 同樣的，由於會有 256 個計算結果，所以原來存放 result 的記憶體位置也要擴大。*/
cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
/* cudaMalloc 和 cudaMemcpy 的用法和一般的 malloc 及 memcpy 類似，不過 cudaMemcpy 則多出一個參數，
指示複製記憶體的方向。 在這裡因為是從主記憶體複製到顯示記憶體，所以使用 cudaMemcpyHostToDevice。
*/
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
/* CUDA 中，要執行一個函式，使用以下的語法:
函式名稱<<<block 數目, thread 數目, shared memory 大小>>>(參數...);
*/
sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
int sum[THREAD_NUM];
clock_t time_used;
cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
int final_sum = 0;
for(int i = 0; i < THREAD_NUM; i++) {
final_sum += sum[i];
}
printf("cuda sum: %d, time %d\n", final_sum, time_used);
final_sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
final_sum += data[i] * data[i];
}
printf("CPU sum: %d\n", final_sum);
return 0;
}
- 比較結果：
- first_cuda2.cu: 69262163 cycles
- first_cuda3.cu: 2586170 cycles
- 大約增進26倍的效能
考慮記憶體的存取模式
顯示卡上的記憶體是 DRAM，因此最有效率的存取方式，是以連續的方式存取。前面的程式，雖然看起來是連續存取記憶體位置（每個 thread 對一塊連續的數字計算平方和），但是我們要考慮到實際上 thread 的執行方式。
前面提過，當一個 thread 在等待記憶體的資料時，GPU 會切換到下一個 thread。也就是說，實際上執行的順序是類似
- thread 0 -> thread 1 -> thread 2 -> ...
因此，在同一個 thread 中連續存取記憶體，在實際執行時反而不是連續了。要讓實際執行結果是連續的存取，我們應該要讓 thread 0 讀取第一個數字，thread 1 讀取第二個數字...依此類推。所以，我們可以把 kernel 程式改寫如下。
- 僅僅是這樣簡單的修改，實際執行的效率就有很大的差別。
- 理論上 256 個 threads 最多只能隱藏 256 cycles 的 latency。但是 GPU 存取 global memory 時的 latency 可能高達 500 cycles 以上。如果增加 thread 數目，就可以看到更好的效率。但若 thread 數目增加太多，那麼在 CPU 端要做的最後加總工作也會變多。 first_cuda4.cu
__global__ static void sumOfSquares(int *num, int* result, clock_t* time) {
const int tid = threadIdx.x;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
- 比較結果：
- first_cuda2.cu: 69262163 cycles
- first_cuda3.cu: 2586170 cycles
- first_cuda4.cu: 613702 cycles
增加block數量
到目前為止，我們都只使用一個 block。究竟 block 是什麼呢？
在 CUDA 中，thread 是可以分組的，也就是 block。一個 block 中的 thread，具有一個共用的 shared memory，也可以進行同步工作。不同 block 之間的 thread 則不行。在我們的程式中，其實不太需要進行 thread 的同步動作，因此我們可以使用多個 block 來進一步增加 thread 的數目。
首先，在 #define DATA_SIZE 的地方，改成如下：
#define DATA_SIZE 1048576 #define BLOCK_NUM 32 #define THREAD_NUM 256這表示我們會建立 32 個 blocks，每個 blocks 有 256 個 threads，總共有 32*256 = 8192 個 threads。
接著，我們把 kernel 部份改成如下：
__global__ static void sumOfSquares(int *num, int* result, clock_t* time) {
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int sum = 0;
int i;
if(tid == 0) time[bid] = clock();
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
sum += num[i] * num[i];
}
result[bid * THREAD_NUM + tid] = sum;
if(tid == 0)
time[bid + BLOCK_NUM] = clock();
}
比較結果：
- first_cuda2.cu: 69262163 cycles
- first_cuda3.cu: 2586170 cycles
- first_cuda4.cu: 613702 cycles
- first_cuda5.cu: 28132 cycles
這個版本的程式，執行的時間減少很多，過，它在 CPU 上執行的部份，需要的時間加長了（因為 CPU 現在需要加總 8192 個數字）。為了避免這個問題，我們可以讓每個 block 把自己的每個 thread 的計算結果進行加總。
Thread的同步
一個 block 內的 thread 可以有共用的記憶體，也可以進行同步。我們可以利用這一點，讓每個 block 內的所有 thread 把自己計算的結果加總起來。把 kernel 改成如下： first_cuda6.cu
__global__ static void sumOfSquares(int *num, int* result, clock_t* time) {
/* 利用 __shared__ 宣告的變數表示這是 shared memory，是一個 block 中每個 thread
都共用的記憶體。它會使用在 GPU 上的記憶體，所以存取的速度相當快，不需要擔心 latency 的問題。
*/
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
/* _syncthreads() 是一個 CUDA 的內建函式，表示 block 中所有的 thread 都要同步到這個點，才能繼續執行。
在我們的例子中，由於之後要把所有 thread 計算的結果進行加總，所以我們需要確定每個 thread 都已經把結果
寫到 shared[tid] 裡面了。
*/
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
result[bid] = shared[0];
}
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
- 可以注意到，現在 CPU 只需要加總 BLOCK_NUM 也就是 32 個數字就可以了。
- 不過，這個程式由於在 GPU 上多做了一些動作，所以它的效率會比較差一些。
當然，效率會變差的一個原因是，在這一版的程式中，最後加總的工作，只由每個 block 的 thread 0 來進行，但這並不是最有效率的方法。理論上，把 256 個數字加總的動作，是可以平行化的。最常見的方法，是透過樹狀的加法first_cuda7.cu
比較結果：
- first_cuda2.cu: 69262163 cycles
- first_cuda3.cu: 2586170 cycles
- first_cuda4.cu: 613702 cycles
- first_cuda5.cu: 28132 cycles
- first_cuda6.cu: 32270 cycles
- first_cuda7.cu: 29328 cycles