CUDA 程式設計(7) -- 來玩泡泡吧 - 顯卡

第七章 來玩泡泡吧

今天我們來玩泡泡排序法吧!

◆ 演算法

為了能夠平行化, 使用交錯的方式對配對的資料進行排序, 如下圖所示, 假設
有 10 筆資料, 我們先對位址 (0 1),(2 3),... 的 0 based 資料對進行排序,
然後再對位址 (1 2),(3 4),... 的 1 based 資料對進行排序, 如此一來,
每次進程可讓最小的往前浮上兩個單位, 總共要執行 N/2 次才會浮到最前面,
每次比對數量為 N-1, 所以是 O(N^2/2) 的泡泡法.

(0 1)(2 3)(4 5)(6 7)(8 9) 0-based
0 (1 2)(3 4)(5 6)(7 8) 9 1-based

若資料量為奇數時, 例如 N=11, 則執行的方式如下,

(0 1)(2 3)(4 5)(6 7)(8 9) 10 0-based
0 (1 2)(3 4)(5 6)(7 8)(9 10) 1-based

所以端點處和偶數時略有不同, 玩泡泡時要注意一下.


◆ 單一區塊範例

以下示範用共享記憶體來做泡泡排序, 每個執行緒負責比較一對相鄰的資料,
因為單一區塊最大的執行緒個數是 512, 所以這個版本最多的輸入資料數為

N = 512*2+1 = 1025

每次進程先比較 0-based 配對, 再比較 1-based 配對, 總共要執行 N/2 次,
也就是 blockDim 次, 程式碼放在附錄 bubble1.cu, 在 GTX 8800 Ultra 上
執行的結果如下 (host 為Intel Q6600)

time[gpu]: 0.72 ms
time[host]: 1.36 ms
ratio: 1.88889 x
------------------------
test[gpu]: pass
test[host]: pass

效能只有 1.8 倍, 因為它只用到一個區塊, 也就是只用到 G80 上面 16 個
多處理器 (MP, multiprocessors) 的其中之一, 如果能夠使用很多個區塊
不僅效能能提昇數倍, 而且也能加大排序資料量.


◆ 練習

(1) 將範例改成用全域記憶體的版本, 並和使用共享記憶體的範例比較效能.

(2) 將泡泡排序改為多區塊版本, 首先將資料分段, 然後使用多區塊對每段
進行排序, 其分段法仿照範例, 例如要排序 10*512 筆資料, 便可設定
10 個區塊, 每個區塊使用 256 個執行緒做區塊的泡泡排序, 整個網格
先做 0-based, 再做 256-based, 每次發出的區塊如下

0-based (0~511), (512~1023), (1024~1535), ...
256-based (256~767), (768~1279), (1280~1791), ...

因為每次最小的群組都會往前步進一個區塊, 所以對 10*512 而言, 只要
讓它 loop 個 10 次, 就可以把資料全部排序好了.
(Note: 網格發出可在 host 上控制, 結果於下禮拜公佈)

(3) 試估計多區塊版本的效能, 並與結果比較. (需考慮 MP 個數)



========================================================================
bubble1.cu
========================================================================
#include <cuda.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//----------------------------------------------
// 排序 N 個 float 元素 (N=1~1025)
// 使用到的記憶體大小為 SIZE 個 bytes
// 只使用單一區塊
// 區塊大小為 N/2
//----------------------------------------------
#define N 1024
#define SIZE (N*sizeof(float))

#define GRID 1
#define BLOCK (N/2)

#define testLoop 1000 //測試效能時的 loop 數

//----------------------------------------------
// 交換函式 (host 和 kernel 都可以使用)
// 因為加了 __host__ 和 __device__ 兩個標籤
//----------------------------------------------
inline __host__ __device__ void swap(float& a, float& b){
float c=a;
a=b;
b=c;
}

//----------------------------------------------
// 泡泡的 kernel (由小到大排列 N 個元素 a->r)
//----------------------------------------------
__global__ void bubble(float *r, float *a){
//*** blockDim=N/2 ***
int j=threadIdx.x; //j=0,1,2,...blockDim-1
int k=2*threadIdx.x; //k=0,2,4,...2*(blockDim-1) 配對的基底索引

//配置共享記憶體
__shared__ float s[N+20];

//載入資料到共享記憶體
__syncthreads(); //同步化執行緒, 加速載入速度 (合併讀取 coalesced)
s[j]=a[j]; //使用全部執行緒一起載入前半段 (0~N/2-1)
s[j+N/2]=a[j+N/2]; //使用全部執行緒一起載入後半段 (N/2~N-1)
if(j==0){
//若 N 為奇數時, 還要多載入一個尾巴, 只使用第 0 個執行緒
s[N-1]=a[N-1];
}

//開始泡泡排序
for(int loop=0; loop<=N/2; loop++){
//排列 0 based 配對資料 (0,1) (2,3) (4,5) ....
__syncthreads(); //同步化確保共享記憶體已寫入
if(s[k]>s[k+1]){
swap(s[k],s[k+1]);
}

//排列 1 based 配對資料 (1,2) (3,4) (5,6) ....
__syncthreads(); //同步化確保共享記憶體已寫入
if(s[k+1]>s[k+2]){
if(k<N-2) //若 N 為偶數時, 最後一個執行緒不作用
swap(s[k+1],s[k+2]);
}
}

//轉出資料到全域記憶體
__syncthreads();
r[j]=s[j];
r[j+N/2]=s[j+N/2];
if(j==0){
r[N-1]=s[N-1];
}

}


//----------------------------------------------
// 泡泡的 host 函數
//----------------------------------------------
void bubble_host(float *r, float *a){
//載入資料
for(int k=0; k<N; k++){
r[k]=a[k];
}

for(int loop=0; loop<=N/2; loop++){
//排列 0 based 配對資料
for(int k=0; k<N-1; k+=2){
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
//排列 1 based 配對資料
for(int k=1; k<N-1; k+=2){
if(r[k]>r[k+1]){
swap(r[k],r[k+1]);
}
}
}
}


//----------------------------------------------
// 主程式
//----------------------------------------------
int main(){
//配置 host 記憶體
float *a=(float*)malloc(SIZE);
float *b=(float*)malloc(SIZE);
float *c=(float*)malloc(SIZE);


//初始化
for(int k=0; k<N; k++){
a[k]=k;
c[k]=0;
}

//對陣列 a 洗牌
srand(time(0));
for(int k=0; k<2*N; k++){
int i=rand()%N;
int j=rand()%N;
swap(a[i],a[j]);
}

//配置 device 記憶體
float *ga, *gc;
cudaMalloc((void**)&ga, SIZE);
cudaMalloc((void**)&gc, SIZE);

//載入 (順便載入 c 來清空裝置記憶體內容)
cudaMemcpy(ga, a, SIZE, cudaMemcpyHostToDevice);
cudaMemcpy(gc, c, SIZE, cudaMemcpyHostToDevice);


//測試 kernel 效能
double t0=(double)clock()/CLOCKS_PER_SEC;
for(int k=0; k<testLoop; k++){
//呼叫 kernel (此為單一 block 的版本)
bubble<<<1,BLOCK>>>(gc,ga);

//同步化執行緒, 避免還沒做完, 量到不正確的時間
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;

//顯示計算時間, 並比較
printf("time[gpu]: %g ms\n",t0*1000);
printf("time[host]: %g ms\n",t1*1000);
printf("ratio: %g x\n",t1/t0);

//讀出 device 資料
cudaMemcpy(c, gc, SIZE, cudaMemcpyDeviceToHost);


//測試 device 結果的正確性
printf("------------------------\n");
bool flag=true;
for(int k=0; k<N; k++){
if(c[k]!=k){
flag=false;
break;
}
}
printf("test[gpu]: %s\n",flag?"pass":"fail");

//測試 host 結果的正確性
flag=true;
for(int k=0; k<N; k++){
if(b[k]!=k){
flag=false;
break;
}
}
printf("test[host]: %s\n",flag?"pass":"fail");

//釋放記憶體
cudaFree(ga);
cudaFree(gc);
free(a);
free(b);
free(c);


return 0;
}
========================================================================

--