改良第一個程式

在上一篇文章中，我們做了一個計算一大堆數字的平方和的程式。不過，我們也提到這個程式的執行效率並不理想。當然，實際上來說，如果只是要做計算平方和的動作，用 CPU 做會比用 GPU 快得多。這是因為平方和的計算並不需要太多運算能力，所以幾乎都是被記憶體頻寬所限制。因此，光是把資料複製到顯示記憶體上的這個動作，所需要的時間，可能已經和直接在 CPU 上進行計算差不多了。

不過，如果進行平方和的計算，只是一個更複雜的計算過程的一部份的話，那麼當然在 GPU 上計算還是有它的好處的。而且，如果資料已經在顯示記憶體上（例如在 GPU 上透過某種演演算法產生），那麼，使用 GPU 進行這樣的運算，還是會比較快的。

剛才也提到了，由於這個計算的主要瓶頸是記憶體頻寬，所以，理論上顯示卡的記憶體頻寬是相當大的。這裡我們就來看看，倒底我們的第一個程式，能利用到多少記憶體頻寬。

程式的平行化
我們的第一個程式，並沒有利用到任何平行化的功能。整個程式只有一個 thread。在 GeForce 8800GT 上面，在 GPU 上執行的部份（稱為 "kernel"）大約花費 640M 個時脈。GeForce 8800GT 的執行單元的時脈是 1.5GHz，因此這表示它花費了約 0.43 秒的時間。1M 個 32 bits 數字的資料量是 4MB，因此，這個程式實際上使用的記憶體頻寬，只有 9.3MB/s 左右！這是非常糟糕的表現。

為什麼會有這樣差的表現呢？這是因為 GPU 的架構特性所造成的。在 CUDA 中，一般的資料複製到的顯示記憶體的部份，稱為 global memory。這些記憶體是沒有 cache 的，而且，存取 global memory 所需要的時間（即 latency）是非常長的，通常是數百個 cycles。由於我們的程式只有一個 thread，所以每次它讀取 global memory 的內容，就要等到實際讀取到資料、累加到 sum 之後，才能進行下一步。這就是為什麼它的表現會這麼的差。

由於 global memory 並沒有 cache，所以要避開巨大的 latency 的方法，就是要利用大量的 threads。假設現在有大量的 threads 在同時執行，那麼當一個 thread 讀取記憶體，開始等待結果的時候，GPU 就可以立刻切換到下一個 thread，並讀取下一個記憶體位置。因此，理想上當 thread 的數目夠多的時候，就可以完全把 global memory 的巨大 latency 隱藏起來了。

要怎麼把計算平方和的程式平行化呢？最簡單的方法，似乎就是把數字分成若干組，把各組數字分別計算平方和後，最後再把每組的和加總起來就可以了。一開始，我們可以把最後加總的動作，由 CPU 來進行。

首先，在 first_cuda.cu 中，在 #define DATA_SIZE 的後面增加一個 #define，設定 thread 的數目：

#define DATA_SIZE    1048576
#define THREAD_NUM   256

接著，把 kernel 程式改成：

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    const int tid = threadIdx.x;
    const int size = DATA_SIZE / THREAD_NUM;
    int sum = 0;
    int i;
    clock_t start;
    if(tid == 0) start = clock();
    for(i = tid * size; i < (tid + 1) * size; i++) {
       sum += num[i] * num[i];
    }
    result[tid] = sum;
    if(tid == 0) *time = clock() - start;
}

程式裡的 threadIdx 是 CUDA 的一個內建的變數，表示目前的 thread 是第幾個 thread（由 0 開始計算）。以我們的例子來說，會有 256 個 threads，所以同時會有 256 個 sumOfSquares 函式在執行，但每一個的 threadIdx.x 則分別會是 0 ~ 255。利用這個變數，我們就可以讓每一份函式執行時，對整個資料不同的部份計算平方和。另外，我們也讓計算時間的動作，只在 thread 0（即 threadIdx.x = 0 的時候）進行。

同樣的，由於會有 256 個計算結果，所以原來存放 result 的記憶體位置也要擴大。把 main 函式中的中間部份改成：

    int* gpudata, *result;
    clock_t* time;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
    cudaMalloc((void**) &time, sizeof(clock_t));
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);
    sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
    int sum[THREAD_NUM];
    clock_t time_used;
    cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM,
        cudaMemcpyDeviceToHost);
    cudaMemcpy(&time_used, time, sizeof(clock_t),
        cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);
    cudaFree(time);

可以注意到我們在呼叫 sumOfSquares 函式時，指定 THREAD_NUM 為 thread 的數目。最後，在 CPU 端把計算好的各組資料的平方和進行加總：

    int final_sum = 0;
    for(int i = 0; i < THREAD_NUM; i++) {
        final_sum += sum[i];
    }
    printf("sum: %d  time: %d\n", final_sum, time_used);
    final_sum = 0;
    for(int i = 0; i < DATA_SIZE; i++) {
        sum += data[i] * data[i];
    }
    printf("sum (CPU): %d\n", final_sum);

編譯後執行，確認結果和原來相同。

這個版本的程式，在 GeForce 8800GT 上執行，只需要約 8.3M cycles，比前一版程式快了 77 倍！這就是透過大量 thread 來隱藏 latency 所帶來的效果。

不過，如果計算一下它使用的記憶體頻寬，就會發現其實仍不是很理想，大約只有 723MB/s 而已。這和 GeForce 8800GT 所具有的記憶體頻寬是很大的差距。為什麼會這樣呢？

記憶體的存取模式
顯示卡上的記憶體是 DRAM，因此最有效率的存取方式，是以連續的方式存取。前面的程式，雖然看起來是連續存取記憶體位置（每個 thread 對一塊連續的數字計算平方和），但是我們要考慮到實際上 thread 的執行方式。前面提過，當一個 thread 在等待記憶體的資料時，GPU 會切換到下一個 thread。也就是說，實際上執行的順序是類似

    thread 0 -> thread 1 -> thread 2 -> ...

因此，在同一個 thread 中連續存取記憶體，在實際執行時反而不是連續了。要讓實際執行結果是連續的存取，我們應該要讓 thread 0 讀取第一個數字，thread 1 讀取第二個數字…依此類推。所以，我們可以把 kernel 程式改成如下：

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    const int tid = threadIdx.x;
    int sum = 0;
    int i;
    clock_t start;
    if(tid == 0) start = clock();
    for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
       sum += num[i] * num[i];
    }
    result[tid] = sum;
    if(tid == 0) *time = clock() - start;
}

編譯後執行，確認結果相同。

僅僅是這樣簡單的修改，實際執行的效率就有很大的差別。在 GeForce 8800GT 上，上面的程式執行需要的時脈是 2.6M cycles，又比前一版程式快了三倍。不過，這樣仍只有 2.3GB/s 的頻寬而已。

這是因為我們使用的 thread 數目還是不夠多的原因。理論上 256 個 threads 最多隻能隱藏 256 cycles 的 latency。但是 GPU 存取 global memory 時的 latency 可能高達 500 cycles 以上。如果增加 thread 數目，就可以看到更好的效率。例如，可以把 THREAD_NUM 改成 512。在 GeForce 8800GT 上，這可以讓執行花費的時間減少到 1.95M cycles。有些改進，但是仍不夠大。不幸的是，目前 GeForce 8800GT 一個 block 最多隻能有 512 個 threads，所以不能再增加了，而且，如果 thread 數目增加太多，那麼在 CPU 端要做的最後加總工作也會變多。

更多的平行化
前面提到了 block。在之前介紹呼叫 CUDA 函式時，也有提到 "block 數目" 這個引數。到目前為止，我們都只使用一個 block。究竟 block 是什麼呢？

在 CUDA 中，thread 是可以分組的，也就是 block。一個 block 中的 thread，具有一個共用的 shared memory，也可以進行同步工作。不同 block 之間的 thread 則不行。在我們的程式中，其實不太需要進行 thread 的同步動作，因此我們可以使用多個 block 來進一步增加 thread 的數目。

首先，在 #define DATA_SIZE 的地方，改成如下：

#define DATA_SIZE   1048576
#define BLOCK_NUM   32
#define THREAD_NUM   256

這表示我們會建立 32 個 blocks，每個 blocks 有 256 個 threads，總共有 32*256 = 8192 個 threads。

接著，我們把 kernel 部份改成：

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    int sum = 0;
    int i;
    if(tid == 0) time[bid] = clock();
    for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
        i += BLOCK_NUM * THREAD_NUM) {
       sum += num[i] * num[i];
    }
    result[bid * THREAD_NUM + tid] = sum;
    if(tid == 0) time[bid + BLOCK_NUM] = clock();
}

blockIdx.x 和 threadIdx.x 一樣是 CUDA 內建的變數，它表示的是目前的 block 編號。另外，注意到我們把計算時間的方式改成每個 block 都會記錄開始時間及結束時間。

main 函式部份，修改成：

    int* gpudata, *result;
    clock_t* time;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result,
        sizeof(int) * THREAD_NUM * BLOCK_NUM);
    cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);
    sumOfSquares<<>>(gpudata, result,
        time);
    int sum[THREAD_NUM * BLOCK_NUM];
    clock_t time_used[BLOCK_NUM * 2];
    cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM * BLOCK_NUM,
        cudaMemcpyDeviceToHost);
    cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
        cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);
    cudaFree(time);
    int final_sum = 0;
    for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) {
        final_sum += sum[i];
    }
    clock_t min_start, max_end;
    min_start = time_used[0];
    max_end = time_used[BLOCK_NUM];
    for(int i = 1; i < BLOCK_NUM; i++) {
        if(min_start > time_used[i])
            min_start = time_used[i];
        if(max_end < time_used[i + BLOCK_NUM])
            max_end = time_used[i + BLOCK_NUM];
    }
    printf("sum: %d  time: %d\n", final_sum, max_end - min_start);

基本上我們只是把 result 的大小變大，並修改計算時間的方式，把每個 block 最早的開始時間，和最晚的結束時間相減，取得總執行時間。

這個版本的程式，執行的時間減少很多，在 GeForce 8800GT 上只需要約 150K cycles，相當於 40GB/s 左右的頻寬。不過，它在 CPU 上執行的部份，需要的時間加長了（因為 CPU 現在需要加總 8192 個數字）。為了避免這個問題，我們可以讓每個 block 把自己的每個 thread 的計算結果進行加總。

Thread 的同步
前面提過，一個 block 內的 thread 可以有共用的記憶體，也可以進行同步。我們可以利用這一點，讓每個 block 內的所有 thread 把自己計算的結果加總起來。把 kernel 改成如下：

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    extern __shared__ int shared[];
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    int i;
    if(tid == 0) time[bid] = clock();
    shared[tid] = 0;
    for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
        i += BLOCK_NUM * THREAD_NUM) {
       shared[tid] += num[i] * num[i];
    }
    __syncthreads();
    if(tid == 0) {
        for(i = 1; i < THREAD_NUM; i++) {
            shared[0] += shared[i];
        }
        result[bid] = shared[0];
    }
    if(tid == 0) time[bid + BLOCK_NUM] = clock();
}

利用 __shared__ 宣告的變數表示這是 shared memory，是一個 block 中每個 thread 都共用的記憶體。它會使用在 GPU 上的記憶體，所以存取的速度相當快，不需要擔心 latency 的問題。

__syncthreads() 是一個 CUDA 的內建函式，表示 block 中所有的 thread 都要同步到這個點，才能繼續執行。在我們的例子中，由於之後要把所有 thread 計算的結果進行加總，所以我們需要確定每個 thread 都已經把結果寫到 shared[tid] 裡面了。

接下來，把 main 函式的一部份改成：

    int* gpudata, *result;
    clock_t* time;
    cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
    cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM);
    cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
    cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
        cudaMemcpyHostToDevice);
    sumOfSquares<<        THREAD_NUM * sizeof(int)>>>(gpudata, result, time);
    int sum[BLOCK_NUM];
    clock_t time_used[BLOCK_NUM * 2];
    cudaMemcpy(&sum, result, sizeof(int) * BLOCK_NUM,
        cudaMemcpyDeviceToHost);
    cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
        cudaMemcpyDeviceToHost);
    cudaFree(gpudata);
    cudaFree(result);
    cudaFree(time);
    int final_sum = 0;
    for(int i = 0; i < BLOCK_NUM; i++) {
        final_sum += sum[i];
    }

可以注意到，現在 CPU 只需要加總 BLOCK_NUM 也就是 32 個數字就可以了。

不過，這個程式由於在 GPU 上多做了一些動作，所以它的效率會比較差一些。在 GeForce 8800GT 上，它需要約 164K cycles。

當然，效率會變差的一個原因是，在這一版的程式中，最後加總的工作，只由每個 block 的 thread 0 來進行，但這並不是最有效率的方法。理論上，把 256 個數字加總的動作，是可以平行化的。最常見的方法，是透過樹狀的加法：

 

把 kernel 改成如下：

__global__ static void sumOfSquares(int *num, int* result,
    clock_t* time)
{
    extern __shared__ int shared[];
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    int i;
    int ffset = 1, mask = 1;
    if(tid == 0) time[bid] = clock();
    shared[tid] = 0;
    for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
        i += BLOCK_NUM * THREAD_NUM) {
       shared[tid] += num[i] * num[i];
    }
    __syncthreads();
    while(offset < THREAD_NUM) {
        if((tid & mask) == 0) {
            shared[tid] += shared[tid + offset];
        }
        offset += offset;
        mask = offset + mask;
        __syncthreads();
    }
    if(tid == 0) {
        result[bid] = shared[0];   
        time[bid + BLOCK_NUM] = clock();
    }
}

後面的 while 迴圈就是進行樹狀加法。main 函式則不需要修改。

這一版的程式，在 GeForce 8800GT 上執行需要的時間，大約是 140K cycles（相當於約 43GB/s），比完全不在 GPU 上進行加總的版本還快！這是因為，在完全不在 GPU 上進行加總的版本，寫入到 global memory 的資料數量很大（8192 個數字），也對效率會有影響。所以，這一版程式不但在 CPU 上的運算需求降低，在 GPU 上也能跑的更快。

進一步改善
上一個版本的樹狀加法是一般的寫法，但是它在 GPU 上執行的時候，會有 share memory 的 bank conflict 的問題（詳情在後面介紹 GPU 架構時會提到）。採用下面的方法，可以避免這個問題：

    ffset = THREAD_NUM / 2;
    while(offset > 0) {
        if(tid < offset) {
            shared[tid] += shared[tid + offset];
        }
        offset >>= 1;
        __syncthreads();
    }

這樣同時也省去了 mask 變數。因此，這個版本的執行的效率就可以再提高一些。在 GeForce 8800GT 上，這個版本執行的時間是約 137K cycles。當然，這時差別已經很小了。如果還要再提高效率，可以把樹狀加法整個展開：

    if(tid < 128) { shared[tid] += shared[tid + 128]; }
    __syncthreads();
    if(tid < 64) { shared[tid] += shared[tid + 64]; }
    __syncthreads();
    if(tid < 32) { shared[tid] += shared[tid + 32]; }
    __syncthreads();
    if(tid < 16) { shared[tid] += shared[tid + 16]; }
    __syncthreads();
    if(tid < 8) { shared[tid] += shared[tid + 8]; }
    __syncthreads();
    if(tid < 4) { shared[tid] += shared[tid + 4]; }
    __syncthreads();
    if(tid < 2) { shared[tid] += shared[tid + 2]; }
    __syncthreads();
    if(tid < 1) { shared[tid] += shared[tid + 1]; }
    __syncthreads();

當然這只適用於 THREAD_NUM 是 256 的情形。這樣可以再省下約 1000 cycles 左右（約 44GB/s）。最後完整的程式檔可以從這裡下載。

 

本文來自CSDN部落格，轉載請標明出處：http://blog.csdn.net/chinacodec/archive/2008/12/23/3583463.aspx

來自 “ ITPUB部落格 ” ，連結：http://blog.itpub.net/22785983/viewspace-619805/，如需轉載，請註明出處，否則將追究法律責任。