文章目錄
- 原子操作簡介
- 計算直方圖
-
- 在 CPU 上計算直方圖
- 在GPU上計算直方圖
-
- 使用全局記憶體原子操作的直方圖核函數
- 使用共享記憶體原子操作和全局記憶體原子操作的直方圖核函數
- 更多相關内容
原子操作簡介
在編寫傳統的單線程應用程式時，程式員通常不需要使用原子操作。即使需要使用原子操作也無需擔心。下面将詳細解釋原子操作是什麼，以及為什麼在多線程程式中需要使用它們。為了解釋原子操作，首先來分析 C 或者 C++ 的基礎知識之一，即遞增運算符：
這是在标準 C 中的一個表達式，在執行完這條語句後，
x
的值應該比執行這條語句之前的值大 1。這條語句中包含如下三個步驟：
- 讀取
x
- 将步驟 1 中讀到的值增加 1
- 将遞增後的結果寫回到
x
有時候，這個過程也稱為讀取-修改-寫入(Read-Modify-Write)操作，其中第 2 步的遞增操作也可以換成其他修改
x
值的操作。
現在考慮一種情況：有兩個線程都需要對
x
的值進行遞增。将這兩個線程分别稱為 A 和 B。A 和 B 在遞增
x
的值時都需要執行上面三個操作。假設
x
的初始值為 7。在理想情況下，我們希望線程 A 和 B 執行下表中的步驟
| 步驟 | 示例 |
|---|---|
1) 線程 A 讀取 | A 從 |
| 2) 線程 A 将讀到的值增加 1 | A 計算得到 8 |
3) 線程 A 将結果寫回到 | |
4) 線程 B 讀取 | B 從 |
| 5) 線程 B 将讀到的值增加 1 | B 計算得到 9 |
6) 線程 B 将結果寫回到 | |
由于
x
的起始值為 7，并且由兩個線程進行遞增，是以在遞增運算完成後，
x
的值變為 9。根據前面的操作順序，這确實是文明得到的結果。遺憾的是，除了這個操作順序外，還存在其他一些操作順序可能導緻錯誤的結果。例如，下表中的順序，其中線程 A 和 B 的操作彼此交叉進行。
| 步驟 | 示例 |
|---|---|
1) 線程 A 讀取 | A 從 |
2) 線程 B 讀取 | B 從 |
| 3) 線程 A 将讀到的值增加 1 | A 計算得到 8 |
| 4) 線程 B 将讀到的值增加 1 | B 計算得到 8 |
5) 線程 A 将結果寫回到 | |
6) 線程 B 将結果寫回到 | |
是以，如果線程的排程方式不正确，那麼最終将得到錯誤的結果。除了上面兩種執行順序外，這 6 個步驟還有許多其他的排序方式，其中有些方式能産生正确的結果，而其他的方式則不能。當把程式從單線程改寫為多線程時，如果多個線程需要對共享值進行讀取或者寫入時，那麼很可能會遇到不可預測的結果。
在前面的示例中，我們需要通過某種方式一次性地執行完讀取-修改-寫入這三個操作，并且在執行過程中不會被其他線程中斷。具體來說，除非已經完成了這三個操作，否則其他的線程都不能讀取或者寫入
x
的值。由于這些操作的執行過程不能分解為更小的部分，是以我們将滿足這種條件限制的操作稱為原子操作。CUDA C 支援多種原子操作，當有數千個線程在記憶體通路上發生競争時，這些操作能夠確定在記憶體上實作安全的操作。
計算直方圖
給定一個包含一組元素的資料集，直方圖表示每個元素出現的頻率。例如，如果計算詞組 “Programming with CUDA C” 中字元頻率的直方圖，那麼将得到下表結果。
| A | C | D | G | H | I | M | N | O | P | R | T | U | W |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 2 | 1 | 2 | 1 | 2 | 2 | 1 | 1 | 1 | 2 | 1 | 1 | 1 |
雖然直方圖的定義很簡單，但卻在計算機科學領域得到了非常廣的應用。在各種算法中都用到直方圖，包括圖像處理、資料壓縮、計算機視覺、機器學習、音頻編碼等等。下面将把直方圖運算作為代碼示例的算法。
在 CPU 上計算直方圖
下面首先給出如何在CPU上計算直方圖。這個示例同時也說明了在單線程應用中計算直方圖是非常簡單的。這個應用程式将處理一個大型的資料流。在實際程式中，這個資料可以是像素的顔色值，或者音頻采樣資料，但在我們的示例程式中，這個資料是随機生成的位元組流。我們可以通過工具函數
big_random_block()
來生成這個随機的位元組流。在應用程式中将生成 100MB 的随機資料。
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#include <windows.h>
#define SIZE (100*1024*1024)
int main(void) {
unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
由于每個随機位元組 (8比特) 都有 256 個不同的可能取值 （從0x00到0xFF），是以在直方圖中需要包含 256 個元素，每個元素記錄相應的值在資料流中出現次數。
下面建立一個包含 256 個元素的數組，并将所有元素的值初始化為 0。
unsigned int histo[256];
for (int i = 0; i < 256; i++)
histo[i] = 0;
在建立了直方圖并将數組元素初始化為 0 後，接下來需要計算每個值在
buffer[]
資料中的出現頻率。
算法的思想是，每當在數組
buffer[]
中出現某個值
z
時，就遞增直方圖數組中索引為
z
的元素。這樣就能計算出值
z
的出現次數。
如果目前看到的值為
buffer[i]
，那麼将遞增數組中索引等于
buffer[i]
的元素。由于元素
buffer[i]
位于
histo[buffer[i]]
，我們隻需一行代碼就可以遞增相應的計數器。
我們在簡單的
for()
循環中對
buffer[]
每個元素執行這個操作。
for (int i = 0; i < SIZE; i++)
histo[buffer[i]]++;
此時，我們已經計算完了輸入資料的直方圖。在實際的應用程式中，這個直方圖可能作為下一個計算步驟的輸入資料。但在這裡的簡單示例中，這就是要執行的所有工作，是以接下來将驗證直方圖的所有元素相加起來是否等于正确的值，然後結束程式。
long histoCount = 0;
for (int i = 0; i < 256; i++) {
histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;
思考一下就會發現，無論輸入數組的值是什麼，這個和值總是相同的。每個元素都将統計相應數值的出現次數，是以所有這些元素值的總和就應該等于數組中元素的總數量。在示例中，這個值就等于
SIZE
。
在執行完運算後需要釋放記憶體并傳回。
free(buffer);
return 0;
}
完整代碼
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#include <windows.h>
#define SIZE (100*1024*1024)
int main(void) {
unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
DWORD start, end;
start = GetTickCount();
unsigned int histo[256];
for (int i = 0; i < 256; i++)
histo[i] = 0;
for (int i = 0; i < SIZE; i++)
histo[buffer[i]]++;
end = GetTickCount();
std::cout << "Time to generate: " << end - start << " ms" << std::endl;
long histoCount = 0;
for (int i = 0; i < 256; i++) {
histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;
free(buffer);
return 0;
}
運作結果
在GPU上計算直方圖
我們把這個直方圖計算示例改在 GPU 上運作。如果輸入的數組足夠大，那麼通過由多個線程來處理緩沖區的不同部分，将節約大量的計算時間。其中，由不同的線程來讀取不同部分的輸入資料是非常容易的。
在計算輸入數組的直方圖時存在一個問題，即多個線程可能同時對輸出直方圖的同一個元素進行遞增。在這種情況下，我們需要通過原子的遞增操作來避免之前提到線程執行順序不同帶來的問題。
main()
函數的開頭與基于 CPU 的版本完全一樣
int main( void ) {
unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
由于要測量代碼的執行性能，是以要初始化計時事件。
cudaEvent_t start, stop;
HANDLE_ERROR(cudaEventCreate(&start));
HANDLE_ERROR(cudaEventCreate(&stop));
HANDLE_ERROR(cudaEventRecord(start, 0));
在設定好輸入資料和事件後，我們需要在 GPU 上為随機輸入資料和輸出直方圖配置設定記憶體空間。在配置設定了輸入緩沖區後，我們将
big_random_block()
生成的數組複制到 GPU 上。同樣，在配置設定了直方圖後，像 CPU 版本中那樣将其初始化為 0。
// 在GPU上為檔案的資料配置設定記憶體
unsigned char* dev_buffer;
unsigned int* dev_histo;
HANDLE_ERROR(cudaMalloc((void**)&dev_buffer, SIZE));
HANDLE_ERROR(cudaMemcpy(dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMalloc((void**)&dev_histo, 256 * sizeof(int)));
HANDLE_ERROR(cudaMemset(dev_histo, 0, 256 * sizeof(int)));
cudaMemset()
這個函數的原型與标準 C 函數
memset()
的原型是相似的，并且這兩個函數的行為也基本相同。二者的差異在于，
cudaMemset()
将傳回一個錯誤碼，而 C 庫函數
memset()
則不是。這個錯誤碼将告訴調用者在設定 GPU 記憶體時發生的錯誤。除了傳回錯誤碼外，還有一個不同之處就是，
cudaMemset()
是在 GPU 記憶體上執行，而
memset()
是在主機記憶體上運作。
在初始化輸入緩沖區和輸出緩沖區後，就做好了計算直方圖的準備。你馬上就會看到如何準備并啟動直方圖核函數。我們暫時假設已經在 GPU 上計算好了直方圖。在計算完成後，需要将直方圖複制回 CPU，是以我們配置設定了一個包含 256 個元素的數組，并且執行從裝置到主機的複制。
unsigned int histo[256];
HANDLE_ERROR(cudaMemcpy(histo, dev_histo, 2
此時，我們完成了直方圖的計算，是以可以停止計時器并顯示經曆的時間。
// 得到停止時間并顯示計時結果
HANDLE_ERROR(cudaEventRecord(stop, 0));
HANDLE_ERROR(cudaEventSynchronize(stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
std::cout << "Time to generate: " << elapsedTime << " ms\n";
此時，我們可以将直方圖作為輸入資料傳遞給算法的下一個步驟。然而，在這個示例中不需要将直方圖用于其他任何操作，而隻是驗證在 GPU 上計算得到的直方圖與在 CPU 上計算得到的直方圖是否相同。首先，我們驗證直方圖的總和等于正确的值。這與 CPU 版本中的代碼是相同的，如下所示：
long histoCount = 0;
for (int i = 0; i < 256; i++) {
histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;
計算出 GPU 直方圖，并在周遊每個數值時，遞減直方圖中相應元素的值。是以，如果完成計算時直方圖每個元素的值都為 0，那麼 CPU 計算的直方圖與 GPU 計算的直方圖相等。從某種意義上來說，我們是在計算這兩個直方圖之間的差異。
// 驗證與基于CPU計算得到的結果是相同的
for (int i = 0; i < SIZE; i++)
histo[buffer[i]]--;
for (int i = 0; i < 256; i++) {
if (histo[i] != 0)
std::cout << "Failure at " << i << "!\n";
}
程式結束前要施放已配置設定的 CUDA 事件，GPU 記憶體和主機記憶體。
HANDLE_ERROR(cudaEventDestroy(start));
HANDLE_ERROR(cudaEventDestroy(stop));
cudaFree(dev_histo);
cudaFree(dev_buffer);
free(buffer);
return 0;
}
由于直方圖包含了 256 個元素，是以可以在每個線程塊中包含 256 個線程，這種方式不僅友善而且高效。但是，線上程塊的數量上還可以有更多選擇。例如，在 100MB 資料中共有 104857600 個位元組。我們可以啟動一個線程塊，并且讓每個線程處理 409600 個資料元素。同樣，我們還可以啟動 409600 個線程塊，并且讓每個線程處理一個資料元素。
最優的解決方案是在這兩種極端情況之間。通過一些性能實驗，我們發現當線程塊的數量為 GPU 中處理器數量的 2 倍時，将達到最優性能。例如，在 GeForce GTX280 中包含了 30 個處理器，是以當啟動 60 個并行線程塊時，直方圖核函數将運作得最快。
如果要基于目前的硬體平台來動态調整線程塊的數量，那麼就要用到其中一個裝置屬性。我們通過以下代碼片段來實作這個操作。
cudaDeviceProp prop;
HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
int blocks = prop.multiProcessorCount;
std::cout << "MultiProcessorCount: " << blocks << std::endl;
histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);
使用全局記憶體原子操作的直方圖核函數
計算直方圖的核函數需要的參數包括：
- 一個指向輸入數組的指針
- 輸入數組的長度
- 一個指向輸出直方圖的指針
核函數執行的第一個計算就是計算輸入資料數組中的偏移。每個線程的起始偏移都是 0 到線程數量減 1 之間的某個值。然後，對偏移的增量為已啟動線程的總數。
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)
// 使用全局記憶體原子操作的直方圖核函數
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
int i = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
while (i < size) {
atomicAdd(&(histo[buffer[i]]), 1);
i += stride;
}
}
函數調用
atomicAdd(addr, y)
将生成一個原子的操作序列，這個操作序列包括讀取位址
addr
處的值，将
y
增加到這個值，以及将結果儲存回位址
addr
。
底層硬體将確定目前執行這些操作時，其他任何線程都不會讀取或寫入位址
addr
上的值，這樣就能確定得到預計的結果。
在這裡的示例中，這個位址就是直方圖中相應元素的位置。如果目前位元組為
buffer[i]
，那麼直方圖中相應的元素就是
histo[buffer[i]]
。原子操作需要這個元素的位址，是以第一個參數為
&(histo[buffer[i])
。由于我們隻是想把這個元素中的值遞增 1，是以第二個參數就是 1。
完整代碼
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)
//使用全局記憶體原子操作的直方圖核函數
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
int i = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
while (i < size) {
atomicAdd(&(histo[buffer[i]]), 1);
i += stride;
}
}
int main( void ) {
unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
cudaEvent_t start, stop;
HANDLE_ERROR(cudaEventCreate(&start));
HANDLE_ERROR(cudaEventCreate(&stop));
HANDLE_ERROR(cudaEventRecord(start, 0));
// 在GPU上為檔案的資料配置設定記憶體
unsigned char* dev_buffer;
unsigned int* dev_histo;
HANDLE_ERROR(cudaMalloc((void**)&dev_buffer, SIZE));
HANDLE_ERROR(cudaMemcpy(dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMalloc((void**)&dev_histo, 256 * sizeof(int)));
HANDLE_ERROR(cudaMemset(dev_histo, 0, 256 * sizeof(int)));
cudaDeviceProp prop;
HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
int blocks = prop.multiProcessorCount;
std::cout << "MultiProcessorCount: " << blocks << std::endl;
histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);
unsigned int histo[256];
HANDLE_ERROR(cudaMemcpy(histo, dev_histo, 256*sizeof(int), cudaMemcpyDeviceToHost));
// 得到停止時間并顯示計時結果
HANDLE_ERROR(cudaEventRecord(stop, 0));
HANDLE_ERROR(cudaEventSynchronize(stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
std::cout << "Time to generate: " << elapsedTime << " ms\n";
long histoCount = 0;
for (int i = 0; i < 256; i++) {
histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;
// 驗證與基于CPU計算得到的結果是相同的
for (int i = 0; i < SIZE; i++)
histo[buffer[i]]--;
for (int i = 0; i < 256; i++) {
if (histo[i] != 0)
std::cout << "Failure at " << i << "!\n";
}
HANDLE_ERROR(cudaEventDestroy(start));
HANDLE_ERROR(cudaEventDestroy(stop));
cudaFree(dev_histo);
cudaFree(dev_buffer);
free(buffer);
return 0;
}
運作結果
使用共享記憶體原子操作和全局記憶體原子操作的直方圖核函數
由于在核函數中隻包含了非少的計算工作，是以很可能全局記憶體上的原子操作導緻性能的降低。
當數千個線程嘗試通路少量的記憶體位置時，将發生大量的競争。為了確定遞增操作的原子性，對相同記憶體位置的操作都将被硬體串行化。這可能導緻未完成操作的隊列非常長，是以會抵消通過并行運作線程而獲得的性能提升。
盡管這些原子操作是導緻這種性能降低的原因，但解決這個問題的方法卻出乎意料地需要使用更多而非更少的原子操作。
這裡的主要問題并非在于使用了過多的原子操作，而是有數千個線程在少量的記憶體位址上發生競争。要解決這個問題，我們将直方圖計算分為兩個階段。
在第一個階段，每個并行線程塊将計算它所處理資料的直方圖。由于每個線程塊在執行這個操作時都是互相獨立的，是以可以在共享記憶體中計算這些直方圖，這将避免每次将寫入操作從晶片發送到DRAM。但是，這種方式仍然需要原子操作，因為線上程塊中多個線程之間仍然會處理相同值的資料元素。然而，現在隻有 256 個線程在 256 個位址上發生競争，這将極大地減少在使用全局記憶體時在數千個線程之間發生競争的情況。
然後，在第一個階段中配置設定一個共享記憶體緩沖區并進行初始化，用來儲存每個線程塊的臨時直方圖。由于随後的步驟将包括讀取和修改這個緩沖區，是以需要調用
__syncthreads()
來確定每個線程的寫入操作線上程繼續前進之前完成。
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
在将直方圖初始化為 0 後，下一個步驟與最初 GPU 版本的直方圖計算非常類似。這裡唯一的差異在于，我們使用了共享記憶體緩沖區
temp[]
而不是全局記憶體緩沖區
histo[]
，并且需要随後調用
__syncthreads()
來確定送出最後的寫入操作。
int i = threadIdx.x + blockIdx.x * blockDim.x;
int offset = blockDim.x * gridDim.x;
while (i < size) {
atomicAdd(&temp[buffer[i]], 1);
i += offset;
}
__syncthreads();
最後一步要求将每個線程塊的臨時直方圖合并到全局緩沖區
histo[]
中。假設将輸入資料分為兩半，這樣就有兩個線程檢視不同部分的資料，并計算得到兩個獨立的直方圖。如果線程 A 在輸入資料中發現位元組 0xFC 出現了 20 次，線程 B 發現位元組 0xFC 出現了 5 次，那麼位元組 0xFC 在輸入資料中共出現了 25 次。同樣，最終直方圖的每個元素隻是線程 A 直方圖中相應元素和線程 B 直方圖中相應元素的加和。
這個邏輯可以擴充到任意數量的線程，是以将每個線程塊的直方圖合并為單個最終的直方圖就是，将線程塊的直方圖的每個元素都相加到最終直方圖中相應位置的元素上。這個操作需要自動完成：
atomicAdd(&(histo[threadIdx.x]), temp[threadIdx.x]);
}
由于我們使用了 256 個線程，并且直方圖中包含了 256 個元素，是以每個線程都将自動把它計算得到的元素隻增加到最終直方圖的元素上。如果線程數量不等于元素數量，那麼這個階段将更為複雜。
注意，我們并不保證線程塊将按照何種順序将各自的值相加到最終直方圖中，但由于整數加法是可交換的，無論哪種順序都會得到相同的結果。
完整代碼
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../../common/book.h"
#include <stdio.h>
#include <iostream>
#define SIZE (100*1024*1024)
// 使用共享記憶體原子操作和全局記憶體原子操作的直方圖核函數
__global__ void histo_kernel(unsigned char* buffer, long size, unsigned int* histo) {
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
int i = threadIdx.x + blockIdx.x * blockDim.x;
int offset = blockDim.x * gridDim.x;
while (i < size) {
atomicAdd(&temp[buffer[i]], 1);
i += offset;
}
__syncthreads();
atomicAdd(&(histo[threadIdx.x]), temp[threadIdx.x]);
}
int main( void ) {
unsigned char* buffer = (unsigned char*)big_random_block(SIZE);
cudaEvent_t start, stop;
HANDLE_ERROR(cudaEventCreate(&start));
HANDLE_ERROR(cudaEventCreate(&stop));
HANDLE_ERROR(cudaEventRecord(start, 0));
// 在GPU上為檔案的資料配置設定記憶體
unsigned char* dev_buffer;
unsigned int* dev_histo;
HANDLE_ERROR(cudaMalloc((void**)&dev_buffer, SIZE));
HANDLE_ERROR(cudaMemcpy(dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMalloc((void**)&dev_histo, 256 * sizeof(int)));
HANDLE_ERROR(cudaMemset(dev_histo, 0, 256 * sizeof(int)));
cudaDeviceProp prop;
HANDLE_ERROR(cudaGetDeviceProperties(&prop, 0));
int blocks = prop.multiProcessorCount;
std::cout << "MultiProcessorCount: " << blocks << std::endl;
histo_kernel << <blocks * 2, 256 >> > (dev_buffer, SIZE, dev_histo);
unsigned int histo[256];
HANDLE_ERROR(cudaMemcpy(histo, dev_histo, 256*sizeof(int), cudaMemcpyDeviceToHost));
// 得到停止時間并顯示計時結果
HANDLE_ERROR(cudaEventRecord(stop, 0));
HANDLE_ERROR(cudaEventSynchronize(stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
std::cout << "Time to generate: " << elapsedTime << " ms\n";
long histoCount = 0;
for (int i = 0; i < 256; i++) {
histoCount += histo[i];
}
std::cout << "Histogram Sum: " << histoCount << std::endl;
// 驗證與基于CPU計算得到的結果是相同的
for (int i = 0; i < SIZE; i++)
histo[buffer[i]]--;
for (int i = 0; i < 256; i++) {
if (histo[i] != 0)
std::cout << "Failure at " << i << "!\n";
}
HANDLE_ERROR(cudaEventDestroy(start));
HANDLE_ERROR(cudaEventDestroy(stop));
cudaFree(dev_histo);
cudaFree(dev_buffer);
free(buffer);
return 0;
}
運作結果
可以看到相比于僅僅使用全局記憶體原子操作的版本要快了 1 倍左右。
更多相關内容
CUDA atomic原子操作 —— -牧野-