CUDA 程式設計(11) -- 速成篇(中) - 顯卡

號外~~ BT 牌的國網光碟已經完成了, 感謝 b 君和 c 君幫忙 ^^
種子下載點 http://www.badongo.com/file/12156676
請大家幫忙傳播~~~

【修正】

(1) 上次忘了說 GRID 的大小雖然是 uint3 的結構，但只能使用 2D 而已
(其 z 成員只能是 1)，BLOCK 才能完整支援 3D 結構。

(2) API 部份 cudaThreadSynchronize() 是用來進行核心和主機程序的同步，
把它錯打成 __syncthreads() 的功能，大概是那天太累了 @_@
詳見第六招。

============================================================================
第五招 記憶體 (主機、裝置、共享、暫存器)
============================================================================
基本的記憶體部份，最重要的是區分主機、裝置、共享記憶體以及暫存器，硬體位置
和存取速度列於下表：

----------------------------------------------------------------
中文名稱 英文名稱 硬體位置 存取速度
----------------------------------------------------------------
主機記憶體 Host Memory PC 上的 DRAM 透過 PCI-E, 很慢
裝置記憶體 Device Memory 顯示卡上的 DRAM 400~600 cycles
共享記憶體 Shared Memory 顯示晶片 4 cycles issue
暫存器 Register 顯示晶片 立即
----------------------------------------------------------------

這些記憶體的大小、功能和使用方式如下

----------------------------------------------------------------
名稱 大小 功能 使用方式
----------------------------------------------------------------
主機記憶體 0.5~10GB 存放傳統的主機資料 透過 API
裝置記憶體 0.5~10GB 存放給 GPU 使用資料 kernel 直接存取
共享記憶體 16KB/BLOCK 執行緒之間交換資料 kernel 直接存取
暫存器 32~64KB/BLOCK 執行緒的區域變數 kernel 直接存取
----------------------------------------------------------------

相關細節如下：

(1) 【主機記憶體】我們都很熟了，就是傳統 C/C++ 使用的變數，可以透過 malloc()、
free()、new、delete 等來配置或釋放。

(2) 【裝置記憶體】的地位和主機記憶體很像，只是主機記憶體是對付 CPU，而
裝置記憶體是對付 GPGPU，兩者間的資料傳送要透過 CUDA 的 API 達成，
可以透過 cudaMalloc()、cudaFree() 等函式來配置或釋放 (詳見第二招)，
在 CUDA 中，這類的記憶體稱為【全域記憶體 (global memory)】。

(3) 【共享記憶體】是比較特殊的，在傳統的序列化程式設計裡沒有直接的對應，
初學者可能要花多一點的時間在這上面，它用在『區塊內執行緒間交換資料』，
只能在 kernel 中使用，宣告時使用 __shared__ 這個標籤，使用時必需同步，
以確保資料讀寫時序上的正確性，現階段在每個區塊上最大的容量為 16KB，
超過便無法執行或編譯，因為 on chip 的關係，它屬於快速記憶體。

(4) 【暫存器】kernel 中大部份的區域變數都是以暫存器的型式存放，不需做額外的
宣告手序，這些暫存器是區塊中執行緒共享的，也就是如果每個執行緒使用到 8 個
暫存器，呼叫這個 kernel 時區塊大小指定為 10，則整個區塊使用到 8*10=80 個
暫存器，若呼叫時指定區塊大小為 50，則整個區塊使用 8*50=400 個暫存器。

(5) 當使用過多的暫存器時 (>32~64KB/BLOCK，看是那個世代的 GPU)，系統會自動把
一些資料置換到全域記憶體中，導致執行緒變多，但效率反而變慢 (類似作業系統
虛擬記憶體的 swap)；另一個會引發這種 swap 的情況是在使用動態索引存取陣列，
因為此時需要陣列的順序性，而暫存器本身是沒有所謂的順序的，所以系統會自動
把陣列置於全域記憶體中，再按索引存取，這種情況建議使用共享記憶體手動避免。

(6) 【暫存器】和【共享記憶體】的使用量會限制執行緒的數目，在開發複雜 kernel 時
宜注意，可使用 nvcc --ptxas-options=-v 這個選項在設計時期監控，或使用
nvcc --maxrregcount 選項限制每個執行緒的暫存器使用量。


//----------------------------------------------------------------------------
//範例(5): 平滑處理 (使用相鄰的三點做加權平均,使資料變平滑)[081119-smooth.cu]
// 執行緒同步 __syncthreads() 和 cudaThreadSynchronize()，詳見第六招
//---------------------------------------------------------------------------
#include<stdio.h>
#include<time.h>
#include<cuda.h>

//設定區塊大小 (shared 版本會用到, 所以先宣告).
#define BLOCK 512

//--------------------------------------------------------
//(1) 對照組 (host 版).
//--------------------------------------------------------
void smooth_host(float* b, float* a, int n){
for(int k=1; k<n-1; k++){
b[k]=(a[k-1]+2*a[k]+a[k+1])*0.25;
}

//邊界為0
b[0]=(2*a[0]+a[1])*0.25;
b[n-1]=(a[n-2]+2*a[n-1])*0.25;
}

//--------------------------------------------------------
//(2) 裝置核心(global 版).
//--------------------------------------------------------
__global__ void smooth_global(float* b, float* a, int n){
int k = blockIdx.x*blockDim.x+threadIdx.x;

if(k==0){
b[k]=(2*a[0]+a[1])*0.25;
}
else if(k==n-1){
b[k]=(a[n-2]+2*a[n-1])*0.25;
}
else if(k<n){
b[k]=(a[k-1]+2*a[k]+a[k+1])*0.25;
}

}


//--------------------------------------------------------
//(3) 裝置核心(shared 版).
//--------------------------------------------------------
__global__ void smooth_shared(float* b, float* a, int n){
//----------------------------------------
//計算區塊的基底
//----------------------------------------
int base = blockIdx.x*blockDim.x;
int t = threadIdx.x;

//----------------------------------------
//宣告共享記憶體.
//----------------------------------------
__shared__ float s[BLOCK+2];

//----------------------------------------
//載入主要資料 s[1]~s[BLOCK]
//----------------------------------------
// s[0] <-- a[base-1] (左邊界)
// s[1] <-- a[base]
// s[2] <-- a[base+1]
// s[3] <-- a[base+2]
// ...
// s[BLOCK] <-- a[base+BLOCK-1]
// s[BLOCK+1] <-- a[base+BLOCK] (右邊界)
//----------------------------------------
if(base+t<n){
s[t+1]=a[base+t];
}

//----------------------------------------
//載入邊界資料 s[0] & s[BLOCK+1] (只用兩個執行緒處理)
//----------------------------------------
if(t==0){
//左邊界.
if(base==0){
s[0]=0;
}
else{
s[0]=a[base-1];
}
}

//*** 使用獨立的 warp 讓 branch 更快 ***
if(t==32){
//右邊界.
if(base+BLOCK>=n){
s[n-base+1]=0;
}
else{
s[BLOCK+1] = a[base+BLOCK];
}
}

//----------------------------------------
//同步化 (確保共享記憶體已寫入)
//----------------------------------------
__syncthreads();

//----------------------------------------
//輸出三點加權平均值
//----------------------------------------
if(base+t<n){
b[base+t]=(s[t]+2*s[t+1]+s[t+2])*0.25;
}

};


//--------------------------------------------------------
//主函式.
//--------------------------------------------------------
int main(){
//--------------------------------------------------
//參數.
//--------------------------------------------------
int num=10*1000*1000;
int loop=130; //測試迴圈次數 (量時間用)

//--------------------------------------------------
//配置主機記憶體.
//--------------------------------------------------
float* a=new float[num];
float* b=new float[num];
float* bg=new float[num];
float* bs=new float[num];

//--------------------------------------------------
//配置裝置記憶體.
//--------------------------------------------------
float *GA, *GB;
cudaMalloc((void**) &GA, sizeof(float)*num);
cudaMalloc((void**) &GB, sizeof(float)*num);

//--------------------------------------------------
//初始化(亂數) & 複製資料到顯示卡的 DRAM.
//--------------------------------------------------
for(int k=0; k<num; k++){
a[k]=(float)rand()/RAND_MAX;
b[k]=bg[k]=bs[k]=0;
}
cudaMemcpy(GA, a, sizeof(float)*num, cudaMemcpyHostToDevice);

//--------------------------------------------------
//Test(1): smooth_host
//--------------------------------------------------
double t_host=(double)clock()/CLOCKS_PER_SEC;
for(int k=0; k<loop; k++){
smooth_host(b,a,num);
}
t_host=((double)clock()/CLOCKS_PER_SEC-t_host)/loop;


//--------------------------------------------------
//Test(2): smooth_global (GRID*BLOCK 必需大於 num).
//--------------------------------------------------
double t_global=(double)clock()/CLOCKS_PER_SEC;
cudaThreadSynchronize();
for(int k=0; k<loop; k++){
smooth_global<<<num/BLOCK+1,BLOCK>>>(GB,GA,num);
}
cudaThreadSynchronize();
t_global=((double)clock()/CLOCKS_PER_SEC-t_global)/loop;

//下載裝置記憶體內容到主機上.
cudaMemcpy(bg, GB, sizeof(float)*num, cudaMemcpyDeviceToHost);

//--------------------------------------------------
//Test(3): smooth_shared (GRID*BLOCK 必需大於 num).
//--------------------------------------------------
double t_shared=(double)clock()/CLOCKS_PER_SEC;
cudaThreadSynchronize();
for(int k=0; k<loop; k++){
smooth_shared<<<num/BLOCK+1,BLOCK>>>(GB,GA,num);
}
cudaThreadSynchronize();
t_shared=((double)clock()/CLOCKS_PER_SEC-t_shared)/loop;

//下載裝置記憶體內容到主機上.
cudaMemcpy(bs, GB, sizeof(float)*num, cudaMemcpyDeviceToHost);

//--------------------------------------------------
//比較正確性
//--------------------------------------------------
double sum_dg2=0, sum_ds2=0, sum_b2=0;
for(int k=0; k<num; k++){
double dg=bg[k]-b[k];
double ds=bs[k]-b[k];

sum_b2+=b[k]*b[k];
sum_dg2+=dg*dg;
sum_ds2+=ds*ds;
}

//--------------------------------------------------
//報告
//--------------------------------------------------
//組態.
printf("vector size: %d \n", num);
printf("\n");

//時間.
printf("Smooth_Host: %g ms\n", t_host*1000);
printf("Smooth_Global: %g ms\n", t_global*1000);
printf("Smooth_Shared: %g ms\n", t_shared*1000);
printf("\n");

//相對誤差.
printf("Diff(Smooth_Global): %g \n", sqrt(sum_dg2/sum_b2));
printf("Diff(Smooth_Shared): %g \n", sqrt(sum_ds2/sum_b2));
printf("\n");

//--------------------------------------------------
//釋放裝置記憶體.
//--------------------------------------------------
cudaFree(GA);
cudaFree(GB);
delete [] a;
delete [] b;
delete [] bg;
delete [] bs;
return 0;
}


//--------------------------------------------------------
//範例(5): 執行結果 (測試 10M 個 float)
//--------------------------------------------------------


vector size: 10000000

Smooth_Host: 36.9231 ms
Smooth_Global: 14.4615 ms
Smooth_Shared: 5.07692 ms

Diff(Smooth_Global): 3.83862e-08
Diff(Smooth_Shared): 3.83862e-08


(1) 這次測試的機器比較爛: P4-3.2 (prescott) vs. 9600GT
不過我們仍可看到共享記憶體使載入的資料量變為 1/3 所得到的增速

(2) 在 smooth_shared() 裡我們用 2 個 warp 使得條件判斷可以獨立,
如果第 2 個 if(t==32) 改成 if(t==16) 或其它小於 32 的值,
也就是和第 1 個 if 使用同一個 warp, 則速度會變慢, 有興趣的朋友
可以去試試看，warp 不打算在新手篇講，之後硬體時才詳細討論。

(3) 測試效能時使用 cudaThreadSynchronize() 同步主機和裝置核心，
以免量到錯誤的時間

(4) 在 smooth_shared() 裡使用 __syncthreads() 同步化執行緒，
以免在計算 output 時仍有共享記憶體還沒完成寫入動作，
卻有執行緒已經需要使用它的資料。


============================================================================
第六招 執行緒同步 (網格、區塊)
============================================================================
同步執行緒有兩個函式，分別是 __syncthreads() 和 cudaThreadSynchronize()

-----------------------------------------------------------------
同步化函式 使用地點 功能
-----------------------------------------------------------------
__syncthreads() 核心程序中 同步化【區塊內的執行緒】
cudaThreadSynchronize() 主機程序中 同步化【核心和主機程序】
-----------------------------------------------------------------

(1) 在 kernel 中，使用 __syncthreads() 來進行區塊內的執行緒的同步，
避免資料時序上的問題 (來自不同 threads)，時常和共享記憶體一起使用，
在範例(5)中示範了使用 __syncthreads() 來隔開共享記憶體的【寫入週期】
和【讀取週期】，避免 WAR 錯誤 (write after read)。

(2) 在主機程序中，使用 cudaThreadSynchronize() 來進行核心和主機程序的同步，
範例(5)中示範了用它來避免量到不正確的主機時間 (kernel仍未完成就量時間)，
因為主機的程序和裝置程序預設是不同步的 (直到下載結果資料之前)，這個 API
可以強迫它們同步。


--
。o O　○。o O　○。o O　○。o O　○。o O　○。o
國網 CUDA 中文教學 DVD 影片 (免費線上版)
請至國網的教育訓練網登入 https://edu.nchc.org.tw

BT 牌的種子下載點
http://www.badongo.com/file/12156676

--