CUDA 程式設計(9) -- 泡泡龍 part 2 - 顯卡

第 9 章 多區塊的泡泡排序

運動完回來了，繼續來玩上次的泡泡吧，這次要做多區塊的版本

=================================================================
◆ SIMT 對傳統演算法帶來的挑戰
=================================================================
上禮拜的演算法，大家可能會覺得不太習慣，所以開始之前，要解釋一下：

//傳統的泡泡排序 (單鍊連續) -----------------------
for(int m=N; m>0; m--){
for(int k=0; k<m-1; k++){
if(a[k]>a[k+1]){
swap(a[k],a[k+1]);
}
}
}

//交錯的泡泡排序 (雙鍊錯開) ------------------------
for(int m=N/2+1; m>0; m--){
__syncthreads();
for(int k=0; k<N-1; k+=2){ //0-based
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
__syncthreads();
for(int k=1; k<N-1; k+=2){ //1-based
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
}
//-------------------------------------------------

這樣的更動的理由，就是在上一篇 (8) SIMT vs OpenMP 中提及的：

(1)程式碼正交性
(2)序列化問題

SIMT 要避免執行緒搞在一起，不然競爭之下，就不知道寫入的是什麼哩，
例如排列 (5,4,3)，兩個 SP 同時作用，平行讀取到的是同一個時段的值，
傳統的泡泡法，(5,4) 要交換，(4,3)也要交換，所以就會變成 (4,?,4)，
中間的就看那個 SP 比較慢寫入，可能是 5 也可能是 3，不過不管是誰，
結果都是錯的，所以才需要錯開，讓執行緒處理的資料沒有重疊 (正交化)，
然後在中間塞入 __syncthreads()，好讓兩個 base 循序的部份得以正確。

很多傳統的好演算法都是基於循序過程，就拿排序來說好了，像 HeapSort、
QuickSort，循序的味道都很重，是否都能用 SIMT 來做，是值得研究的問題，
但軟體和硬體總是交互成長，我想 SIMT 應該會刺激演算法的革新，而演算法
發展所遇到的困難，也會帶動 SIMT 架構上的改進。


=================================================================
◆ 多區塊的泡泡排序
=================================================================
多區塊的泡泡，要分兩層來做

(1) BLOCK 層的泡泡 -- 就是上禮拜的單一區塊版本
(2) GRID 層的泡泡 -- 需要用 HOST 來派遣

例如下面要排序 16 個數字 9,5,8,4, 7,6,3,2, 6,8,7,4, 3,2,1,0 我們先把它分成
4 個區塊，每個區塊就用上次的單區塊版本來排序 (1a)，這時可以設定 gridDim=4，
讓這些區塊使用多個 MP 同時排序，完成後再把每個區塊的結果切成兩半，跟另一個
相鄰的半區塊結合，形成一個完整的區塊，設定 gridDim=3 對這三塊進行排序 (1b)，
我們留意到每經過一個 ab 流程，最小的半區塊 (尾巴的01) 就會往前步進一個區塊，
令 G=N/BLOCK，只需經過 G-1 個 ab 週期，便可將它推到最前面的區塊，然後再
過一個週期，就可完成區塊內排序，整個工作就完成了，需要跑 G 次。


---------------------
9584 7632 6874 3210 ------ (1a) 排列 4 個區塊
4589 2367 4678 0123

45 8923 6746 7801 23 ------ (1b) 排列 3 個區塊
45 2389 4667 0178 23
---------------------
4523 8946 6701 7823 ------ (2a)
2345 4689 0167 2378

23 4546 8901 6723 78 ------ (2b)
23 4456 0189 2367 78
---------------------
2344 5601 8923 6778 ------ (3a)
2344 0156 2389 6778

23 4401 5623 8967 78 ------ (3b)
23 0144 2356 6789 78
---------------------
2301 4423 5667 8978 ------ (4a)
0123 2344 5667 7889

01 2323 4456 6778 89 ------ (4b)
01 2233 4456 6778 89
---------------------

程式碼列於最後面的附錄中，值得注意的是不同區塊間無法在 device 層級同步化，
而這種交錯的方式在各進程之間有循序性，需要前面的資料正確寫入，所以才把它
拆成多個 host 控制的 kernel 呼叫，中間塞入 cudaThreadSynchronize() 這個 API
來同步。


=================================================================
◆ 演算法效能評估
=================================================================

(1) 單區塊排序 B 個元素，需經過 B/2 個進程，每個進程比較 B 個元素，
所以時間約為

T1 = O(B^2/2) = alpha * B^2/2 (令系數為 alpha)

(2) 排序 N=G*B 筆資料, 要排 G 個區塊, 每個區塊排序 B 個, 每次的
ab進程要發出 2G-1 個區塊, 迴圈 G 次, 若全部用單區塊版本來排序，
所需的時間約為

TG = (2G-1)*G*T1 = O((G*B)^2) = O(N^2)

看起來比 Host 要差, 因為由(1)可知，Host 所需的時間約為

TH = O(N^2/2) = beta * N^2/2 (令系數為 beta)

多區塊的版本引入 factor 2

(3) 考慮多處理器 MP 帶來的增速 (工作量分擔), 所需的時間變為

TM = TG/MP = O(N^2/MP) = alpha * N^2/MP

(4) 考慮單區塊和 HOST 的效能比 R1 = beta/alpha，估計多區塊和 HOST 的
效能比約為

RG = TH/TM = R1*MP/2

(5) 上次執行的結果 R1 約為 1.9 倍，所以多區塊版本和 HOST 的效能比約為

RG ~ MP

可以用來測試顯示卡有幾個多處理器 (MP, multiprocessors), 當然這個
1.9 倍是以 Intel Q6600 做準，其它 CPU 還是用 TM 來看會比較單純


Note: 泡泡排序的核心每載入 & 輸出 N 個元素，就要對它進行 N^2/2 個操作,
平均每傳送一個元素，所搭配的運算量為 N/4，當 N 是 512 時 N/4=128
所以它的瓶頸是在運算(或者說是 branch)，不是記憶體傳輸。


=================================================================
◆ 測試結果
=================================================================
以下是一些機器上的測試結果，大家可以看到跟我們上面的估計差不多，
和主機的效能比約等於 multiprocessor 的個數，而且因為它的瓶頸在計算，
所以透過 PCI-E 傳輸資料的時間跟本不重要 (詳見 code)

8800 GTX Ultra vs. Intel(R) Core(TM)2 Quad CPU Q6600 @ 2.40GHz
-------------------------------------
NUM=131072 GRID=128 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 1330 ms
time[host]: 22940 ms
ratio: 17.2481 x

-------------------------------------
NUM=262144 GRID=256 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 5330 ms
time[host]: 91850 ms
ratio: 17.2326 x


GTX 280 vs. Intel(R) Core(TM)2 Quad CPU Q9400 @ 2.66GHz
-------------------------------------
NUM = 262144 GRID=256 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 2920 ms
time[host]: 84070 ms
ratio: 28.7911 x

-------------------------------------
NUM = 524288 GRID=512 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 11470 ms
time[host]: 348880 ms
ratio: 30.4167 x


C1060 vs. Intel(R) Core(TM)2 Quad CPU Q9300 @ 2.50GHz
-------------------------------------
NUM=131072 GRID=128 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 830 ms
time[host]: 22450 ms
ratio: 27.0482 x


-------------------------------------
NUM=262144 GRID=256 BLOCK=512
-------------------------------------
test[gpu]: pass
test[host]: pass
time[gpu]: 2900 ms
time[host]: 90260 ms
ratio: 31.1241 x

Note: 網格數若設太小，效能會比較差，因為 MP 無法吃盡網格，會閒置，
而且這些額外的時間比較沒有機會被其它成功填滿的網格行程分擔掉，
還有無法整除的情況，如 GTX 280 系列的 MP 有 30 個，網格數若設 128，
則 12830=4…8 餘下的 8 個區塊在執行時其它的 22 個 MP 閒置著，
所以平均效能就會比較弱一點 (這邊測試的 8800 GTX 網格數都可以整除，
所以效能會比預計的稍微高一點，這是其中一個原因)


=================================================================
◆ 附錄: 多區塊的泡泡程式碼
=================================================================
#include <cuda.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//----------------------------------------------
// 請用 BLOCK 和 GRID 設定 NUM 的大小
//----------------------------------------------
#define BLOCK 512 //區塊大小
#define GRID 32 //網格大小
#define NUM (2*BLOCK*GRID) //需要排列的元素個數

#define SIZE (NUM*sizeof(float))
#define testLoop 1

//----------------------------------------------
inline __host__ __device__ void swap(float& a, float& b){
float c=a;
a=b;
b=c;
}

//----------------------------------------------
// 單一區塊的泡泡排序 (偶數版, 只排 BLOCK*2 個元素, 要對 blockIdx 定位)
//----------------------------------------------
__global__ void bubble(float *r, float *a){
const int j=threadIdx.x;
const int k=2*threadIdx.x;
a=a+blockIdx.x*BLOCK*2;
r=r+blockIdx.x*BLOCK*2;
__shared__ float s[BLOCK*2+10];

__syncthreads();
s[j]=a[j];
s[j+BLOCK]=a[j+BLOCK];

for(int loop=0; loop<=BLOCK; loop++){
__syncthreads();
if(s[k]>s[k+1]){
swap(s[k],s[k+1]);
}
__syncthreads();
if(j<BLOCK-1)
if(s[k+1]>s[k+2]){
swap(s[k+1],s[k+2]);
}
}

__syncthreads();
r[j]=s[j];
r[j+BLOCK]=s[j+BLOCK];

}


//----------------------------------------------
// Host 的泡泡排序
//----------------------------------------------
void bubble_host(float *r, float *a){
for(int k=0; k<NUM; k++){
r[k]=a[k];
}

for(int loop=0; loop<=NUM/2; loop++){
for(int k=0; k<NUM-1; k+=2){
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
for(int k=1; k<NUM-1; k+=2){
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
}
}

//----------------------------------------------
int main(){

srand(time(0));
float *a=(float*)malloc(SIZE);
float *b=(float*)malloc(SIZE);
float *c=(float*)malloc(SIZE);

printf("------------------------\n");
printf("NUM = %d\n",NUM);
printf("------------------------\n");


for(int k=0; k<NUM; k++){
a[k]=k;
c[k]=0;
}

for(int k=0; k<2*NUM; k++){
int i=rand()%NUM;
int j=rand()%NUM;
swap(a[i],a[j]);
}

float *gc;
cudaMalloc((void**)&gc, SIZE);
cudaMemcpy(gc, a, SIZE, cudaMemcpyHostToDevice);


//測試 GPU 多區塊泡泡效能
double t0=(double)clock()/CLOCKS_PER_SEC;
for(int k=0; k<testLoop; k++){
//連 PCI-E 載入的時間算進去都划算
cudaMemcpy(gc, a, SIZE, cudaMemcpyHostToDevice);
cudaThreadSynchronize();

for(int g=0; g<GRID; g++){

//進程 a
bubble<<<GRID,BLOCK>>>(gc,gc);
cudaThreadSynchronize();

//進程 b: 區塊數少1, 工作點平移半個區塊
bubble<<<GRID-1,BLOCK>>>(gc+BLOCK,gc+BLOCK);
cudaThreadSynchronize();

}
cudaMemcpy(c, gc, SIZE, cudaMemcpyDeviceToHost);
cudaThreadSynchronize();
}
t0=((double)clock()/CLOCKS_PER_SEC-t0)/testLoop;

//測試 Host 泡泡效能
double t1=(double)clock()/CLOCKS_PER_SEC;
for(int k=0; k<testLoop; k++){
bubble_host(b,a);
}
t1=((double)clock()/CLOCKS_PER_SEC-t1)/testLoop;




//測量準確性
bool flag=true;
for(int k=0; k<NUM; k++){
if(c[k]!=k){
flag=false;
break;
}
}
printf("test[gpu]: %s\n",flag?"pass":"fail");

flag=true;
for(int k=0; k<NUM; k++){
if(b[k]!=k){
flag=false;
break;
}
}
printf("test[host]: %s\n",flag?"pass":"fail");


printf("time[gpu]: %g ms\n",t0*1000);
printf("time[host]: %g ms\n",t1*1000);
printf("ratio: %g x\n",t1/t0);



printf("\n");

cudaFree(gc);
free(a);
free(b);
free(c);


return 0;
}







--