CUDA 程式設計(12) -- 速成篇(下) - 顯卡

============================================================================
第七招 合併存取
============================================================================

【合併存取】(coalesced I/O) 是 CUDA 中最基本且最重要的最佳化手段, 因為
GPU 的計算能力太強, 使得效能瓶頸卡在顯示記憶體到 GPU 之間的 I/O 上,
合併存取可讓多個顯示記憶體的交易合併成一次, 而加速記憶體的存取.

現階段 GPU 在合併存取上是自動發生的, 以半個 warp 為單位(16 個相鄰的執行緒),
如果它們的資料位址是連續的, 就會被合併, 所以使用上很簡單, 只要 threadIdx
對齊即可, 它可以合併 4, 8, 16 bytes 的資料, 成為一次 or 兩次的交易,
合併成的最大封包長度可為 32, 64, 128 bytes (其中 32 bytes 的封包只有在
版本 1.2 以後支援, 避免位址分散情況下的 overhead)，以下為連續的資料位址
的合併情況

16 x 4 bytes -> 64 bytes
16 x 8 bytes -> 128 bytes
16 x 16 bytes -> 2 x 128 bytes


【範例: 矩陣的 transpose】
========================================================================
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <cuda.h>

//transpose m*n matrix a to n*m matrix b on host ---------------------
void transpose_host(float* b, float* a, int m, int n){
for(int y=0; y<m; y++)
for(int x=0; x<n; x++){
b[x*m+y]=a[y*n+x];
}
}

//transpose naive (讀取合併). ----------------------------------------
__global__ void transpose_naive_cr(float* b, float* a, int m, int n){
int x=blockIdx.x*blockDim.x + threadIdx.x;
int y=blockIdx.y*blockDim.y + threadIdx.y;
if(y<m && x<n){
b[x*m+y]=a[y*n+x];
}
}

//transpose naive (寫入合併). ----------------------------------------
__global__ void transpose_naive_cw(float* b, float* a, int m, int n){
//x,y 的 threadIdx 足標對調.
int x=blockIdx.x*blockDim.x + threadIdx.y;
int y=blockIdx.y*blockDim.y + threadIdx.x;
if(y<m && x<n){
b[x*m+y]=a[y*n+x];
}
}

//transpose shared (讀取 & 寫入合併). -------------------------------
__global__ void transpose_shared(float* b, float* a, int m, int n){
//宣告共享記憶體.
__shared__ float s[256];

//讀取合併.
int x=blockIdx.x*blockDim.x + threadIdx.x;
int y=blockIdx.y*blockDim.y + threadIdx.y;
if(y<m && x<n){
int t=threadIdx.y*blockDim.x + threadIdx.x;
int i=y*n+x;
s[t]=a[i];
}
__syncthreads();

//寫入合併 (x,y 的 threadIdx 足標對調).
x=blockIdx.x*blockDim.x + threadIdx.y;
y=blockIdx.y*blockDim.y + threadIdx.x;
if(y<m && x<n){
//共享記憶體中的 threadIdx 足標亦要對調.
int t=threadIdx.x*blockDim.y + threadIdx.y;
int o=x*m+y;
b[o]=s[t];
}
}

//計算相對誤差. ----------------------------------------------------
double rd(float*a, float*b, int size){
double s=0, d=0;
for(int k=0; k<size; k++){
double w=a[k]-b[k];
s+=a[k]*a[k];
d+=w*w;
}
return sqrt(d/s);
}

//timer functions --------------------------------------------------
time_t timer[10];
void set_timer(int k=0){
timer[k]=clock();
}
double get_timer(int k=0){
return (double)(clock()-timer[k])/CLOCKS_PER_SEC;
}

//測試程式. (輸入 m*n 矩陣) ----------------------------------------
void test(int m, int n, int loop=100, int loop_host=10){
int size=m*n;
printf("matrix size: %d x %d\n", m,n);


//配置主記憶體, 並設定初始值.
float *a=new float[size];
float *b=new float[size];
float *c=new float[size];

for(int k=0; k<size; k++){
a[k]=(float)rand()*2/RAND_MAX-1;
b[k]=0;
}


//配置顯示記憶體, 載入資料.
float *ga, *gb;
cudaMalloc((void**)&ga, size*sizeof(float));
cudaMalloc((void**)&gb, size*sizeof(float));
cudaMemcpy(ga, a, size*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(gb, b, size*sizeof(float), cudaMemcpyHostToDevice);


//網格區塊設定.
dim3 grid(n/16+1,m/16+1,1); //網格要含蓋所有範圍, 所以除完要加 1.
dim3 block(16,16,1); //區塊設定 16x16.
printf("grid(%d,%d,%d)\n",grid.x,grid.y,grid.z);
printf("block(%d,%d,%d)\n",block.x,block.y,block.z);



//測試 transpose_host 函數.
set_timer();
for(int k=0; k<loop_host; k++){
transpose_host(b,a,m,n);
}
double t0=get_timer()/loop_host;
printf("host time: %g ms\n",t0*1000);



//測試 transpose_naive_cr 函數.
cudaMemset(gb,0,size*sizeof(float));
transpose_naive_cr<<<grid,block>>>(gb,ga,m,n);
cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost);

set_timer();
cudaThreadSynchronize();
for(int k=0; k<loop; k++){
transpose_naive_cr<<<grid,block>>>(gb,ga,m,n);
}
cudaThreadSynchronize();
double t1=get_timer()/loop;
printf("naive(r) time: %g ms (%dx)\t error: %g\n",t1*1000,
(int)(t0/t1),rd(b,c,size));

//測試 transpose_naive_cw 函數.
cudaMemset(gb,0,size*sizeof(float));
transpose_naive_cw<<<grid,block>>>(gb,ga,m,n);
cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost);

set_timer();
cudaThreadSynchronize();
for(int k=0; k<loop; k++){
transpose_naive_cw<<<grid,block>>>(gb,ga,m,n);
}
cudaThreadSynchronize();
double t2=get_timer()/loop;

printf("naive(w) time: %g ms (%dx)\t error: %g\n",t2*1000,
(int)(t0/t2),rd(b,c,size));

//測試 transpose_shared 函數.
cudaMemset(gb,0,size*sizeof(float));
transpose_shared<<<grid,block>>>(gb,ga,m,n);
cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost);

set_timer();
cudaThreadSynchronize();
for(int k=0; k<loop; k++){
transpose_shared<<<grid,block>>>(gb,ga,m,n);
}
cudaThreadSynchronize();
double t3=get_timer()/loop;
printf("shared(r/w)time: %g ms (%dx)\t error: %g\n",t3*1000,
(int)(t0/t3),rd(b,c,size));


//釋放記憶體
cudaFree(ga);
cudaFree(gb);

delete [] a;
delete [] b;
delete [] c;
}

//主函數 ------------------------------------------------------------------
int main(){
srand(time(0));

printf("----------------------------\n");
test(2047,4000);
printf("----------------------------\n");
test(2048,4000);
printf("----------------------------\n");
test(2049,4000);
return 0;
}


【說明】
==========================================================================
(1) 程式中使用 dim3(16,16,1) 的區塊設定, 所以 threadIdx.x=0~15 形成所謂的
「半個 warp」.

(2) 在 transpose_naive_cr() 中 int x=blockIdx.x*blockDim.x + threadIdx.x;
包含半個 warp (threadIdx.x), 所以 x 在半個 warp 中連續, 讀取 a[y*n+x]
被合併, 形成所謂的【讀取合併】(coalesced read)

(3) 在 transpose_naive_cw() 中, 把 naive_cr 版本 (x,y) 的 threadIdx 足標對調,
換成 y 在半個 warp 上是連續的, 使得寫入 b[x*m+y] 被合併

int x=blockIdx.x*blockDim.x + threadIdx.y;
int y=blockIdx.y*blockDim.y + threadIdx.x;

形成所謂的【寫入合併】(coalesced write).

(4) 在 transpose_shared() 中, 上半週期應用 half warp 到讀取位址的連續

int x=blockIdx.x*blockDim.x + threadIdx.x;
int y=blockIdx.y*blockDim.y + threadIdx.y;
int i=y*n+x; <--- (OOO) + threadIdx.x

形成【讀取合併】, 並把資料存於快速的共享記憶體中.

下半週期要對調 threadIdx 足標, 使得 half warp 在寫入位址連續

int x=blockIdx.x*blockDim.x + threadIdx.y;
int y=blockIdx.y*blockDim.y + threadIdx.x;
int o=x*m+y; <--- (XXX) + threadIdx.x

形成【寫入合併】

(5) transpose_shared() 中使用共享記憶體的理由是資料的讀取和寫入使用不同的
執行緒 (threadIdx 足標對調), 所以執行緒之間必需交換資料, 注意足標對調時
共享記憶體中的 threadIdx 足標亦要對調, 以對應到相對的元素

讀取 int t=threadIdx.y*blockDim.x + threadIdx.x;
寫入 int t=threadIdx.x*blockDim.y + threadIdx.y;

【執行結果】
=========================================================================
這次測試的機型為 GTX260 vs Intel E8400(3.00GHz),
我們可以看出, 合併存取的另一個好處是減少 bank conflict 對它的影響,
當 host 因為 bank conflict 變慢時, GTX260 還是維持差不多的水準.
平均 shared 版的 transpose 的效能大概是 host 的 20 多倍

----------------------------
matrix size: 2047 x 4000
grid(251,128,1)
block(16,16,1)
host time: 47 ms
naive(r) time: 8 ms (5x) error: 0 //讀取合併 (naive_cr)
naive(w) time: 4 ms (11x) error: 0 //寫入合併 (naive_cw)
shared(r/w)time: 2.1 ms (22x) error: 0 //讀寫合併 (shared)
----------------------------
matrix size: 2048 x 4000 (這是可以產生 bank conflict 的情況)
grid(251,129,1)
block(16,16,1)
host time: 305 ms
naive(r) time: 8.5 ms (35x) error: 0
naive(w) time: 4.1 ms (74x) error: 0
shared(r/w)time: 2.1 ms (145x) error: 0 //因為 host 太遜了
----------------------------
matrix size: 2049 x 4000
grid(251,129,1)
block(16,16,1)
host time: 46 ms
naive(r) time: 8.1 ms (5x) error: 0
naive(w) time: 4.2 ms (10x) error: 0
shared(r/w)time: 2.2 ms (20x) error: 0


【GT200 的改進】
==========================================================================
在 G80/G90 系列 (compute 1.0 & 1.1), 合併存取必需要符合許多要求, 例如
半個 warp 中的資料位址要按照 threadIdx 的順序, 封包的邊界位址必需滿足
64 bytes 或 128 bytes 的對齊要求, 否則會放棄合併, 產生 16 個記憶體要求.

在 GT200 系列 (compute 1.2 之後), 合併存取有大幅度改進, 容許半個 warp 的
資料位址不按照順序排列, 不需 16 個執行緒全部連續, 容許部份連續的情況
(此時會進行部份位址合併, 並拆成數個封包發出), 在節省頻寬方面也會選擇
最有效率的方式分解封包, 並引進 32 bytes 的小型封包, 避免較小的連續區塊
仍需使用大的合併封包, 整體記憶體效能比 G80/G90 系列好甚多.


【範例：GTX260 半個 warp 的存取順序】
==========================================================================
(1) 將範例中的 transpose_naive_cr 的 x 反序
__global__ void transpose_naive_cr(float* b, float* a, int m, int n){
int x=blockIdx.x*blockDim.x + (15-threadIdx.x);
int y=blockIdx.y*blockDim.y + threadIdx.y;
if(y<m && x<n){
b[x*m+y]=a[y*n+x];
}
}

(2) 使用人工亂數, 將範例中的 transpose_naive_cw 改為隨機順序
__global__ void transpose_naive_cw(float* b, float* a, int m, int n){
int x=blockIdx.x*blockDim.x + threadIdx.y;
int y=blockIdx.y*blockDim.y + (threadIdx.x*7+5)%16;
if(y<m && x<n){
b[x*m+y]=a[y*n+x];
}
}


半個 warp 的存取順序對應如下
-----------------------
tid.x (tid.x*7+5)%16
-----------------------
0 5
1 12
2 3
3 10
4 1
5 8
6 15
7 6
8 13
9 4
10 11
11 2
12 9
13 0
14 7
15 14
-----------------------


【執行結果】
============================================================================
在 GTX260 上執行的時間仍相同, 基本上不太受存取順序的影響。

----------------------------
matrix size: 2047 x 4000
grid(251,128,1)
block(16,16,1)
host time: 47 ms
naive(r) time: 7.8 ms (6x) error: 0 //naive_cr 版
naive(w) time: 4.1 ms (11x) error: 0 //naive_cw 版
shared(r/w)time: 2 ms (23x) error: 0 //shared 版
----------------------------
matrix size: 2048 x 4000 (這是可以產生 bank conflict 的情況)
grid(251,129,1)
block(16,16,1)
host time: 308 ms
naive(r) time: 8.5 ms (36x) error: 0
naive(w) time: 4.1 ms (75x) error: 0
shared(r/w)time: 2 ms (153x) error: 0
----------------------------
matrix size: 2049 x 4000
grid(251,129,1)
block(16,16,1)
host time: 46 ms
naive(r) time: 8.1 ms (5x) error: 0
naive(w) time: 4.1 ms (11x) error: 0
shared(r/w)time: 2 ms (22x) error: 0
----------------------------


【補充: warp】
==========================================================================
warp 是 GPU 硬體上一次執行的實際單位, 一個區塊可分成數個 warp,
分時在 multiprocessor 中執行, 例如

128 threads 的區塊 => 4 warps
200 threads 的區塊 => 6+1 warps (多出的未填滿仍算 1 個 warp)

詳細情形在之後硬體篇會再 review, warp 和 threadIdx 的關係如下

-------------------
warp threadIdx
-------------------
0 0~31
1 32~63
2 64~95
3 96~127
...
-------------------

半個 warp 是記憶體合併的單位, 也就是一個 warp 的記憶體讀寫可分成兩組
獨立進行合併, 和 threadIdx 的關係如下

-----------------------
half warp threadIdx
-----------------------
0 0~15
1 16~31
2 32~47
3 48~63
...
-----------------------


【補充: threadIdx 順序】
=========================================================================
若區塊中的執行緒使用 3D 結構安排時, 其順序是 threadIdx.x 在最裡面,
然後是 threadIdx.y, 最外層是 threadIdx.z, 結構類似 C/C++ 的三維陣列
在記憶體中的順序

threads[z][y][x];

打平成一維

tid = threadIdx.z * (blockDim.y*blockDim.x)
threadIdx.y * (blockDim.x)
threadIdx.x;

所以在範例中使用 dim(16,16,1) 區塊, 其 threadIdx.x 維度包含 16 個執行緒,
剛好每一個 x 維度形成半個 warp, 整個區塊有 8 個 warps.



--
。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

--