CUDA 程式設計(5) -- 第一支程式 (向量加法) - 顯卡

※ 第五章 第一支程式 (向量加法)

這單元以向量加法為例, 示範 CUDA 的寫作, 總共有四個版本, 使用隨機向量,
測試準確度和效能, 其中準確度我們拿另一個 host 的版本計算結果做為比較,
效能我們測量 kernel 和 host 函數執行數次的平均時間來比較, 分成以下四個範例
(ps. 副檔名記得要 .cu 而非 .cpp)

(1) 單一區塊, 單一執行緒
可以用來測試單一串流處理器 (SP, stream processor) 的效能, 並讓我們容易
和傳統 C 接軌, 因為只使用單一處理器, 程式碼完全相同.

(2) 單一區塊, 多執行緒
可以用來測試單一多處理器 (MP, multiprocessor) 的效能, 一般而言單一區塊
不論包含多少個執行緒, 都使用一個 MP 來執行, 也就是 MP 裡的 8 個 SP
輪流在這些執行緒間切換, 硬體上以 warp (32 個執行緒)為單位進行群組切換,
當有多個區塊存在且資源足夠時(暫存器&共享記憶體), 可能會有多個區塊合併
到同一個 MP 執行.

(3) 多區塊, 多執行緒 (不使用迴圈, 用網格與區塊設定代替)
當迴圈數目小於 max(gridDim)*max(blockDim) 而且前後無相依性時, 可以
將迴圈打平, 使用網格與區塊設定代替, 每一個執行緒只計算一次, 整個核心
就是一個平行化的大迴圈, 從 compute 1.0 到 1.3, 這些數目沒什麼更動

每個網格最大的區塊數量 max(gridDim) 是 65,535
每個區塊最大的執行緒數量 max(blockDim)是 512

這種做法是為了避免迴圈拖慢效能, 因為 GPGPU 的分支指令(branch, 條件判斷)
是比較弱的, 理由是其執行上是以 warp 為單位, 當 32 個執行緒條件不同時,
就會導致部份的執行緒停擺, 等到下一執行週期再執行, 這種變成循序執行的
動作稱為發散 (divergence), 會造成這段指令執行需要兩倍時間而拖慢效能.

這種去除條件判斷和迴圈的程式設計法我稱為「乾式(dry)」的設計法, 因為
缺乏流程控制 (flow control), 其缺點為在比較複雜的程式中容易失去彈性,
而且必需付出計算資料位址的額外成本 (每個執行緒都必需計算一次).

(4) 多區塊, 多執行緒 (使用迴圈)
一般而言, 迴圈會造成條件發散的只有最後一個, 也就是大部份的 warp 執行
仍是良好的, 所以我們也不用那麼緊張的將迴圈弄成乾式, 即使存在時常發生
發散的條件判斷也不用擔心, 只要把條件所附帶的程式碼儘量簡短就行了 (讓
發散變成兩倍的程式碼執行時間縮短, 額外的時間便能縮短), 融合條件判斷
與平行化可使程式兼具彈性與效能。


===========================================================
test1.cu (單一區塊, 單一執行緒)
===========================================================
以下是範例(1)的程式, 只使用到單一處理器, 核心 gpu_add() 和對照組的 host_add()
完全相同, 將它存成 test1.cu 後, 在 linux 下使用 nvcc test1.cu -O3 即可編譯,
Windows 上可以修改 sdk 的 template project 來編譯, 其行為如下圖

執行緒 (一個一個執行)
(0) --> 記憶體
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo


在 GTX 280 上執行的結果和 Intel Q9400 的比較

vector add N(1048576) elements, diff = 0
host time: 3.3 ms
gpu time: 505.8 ms

也就是單一核心 GPU 慢上 150 倍, 原因除了單核 GPU 本身計算力較弱外, 主要還是
單核心在記憶體 I/O 上的問題 (將在後續章節再討論).

===========================================================
test1.cu 程式碼
===========================================================
#include <cuda.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//----------------------------------------------
//向量加法的運算核心 (GPU) **函式前加 __global__ 即為核心, 核心只傳回 void**
__global__ void gpu_add(float* c, float* a, float* b, int n){
for(int k=0; k<n; k++){
c[k]=a[k]+b[k];
}
}

//----------------------------------------------
//向量加法的一般函式 (Host)
void host_add(float* c, float* a, float* b, int n){
for(int k=0; k<n; k++){
c[k]=a[k]+b[k];
}
}


//----------------------------------------------
//計算誤差用的函式
double diff(float* a, float* b, int n){
double s=0, r=0;
for(int k=0; k<n; k++){
double w=a[k]-b[k];
s+=w*w;
r+=a[k]*a[k];
}
return sqrt(s/r); //相對誤差
}

//----------------------------------------------
//時間函數 (傳回單位:千分之一秒)
double ms_time(){
return (double)clock()/CLOCKS_PER_SEC*1000.0;
}

//----------------------------------------------
//主程式
int main(){
//設定向量大小
int n=1024*1024;
int size=n*sizeof(float);

//網格與區塊設定
int grid=1; //gridDim (每個網格具有的區塊數)
int block=1; //blockDim (每個區塊具有的執行緒數)

//設定呼叫次數 (測量平均效能)
int loop=100;

//配置主機記憶體
float *a,*b,*c,*d;
a=(float*)malloc(size);
b=(float*)malloc(size);
c=(float*)malloc(size);
d=(float*)malloc(size);

//設定亂數的輸入向量
srand(time(0));
for(int k=0; k<n; k++){
a[k]=(float)rand()/RAND_MAX*2-1;
b[k]=(float)rand()/RAND_MAX*2-1;
}

//配置顯示卡記憶體
float *ga,*gb,*gc;
cudaMalloc((void**)&ga, size);
cudaMalloc((void**)&gb, size);
cudaMalloc((void**)&gc, size);

//載入向量 a,b 到顯示卡記憶體中
cudaMemcpy(ga, a, size, cudaMemcpyHostToDevice);
cudaMemcpy(gb, b, size, cudaMemcpyHostToDevice);

//---- part 1 : 測量精確度 --------

//呼叫核心來運算 (GPU)
gpu_add<<<grid, block>>>(gc, ga, gb, n);

//呼叫一般函數來運算 (Host)
host_add(d, a, b, n);

//把計算結果存回主機
cudaMemcpy(c, gc, size, cudaMemcpyDeviceToHost);

//比較兩者差異
printf("vector add N(%d) elements, diff = %g\n", n, diff(c,d,n));



//---- part 2 : 測量效能 --------

//測量 GPU 核心效能
double gpu_dt = ms_time();
for(int w=0; w<loop; w++){
gpu_add<<<grid, block>>>(gc, ga, gb, n);
cudaThreadSynchronize(); //避免核心執行不完全
}
gpu_dt = (ms_time()-gpu_dt)/loop; //平均時間


//測量 Host 函數效能
double host_dt = ms_time();
for(int w=0; w<loop; w++){
host_add(d, a, b, n);
}
host_dt = (ms_time()-host_dt)/loop; //平均時間


//輸出平均執行時間
printf("host time: %g ms\n", host_dt);
printf("gpu time: %g ms\n", gpu_dt);


//釋放主機記憶體
free(a);
free(b);
free(c);
free(d);

//釋放顯示卡記憶體
cudaFree(ga);
cudaFree(gb);
cudaFree(gc);

return 0;
}

===========================================================
test2.cu (單一區塊, 多執行緒)
===========================================================
以下是範例(2)的需要修改的部份, 其中我們用到兩個內建唯讀變數

threadIdx.x 區塊中的執行緒索引
blockDim.x 區塊中的執行緒數量

假設我們有 8 個執行緒, 其 id 為 0~7, 這些平行的執行緒可以想向成一排印章,
沿著記憶體一塊一塊的印過去, 每次步進 8 個單位, 如下圖所示


執行緒區塊
(01234567) --> 記憶體
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

(01234567) 步進 8 個單位的距離 (blockDim.x)
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

(01234567) 再步進 8 個單位
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

...

(01234567)
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

執行的結果為

vector add N(1048576) elements, diff = 0
host time: 3.7 ms
gpu time: 12.2 ms

還是比 CPU 慢上一些, 因為還有很多 MP 沒有發揮.

===========================================================
test2.cu 程式碼在 test1.cu 要修改的部份
===========================================================
//運算核心
__global__ void gpu_add(float* c, float* a, float* b, int n){
for(int k=threadIdx.x; k<n; k+=blockDim.x){
c[k]=a[k]+b[k];
}
}

int main(){
...
//網格與區塊設定
int grid=1; //只使用一個區塊
int block=512;
...
}


===========================================================
test3.cu (多區塊, 多執行緒) 不使用迴圈,用網格與區塊設定代替
===========================================================

再來是使用乾式設計法的範例(3), 先定出 global idx, 然後一次印完整個向量,
使用到三個內建唯讀變數

blockIdx.x 網格中的區塊 id
blockDim.x 區塊中的執行緒數量
threadIdx.x 區塊中的執行緒 id

其中全域索引可用下式定址

global idx = blockIdx.x*blockDim.x + threadIdx.x

行為如下

global idx
(012345678...............................................N)
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
記憶體

執行結果為
vector add N(1048576) elements, diff = 0
host time: 3.56 ms
gpu time: 0.12 ms

所以 GPU 的能力被釋放出來了, 比 CPU 快上 30 倍.

===========================================================
test3.cu 程式碼在 test1.cu 要修改的部份
===========================================================
__global__ void gpu_add(float* c, float* a, float* b, int n){
int j=blockIdx.x*blockDim.x+threadIdx.x;
c[j]=a[j]+b[j];
}

int main(){
...
//網格與區塊設定
int block=512;
int grid=n/block; //必需能被整除
...
}



===========================================================
test4.cu (多區塊, 多執行緒) 使用迴圈,較為一般性
===========================================================
再來是使用迴圈的範例(4), 先定出 global id, 但是它並不像範例(3)一樣包含整個
向量大小, 然後像範例(2) 一塊一塊印過去

行為如下


global id (M = gridDim*blockDim)
(01234567...M-1) --> 記憶體
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

(01234567...M-1) 步進 M 個單位
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

...
(01234567...M-1)
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo


執行結果為

vector add N(1048576) elements, diff = 0
host time: 3.64 ms
gpu time: 0.12 ms

和乾式的效能差不多, 但向量長度變得更有彈性, 可以和 blockDim & gridDim 無關.

===========================================================
test4.cu 程式碼在 test1.cu 要修改的部份
===========================================================
__global__ void gpu_add(float* c, float* a, float* b, int n){
int j=blockIdx.x*blockDim.x+threadIdx.x;
int m=gridDim.x*blockDim.x;
for(int k=j; k<n; k+=m){
c[k]=a[k]+b[k];
}
}

int main(){
...
//網格與區塊設定
int block=256;
int grid=30;
...
}


===========================================================



※ 使用到的 CUDA API ※

---------------------------------------------------------------
(1)顯示卡記憶體配置 cudaMalloc()

cudaError_t cudaMalloc(void** ptr, size_t count);

ptr 指向目的指位器之位址
count 欲配置的大小(單位 bytes)

傳回值 cudaError_t 是個 enum, 執行成功時傳回 0, 其它的錯誤代號可用
cudaGetErrorString() 來解譯.

---------------------------------------------------------------
(2)顯示卡記憶體釋放 cudaFree()

cudaError_t cudaFree(void* ptr);

ptr 指向欲釋放的位址 (device memory)


---------------------------------------------------------------
(3)記憶體拷貝 cudaMemcpy()

cudaError_t cudaMemcpy(void* dst, const void* src, size_t count,
enum cudaMemcpyKind kind);

dst 指向目的位址
src 指向來源位址
count 拷貝區塊大小 (單位 bytes)
kind 有四種拷貝流向
cudaMemcpyHostToHost 主機 -> 主機
cudaMemcpyHostToDevice 主機 -> 裝置
cudaMemcpyDeviceToHost 裝置 -> 主機
cudaMemcpyDeviceToDevice 裝置 -> 裝置

---------------------------------------------------------------
(4) 錯誤字串解譯 cudaGetErrorString()

const char* cudaGetErrorString(cudaError_t error);

傳回錯誤代號(error)所代表的字串


---------------------------------------------------------------

※名詞解釋※

(1) 執行纜 (warp): 包含 32 個執行緒的單元, 多處理器(MP) 的執行緒切換以此為單位.
(2) 乾式 (dry): 相對於 flow 而言, 缺乏流程控制 (flow control),使用上缺乏彈性.
(3) 發散 (divergence): 同一 warp 的執行緒條件不同時, 就會導致部份的執行緒停擺,
等到下一執行週期再循序執行.
(4) 全域索引 (global idx): 由 (blockIdx,threadIdx) 合成的整體索引.



--