本周主要內容如下：
- 如何分析GPU算法的速度和效率
- 3個新的基本算法：歸約、掃描和直方圖(Reduce、Scan、Histogram)
一、評估標准
首先介紹用於評估GPU計算的兩個標准：
- step ：完成某特定計算所需時間--挖洞操作（Operation Hole Digging）
- work：工作總量
如下圖示，第一種情況只有一個工人挖洞，他需要8小時才能完成，所以工作總量（Work）是8小時。第二種情況是有4個工人，它們2個小時就能完成挖洞任務，此時工作總量是8小時。第三種情況同理不加贅述。
另一種直觀理解如下圖示。該計算任務步驟復雜度是3，工作復雜度是7。而接下來的課程的目的則是學會如何優化GPU算法。
二、3個新的基本算法
2.1 Reduce
2.1.1 Reduce運算基本介紹
下圖展示的是reduce運算。
reduce運算由兩部分組成：輸入值和 歸約運算符（Reduction Operation）。
其中歸約運算符需要滿足下面兩個條件：
- Binary 二元運算符：即對兩個輸入對象進行操作得到一個結果
- Associative 結合性：直白解釋就是運算符需要與順序無關，例如加號就滿足這一性質，而減號不滿足。
2.1.2 串行方式實現歸約和並行方式實現歸約
- 串行方式實現歸約
reduce運算和map運算有點類似，但不同的地方在於map是可以對多個元素做單獨操作的，而reduce每次計算都需要依賴上一次運算的結果。示意圖如下：
- 並行方式實現歸約
並行方式與串行方式的不同點在於計算的方式不同，如果以穿行方式對a,b,c,d進行加法歸約運算，那么計算方式則為\(((a+b)+c)+d\)。而如果以並行方式的話，則為\((a+b)+(c+d)\)。
兩種方式的復雜度如下圖示：
上面的例子貌似並不能很明顯的看出串行和並行的差異，下面我們通過大致計算來更直觀的看出不同：
首先假設使用的並行的方式計算，當輸入對象有兩個，那么step=1。當有4個輸入對象，那么就有step=2。以此類推，當有N個輸入對象，那么step=logN。而串行的step=N-1。因此當N一定大的時候，並行方式的優勢就很明顯了。
下面通過具體示例來看看使用 Global Memory和 Shared Memory進行並行歸約運算的區別。
如下圖示，一共有1024個線程塊，其中每個線程塊有1024個線程需要進行歸約運算。運算分成兩個步驟，首先塊內運算得到1個值，之后再將1024個值再做歸約運算。
下面給出了使用全局變量的代碼。實現思路是將線程分成兩份，然后等距相加。
具體來說就是首先將后512個線程值添加對應添加到前512個線程中去，之后再折半添加，即將前面512個線程拆成前256和后256個線程，以此類推，直到最后計算得到一個值。
下面給出了使用 Shared Memory的代碼，它的思路與上面類似，但是因為使用共享內存，速度得到了提升
理論上來說使用Global Memory需要進行2n次讀操作，n次寫操作，而Shared Memory需要進行n次讀操作，1次寫操作，所以后者應該比前者快3倍左右。
但是因為該示例並沒有使得計算飽和，所以最終測試結果並沒有3倍，具體時間效果如下圖。
2.2 Scan
2.2.1 Scan基本介紹
下圖給出了Scan的例子，輸入為1,2,3,4，運算符是+，輸出結果是當前輸入值與前面值的總和。
咋看貌似並不像是並行計算，但是Scan運算對於並行計算具有很大的作用。
下圖給出了Scan的在實際生活中的例子，即銀行存款賬戶余額情況，左邊表示存錢，取錢數，右邊表示余額。
2.2.2 Scan運算組成
Scan運算由輸入向量，二元結合運算符(與Reduce運算類似)以及標識值組成。
基本上該課程中提到的運算符都需要具有Associative(結合性)，這樣更加符合並行計算的特點。
為方便說明，標識值(Identify element)用I表示，它需要滿足\(I\,\,op \,\,a = a\),其中op表示二元結合運算符，a表示任意值。舉例來說可能更好理解：
例如如果運算符是and，那么標識值就是1。如果運算符是or，那么標識值就是0。
2.2.3 Exclusive Scan & Inclusive Scan
Scan運算分為兩種：Exclusive和 Inclusive。區別是前者的輸出不包含當前的輸入值，后者則包括。下圖給出了具體例子，使用的運算符是add。
下面以Inclusive Scan為例來計算Step復雜度和Work復雜度。
- Step復雜度
它的Steps復雜度和Reduce運算相同都是\(O(logn)\)，不再詳細解釋。
- Work復雜度
Work復雜度是\(o(n^2)\)。如下圖示，第一個輸出值的計算量是0，第二個是1，第三個是2，以此類推，最終復雜度是\(0+1+2+...+(n-1)≈\frac{n*(n-1+0)}{2}=O(n^2)\)
2.2.4 改進Scan運算的實現方法
如果采用上面的實現方式，Work復雜度太大，在實際運用中消耗太多資源，為了改善這一問題，出現了如 Hillis Steele和 Blelloch Scan方法，下面分別進行介紹。
- 1. Hillis Steele Inclusive Scan
還是從1到8做scan運算。該方法的主要實現思路是對於step i，每個位置的輸出值是當前值加上其\(2^i\)左邊的值。
Step 0：每個輸出值是當前值加上前一個值。（綠色）
Step 1：每個輸出值是當前值加上前兩個值（藍色）。
Step 2：每個輸出值是當前值加上前4 (\(=2^2\)) 個值（紅色）。
總結起來計算過程就是：
step 0: out[i] = in[i] + in[i-\(2^0\)];
step 1: out[i] = in[i] + in[i-\(2^1\)];
step 2: out[i] = in[i] + in[i-\(2^2\)];
通過下面完整的計算過程，我們可以明顯看到Step復雜度依舊是\(logn\).
Work復雜度不能很明顯的算出來，我們可以把所有的計算過程看成一個矩陣，縱向長度是\(logn\),橫向長度是\(n\),所以Work復雜度是\(nlogn\)，該算法相對於之前的\(o(n^2)\)Work復雜度有了一定的提升。
- 2. Blelloch Scan
這是一種優化Work complexity的算法，比上面的要復雜一些。主要分為reduce和downsweep兩步。
該算法由兩部分組成：Reduce+Downsweep。
以exclusive scan加法運算為例，如下圖示。
為更好說明用下圖進行解釋說明(圖片來源：CUDA系列學習（五）GPU基礎算法: Reduce, Scan, Histogram
如圖示，上部分是Reduce，不再贅述。
下部分的Downsweep其實可以理解成Reduce的鏡像操作，對於每一組輸入in1 & in2有兩個輸出，左邊輸出out1 = in2，右邊輸出out2 = in1 op in2 （這里的op就是reduce部分的op），如圖：
該算法有點繞，是否明白了呢？下面做個題來看看是不是真的明白了：
運算符是Max，將框中答案補充完整。
Ans：
介紹了這么多，那該方法的復雜度如何呢？
因為我們已經知道Reduce的Step復雜度是\(O(logn)\),Work復雜度是\(O(n)\)。
而Downsweep其實可以理解成Reduce的鏡像運算，所以復雜度與之相同。
所以該算法整體的Step復雜度是\(O(2logn)\),Work復雜度是\(O(n)\)。這里對於Work復雜度存疑，不明白為什么不是\(O(2n)\)，但是網上好幾篇博文都說是\(O(n)\)）
2.2.5 如何選擇合適的算法
上面已經分別介紹了 Hillis Steele和Blelloch Scan算法，下面將二者的復雜度總結如下：
| Method | Step complexity | Work Complexity | type |
|---|---|---|---|
| Hillis Steele | \(O(logn)\) | \(O(nlogn)\) | SE (Step Effient) |
| Blelloch | \(O(2logn)\) | \(O(n)\) | WE (Work Efficient) |
在實際中可能會遇到如下圖示的情況，隨着工作量的變化我們需要動態改變算法策略。
例如最開始可能任務中地Work要遠多於已有的處理器數量，這種情況執行速度受到了處理器數量的限制，此時我們則需要選擇 更具高效工作性的算法(如Blelloch)。
隨着工作量不斷完成，漸漸地可用處理器數量比工作更多，此時我們則可以選擇 更高效步驟的算法實現(如Hillis Steele算法)。
下面來做個題看是否已經理解了如何選擇合適的算法。
| 情況 | Serial(串行) | Hillis Steele | Blelloch |
|---|---|---|---|
| 128k的vector， 1個processor | √ | ||
| 512個元素的vector， 512個processor | √ | ||
| 一百萬的vector， 512個processor | √ |
解析：
- 第一種情況只有一個處理器，難不成還能選擇並行算法？
- 第二種情況處理器數量剛好合適，工作量也不大，所以可以選擇Hillis Steele
- 第三種情況工作量遠大於處理器數量，所以選擇Blelloch算法來提高工作效率。
2.3 Histogram
2.3.1 Histogram是啥
顧名思義，統計直方圖就是將一個統計量在直方圖中顯示出來。
2.3.2 Histogram 的 Serial 實現
分兩部分：1. 初始化，2. 統計
for(i = 0; i < bin.count; i++)
res[i] = 0;
for(i = 0; i<nElements; i++)
res[computeBin(i)] ++; // computeBin(i)函數是用來判斷第i個元素屬於哪一類
2.3.3 Histogram 的 Parallel 實現
直接實現：
__global__ void naive_histo(int* d_bins, const int* d_in, const in BIN_COUNT){
int myID = threadIdx.x + blockDim.x * blockIdx.x;
int myItem = d_in[myID];
int myBin = myItem % BIN_COUNT;
d_bins[myBin]++;
}
仔細看上面的代碼，我們很容易看出d_bins[myBin]++;語句在並行計算過程中會出現read-modify-write（從全局內存中讀取數據，修改數據，寫回數據） 沖突問題，進而造成結果出錯。
而serial implementation不會有這個問題，那么想實現parallel histogram計算有什么方法呢？
1. accumulate using atomics
下圖給出了read-modify-write介紹，可以看出讀取，修改，寫入是3個獨立的原子運算，但是如果我們將這3個操作整合成1個原子操作那么就可以很好地解決上述問題。而且現如今的GPU能夠鎖定特定的內存地址，因此其他的線程就無法訪問該地址。
具體實現代碼只需要將 d_bins[myBin]++;修改成atomicAdd(&(d_bins[myBin]), 1); 即可。
__global__ void simple_histo(int* d_bins, const int* d_in, const in BIN_COUNT){
int myID = threadIdx.x + blockDim.x * blockIdx.x;
int myItem = d_in[myID];
int myBin = myItem % BIN_COUNT;
atomicAdd(&(d_bins[myBin]), 1);
}
但是對於atomics的方法而言，不管GPU多好，並行線程數都被限制到histogram個數N，也就是最多只有N個線程並行。
2. local memory + reduce
思路原理：設置n個並行線程，每個線程都有自己的local histogram（一個長為bin數的vector）；即每個local histogram都被一個thread順序訪問，所以這樣沒有shared memory，即便沒有用atomics也不會出現read-modify-write問題。然后，我們將這n個histogram進行合並（即加和），可以通過reduce實現。
舉例：
如下圖示，假設有128個item，有8個線程可以使用，然后需要分成3類。
這樣每個線程理論上需要處理16個item，即將這16個item分成3類。
通過使用並行線程，我們無需使用上面那個方法將read-modify-write原子化，因為這8個線程使用的是自己的local memory，而不是shared memory。如此一來我們只需要在所有線程計算結束后，再使用Reduce運算將8個線程的計算結果對應相加即可。
3. sort then reduce by key
該方法的主要思路是采用鍵值對(類似Python中的字典)的方式來記錄數據，之后對鍵值對按鍵的大小排序得到下圖中藍色鍵值對。
如此一來鍵相同的挨在一起，之后便可用Reduce運算將相同類別的值相加即可。
3.作業應用
Tone Mapping(色調映射) 是轉換一組顏色到另一組顏色的過程。之所以要色調映射是因為真實世界中的光譜強度是非常廣的，例如白天太陽的光亮程度到晚上月亮的光亮程度，這之間的范圍十分大。而我們的設備，如電腦顯示屏，手機所能顯示的光亮程度遠小於真實世界，所以我們需要進行色調映射。我們手機上拍照功能中的HDR模式就是色調映射的作用。
該過程如果沒有處理好，那么得到的效果則要么變得過於暗淡，要么過於明亮。
下圖很好地展示了這兩種情況。
下圖左曝光較低，很多細節丟失。下圖右雖然細節保留更多，但是窗戶部分因為曝光過多基本模糊掉了。兩圖片下面的直方圖也很直觀的呈現了像素分布情況。
而色調映射就是為了解決上述情況，以期望達到如下圖的效果。
本次作業需要結合所學的三種運算Reduce，Scan，Histogram。
