GPU基本演算法（Part II）

Scan應用

壓縮（Compact）

Compact實際上是在一組資料中把我們需要的部分挑出來的一種方法，具體步驟如下：第一步對資料進行一個predicate，將我們需要的資料標為true，其他的資料標為false；第二步開闢一個數組與原陣列對應，將prdicate結果為true對應的位置存入1，其他的存入0；第三步，對這個陣列進行exclusive scan，就可以得到這些資料在新陣列中的地址了，我們把這個地址叫做scatter address；第四步就是將輸入元素對映到輸出元素中。

分配（Allocate）

Allocate是類似compact的操作，其中輸出的項數可以動態地從輸入項計算出。例如下面這個例子，輸入的元素個數不固定，然後需要將其對映到輸出，那麼這種情況我們怎麼獲取scatter address呢？與compact不同的就是在做predicate的時候我們不再用0和1表示，而是寫出輸入元素的個數，之後進行exclusive scan，這樣我們就能得到地址的值了。

稀疏矩陣向量乘法（SpMv）

當一個矩陣0的個數遠大於有值位置的個數的時候，我們把這個矩陣稱為稀疏矩陣，對於稀疏矩陣與向量相乘，效率是非常低的，因為會有大量的0*0，而這種計算是沒有意義的。所以我們需要一種更加高效的演算法。表示稀疏矩陣向量的方法叫做壓縮稀疏行（compressed sparse row），我們對於一個稀疏矩陣用三個向量表示：value、column與rowptr。其中value是按順序寫下矩陣中的每個元素，column是記錄下value中的值分別對應哪一列，rowptr記錄下每一行開始元素對應其在value中的索引。

下面讓我們來看怎麼用CSR來進行稀疏矩陣向量乘法，首先我們根據value與rowptr的值生成分段的元素，其次我們根據column來計算元素對應與向量中哪個元素相乘，進行乘法運算得到一個新的結果，在對這個結果進行分段掃描（segmented scan）就可以得到最終的結果了。在這個過程中我們可以看出scan仍是計算最主要的部分，利用平行計算可以大大地加速。

排序（Sort）

奇偶排序（odd-even sort/brick sort）

在並行世界中的演算法跟序列世界中的演算法還是有一些區別的，在並行世界中最簡單的演算法就是奇偶排序（odd-even sort/brick sort）。這種演算法不斷交換相鄰元素的位置最終使元素達到有序，這種相鄰位置元素的交換跟bubble sort很像，因此這種演算法的複雜度並不是很讓人滿意，step complexity為O（n），work complexity為O（n^2）。

歸併排序

如果我們想有什麼序列演算法最適合平行計算的，那無疑就是歸併排序了，但是如果我們仔細想想就會發現歸併排序中也存在並行不是特別好的部分，下面我們來解決如何將這部分用平行計算來表示。對於一個數據量很大的排序（如下圖），我們可以把它分為三個部分：第一個部分是由許多個小的排序構成，我們往往將每一個任務分給它一個執行緒，這個任務就足夠得以解決。第二部分與第一部分有一個臨界值，是shared memory的大小值，因為我們在之前學到如果把需要頻繁訪問的資料放到shared memory中可以提高效率，因此對於第一部分資料量小的時候我們可以把資料全放進shared memory中。而這部分資料我們往往不採用歸併排序，而是有更高效的排序網路；然後到第二部分，這一部分資料量多了起來，一個排序任務的資料量比較多，因此我們將每個任務分給一個block處理，因此這個階段在執行的SM是與任務數相當的；到了第三個階段我們每次排序的資料量非常多，而任務數很少，這種情況我們將一個大任務分成許多個小任務，交由不同的SM處理，這樣子我們就讓每一個SM都有任務可以做，大大提高了GPU的效率。

第二部分

我們先來看第二階段歸併排序的並行版本改進。我們知道對於兩個有序序列要合成一個有序序列需要一個序列處理器，在兩個有序序列的頭部挑一個最小的元素作為輸出元素，如此反覆，知道兩個有序序列全部讀完，這時候我們得到的輸出序列就是一個有序序列。

我們可以看出這種序列演算法很顯然不適合平行計算，對於較長的序列，需要很久的處理時間，那我們怎麼講其分解為並行的任務呢？答案就是我們對於每一個元素都算出它在輸出序列中的地址，然後由輸入到輸出對映即可。那麼如何得到這個地址呢？我們知道這個地址是由當前元素在本序列的位置加上在另外一個序列中按大小排序的位置相加得到的。例如下圖中的12在本序列中的索引是2，而將這個數放到序列2中按大小排序得到的位置索引也是2，因此這個元素在輸出序列中的位置索引是4。那麼這個數值應該如何計算呢？首先這個元素在自己序列的索引很容易知道，那麼我們需要知道的就是它在另一個序列中的位置，計算方法就是把這個元素放到另一個序列中進行二分搜尋，就能得到第二個位置索引了。而二分搜尋的過程就可以採用平行計算了。

第三部分

我們知道當資料量很多，任務數很少的時候，如果我們還是讓一個block負責一個任務的話，就會有很多SM處於空閒狀態，為了解決這個問題我們可以將大任務分解成很多小任務。首先我們在大序列中取出一些元素稱為（splitter），然後可以利用上述的方法對其進行排序，也可以求出當前元素在另外一個序列中的位置，這樣我們就將本來很大序列分為了許多小塊，為了使最終排序的高效性，我們可以將splitter的間隔限制在一定的範圍，以保證最終排序的元素都能放入shared memory中。當我們將序列分為許多區間以後，我們就可以對其進行排序了。例如FC之間，我們已經確定了F與C的位置，那麼輸出序列中FC之間的元素都在下圖的紅色區域中，因此我們只需要對這兩個子序列進行排序就可以得到最終輸出序列FC之前的元素順序。這樣我們就將本來兩個大的有序序列分為了很多個成對的小有序序列，之後我們採用第二部分的方法就可以對子任務進行計算了。

第一部分

在這一部分我們需要考慮與輸入無關的演算法，也就是不管輸入資料是以什麼順序排列的，其計算複雜度是一樣的。我們把這一部分稱為排序網路（Sorting network），而最常用的演算法就是雙調排序，我們把序列先上升後下降（或先下降後上升的序列）稱為雙調序列，我們將雙調序列內部對應元素進行比較，將大的元素放在一個數組中，小的元素放在另一個數組中，我們可以得到兩個雙調序列，而且第一個雙調序列中的元素全部大於第二個雙調序列中的元素。這樣重複進行最終就能得到排序完成的序列了。而我們剛才談的是對雙調序列的排序，那麼對於任意一個普通的序列我們如何生成它的雙調序列的，過程就是上述的逆過程，從相鄰的兩個元素開始生成小的雙調序列，然後對於小雙調序列排序，得到有序的小序列，將兩個小的有序序列以相反的方向拼接在一起就得到一個大的雙調序列了，這樣重複進行最終就能得到原序列對應的雙調序列了。之前浙大考核的時候我也做過雙調序列更詳細的學習，具體內容：戳我。對於這個過程在GPU的實現，我們只需要給每個元素分配一個執行緒來儲存比較後元素的值，這樣每次進行一組比較後我們進行一次同步就可以開始下一次比較了。

如果能將輸入元素全都放進shared memory中，那麼排序網路是一種非常有效的方法，需要注意的是雙調排序並不是排序網路的唯一演算法，還有奇偶歸併排序（odd-even merge sort）等演算法也可以用作排序網路。

基數排序

最後要講的是在GPU中效率最高的方法：基數排序。到目前為止我們講的排序方法都是比較排序，也就是交換元素的順序，而比較排序是依賴於數字位置的比較方法。我們將一個數用二進位制位來表示，從最低位到最高位開始掃描，把每一位是0的放在上面，是1的放在下面，之後開始掃描第二位，重複這個過程直到最高位掃描結束，得到的序列就是排序好的序列。這個演算法在GPU上執行流暢的原因有兩個，第一個就是其優越的複雜度，為O（kn），k是表示一個數的位元位。第二個就是每一次排序過程都可以用我們之前學過的演算法表示，也就是compact，我們對一個特定的位元位進行掃描就能得到元素的在新序列中的地址索引。通過每次掃描多個位元位我們還能提高這個演算法的效率。

快速排序

快速排序也是我們在序列世界中用的很多的演算法，我們選用一個參考元素（pivot element），之後把元素分為小於該元素，等於該元素，大於該元素的三個陣列，再分別對這三個陣列用同樣的方法重複進行，直到排序完畢。

在並行實現的時候我們可以通過distribute、map、compact的操作將原陣列分段，然後再段內開啟另外的執行緒來達到與遞迴類似的效果。