# Natural merge sort 在特定硬體的加速
contributed by <`petermouse`>, <`heathcliffYang`>, <`janetwei`>
程式碼: [natural-mergesort](https://github.com/petermouse/natural-mergesort)
[Youtube](https://youtu.be/NiAH4nSldcs)
## 目標
- 改良 Natural merge sort 的 code，注意 branch prediction 的影響
- 嘗試使用硬體提供的 clz 指令來改寫
## 討論
### `__builtin_clz` 與 `LZCNT` ？
我們看到 x86 有提供組語的指令 `LZCNT` 可以算 clz，但是手冊上提供的 C 版函式需要用 intel 的 compiler， gcc 並不提供，後來又發現給 gcc 的函式`__builtin_clz`，因此測試 gcc 的 compiler 在編譯的時候，到底會不會編成`LZCNT`?
> 從 gcc 的參數`-g -Wa,-adhl` 可以對照組語
>> 用 `objdump -d` 觀察 [name=jserv]
>> OK！謝謝老師！ [name=heathcliffYang]
### `__builtin_clz` 與 `LZCNT` 測試
`objdump` 可以顯示有關於目的檔的資訊，其中的 `-d` 參數代表對檔案反組譯( disassemble )
因此我寫了個簡單的程式測試 gcc 內建的 `__builtin_clz` 是否會變成 `LZCNT`
```c=
#include <stdio.h>
int main()
{
unsigned int x=100;
int clz;
clz = __builtin_clz(x);
printf("x:%d clz(x):%d\n", x, clz);
return 0;
}
```
編譯過後對檔案反組譯:因為很長只列一小段
```
>0000000000400526 <main>:
400526: 55 push %rbp
400527: 48 89 e5 mov %rsp,%rbp
40052a: 48 83 ec 10 sub $0x10,%rsp
40052e: c7 45 f8 64 00 00 00 movl $0x64,-0x8(%rbp)
400535: 0f bd 45 f8 bsr -0x8(%rbp),%eax
400539: 83 f0 1f xor $0x1f,%eax
40053c: 89 45 fc mov %eax,-0x4(%rbp)
40053f: 8b 55 fc mov -0x4(%rbp),%edx
400542: 8b 45 f8 mov -0x8(%rbp),%eax
400545: 89 c6 mov %eax,%esi
400547: bf e4 05 40 00 mov $0x4005e4,%edi
40054c: b8 00 00 00 00 mov $0x0,%eax
400551: e8 aa fe ff ff callq 400400 <printf@plt>
400556: b8 00 00 00 00 mov $0x0,%eax
40055b: c9 leaveq
40055c: c3 retq
40055d: 0f 1f 00 nopl (%rax)
```
* 整份沒有任何關於`lzcnt`的指令， gcc 的 `__buildin_clz` 並沒有使用到 `lzcnt`
* 可以注意到` 400535: 0f bd 45 f8 bsr -0x8(%rbp),%eax`
* [bsr](http://x86.renejeschke.de/html/file_module_x86_id_20.html) 是 bit scan reverse 的縮寫，在 x86 指令集中代表找尋 most significant set bit (也就是第一個1)的位置，所以 clz 是根據 `bsr` 求得
#### 關於`lzcnt` 與 `bsr` 的差異
* `bsr`:
>==Searches the source operand (second operand) for the most significant set bit (1 bit). If a most significant 1 bit is found, its bit index is stored in the destination operand (first operand)==. The source operand can be a register or a
memory location; the destination operand is a register. The bit index is an unsigned offset from bit 0 of the source operand. ==If the content source operand is 0, the content of the destination operand is undefined==.
* `lzcnt`:
>==Counts the number of leading most significant zero bits in a source operand (second operand) returning the result
into a destination (first operand).==
LZCNT differs from BSR. For example, ==LZCNT will produce the operand size when the input operand is zero.== ==It should be noted that on processors that do not support LZCNT, the instruction byte encoding is executed as BSR.==
In 64-bit mode 64-bit operand size requires REX.W=1.
* `lzcnt` 的 flag affect:
>ZF flag is set to 1 in case of zero output (most significant bit of the source is set), and to 0 otherwise, CF flag is set to 1 if input was zero and cleared otherwise. OF, SF, PF and AF flags are undefined.
[SSE4的Wikipedia](https://en.wikipedia.org/wiki/SSE4) (註: `lzcnt` 雖然不屬於 SIMD 指令，卻是在與 SSE4 同時期引入的)
>The result of lzcnt is ==31 minus the result of the bsr (bit scan reverse), except when the input is 0==. ==lzcnt produces a result of 32, while bsr produces an undefined result==(and sets the zero flag). The encoding of lzcnt is similar enough to bsr that if lzcnt is performed on a CPU not supporting it such as Intel CPU's prior to Haswell, it will perform the bsr operation instead of raising an invalid instruction error.
另外在很多[硬體指令或軟體提供的函式](https://en.wikipedia.org/wiki/Find_first_set#Hardware_support)(~~其實幾乎都是~~ 更正:部份intel的處理器)，如果輸入為0都是未定義行為，不太懂硬體在發展上有什麼考量或困難，使得0沒辦法做處理呢?
>> 因為硬體不會知道你預期的資料表示法 (data representation)，像是 C 語言的 `uint32_t` 和 `uint64_t` 編譯為機械碼之後，就喪失了原始資訊，致使 `clz(0)` 的輸入無法確認是 16 個 0、32 個 0，還是 64 個 0 [name=jserv]
>>老師我原本看錯，以為`bsr`與`lzcnt`都不行處理 0，後來發現只有`bsr`是 undefined behavior ，`lzcnt`會回傳input size，`lzcnt`的狀況好像跟老師說的不太一樣嗎?倒是軟體提供的幾乎沒辦法提供正確結果，也是類似老師所說的嗎?[name=petermouse]
>> 不要隨便說「好像」，不要跟文組的人一樣「憑感覺」。這些指令都有明確定義，要先確定 gcc 的 `__builtin_clz` 到底用產生哪些指令 (去檢驗手冊和原始碼)，然後再去找 Intel 指令說明，像是 [LZCNT](http://www.felixcloutier.com/x86/LZCNT.html) 和 [Bit Manipulation Instruction Sets](https://en.wikipedia.org/wiki/Bit_Manipulation_Instruction_Sets) [name=jserv]
### 改成 __builtin_clz 的測試
在只有改`__builtin_clz`的狀況下，其實根本沒有改善的差別，因為他在natural-mergesort裡面佔的部份很少。
### 分析程式熱點
用`perf record ./natural-mergesort`
在用 `perf report` 觀察結果
```
Samples: 1K of event 'cycles:pp', Event count (approx.): 1488328678
Overhead Command Shared Object Symbol
31.95% natural-mergeso natural-mergesort [.] merge
14.49% natural-mergeso natural-mergesort [.] integer_cmp
13.41% natural-mergeso libc-2.23.so [.] _IO_vfscanf
6.45% natural-mergeso libc-2.23.so [.] vfprintf
2.66% natural-mergeso libc-2.23.so [.] 0x000000000014d45b
2.12% natural-mergeso libc-2.23.so [.] 0x000000000003b6da
2.03% natural-mergeso libc-2.23.so [.] __isoc99_fscanf
1.98% natural-mergeso libc-2.23.so [.] _IO_file_xsputn
1.72% natural-mergeso libc-2.23.so [.] 0x000000000003b6f0
1.47% natural-mergeso natural-mergesort [.] build_run_length_queue
1.24% natural-mergeso natural-mergesort [.] memcpy@plt
1.24% natural-mergeso natural-mergesort [.] main
1.08% natural-mergeso libc-2.23.so [.] 0x000000000014d2c3
1.08% natural-mergeso libc-2.23.so [.] 0x000000000014d45d
1.02% natural-mergeso libc-2.23.so [.] 0x000000000014d433
0.85% natural-mergeso libc-2.23.so [.] _IO_sputbackc
0.85% natural-mergeso libc-2.23.so [.] strchrnul
0.79% natural-mergeso natural-mergesort [.] run_length_queue_dequeue
0.74% natural-mergeso natural-mergesort [.] run_length_queue_enqueue
0.74% natural-mergeso libc-2.23.so [.] 0x000000000003b6d5
```
可以發現到 `merge` 是最花時間的地方，再來是 comparator
- [ ] merge()
**1. 連續小於的偵測**
在這個函式中，是逐個 element ，也就是以 integer 為單位比較，並且做 memcpy 。但是這會叫用很多次 memcpy ，所以可以試著讓 copy 的區塊大一點。
**2. Segment 分級**
畢竟 merge 要用到的比較很多，小到一定程度的 segment 就不必一定要用到 merge 來做，而可以用簡單的 qsort 來代替
>quick sort是不穩定排序 merge sort是穩定排序，如果需要穩定性替換就不能用quick sort了 [name=petermouse]
testcase : random 1000000 10000
```c
int conti_times_right = 0, conti_times_left = 0;
while (left + conti_times_left < left_bound && right + conti_times_right < right_bound) {
if (cmp(((char*) source) + size * (right + conti_times_right), ((char*) source) + size * (left + conti_times_left)) < 0) { //cmp() means? //A: user-defined comparasion function
if (conti_times_left) {
memcpy(((char*) target) + size * target_index,
((char*) source) + size * left,
size * conti_times_left);
left += conti_times_left;
target_index += conti_times_left;
conti_times_left = 0;
}
conti_times_right++;
} else {
if (conti_times_right) {
memcpy(((char*) target) + size * target_index,
((char*) source) + size * right,
size * conti_times_right);
right += conti_times_right;
target_index += conti_times_right;
conti_times_right = 0;
}
conti_times_left++;
}
}
if (conti_times_left) {
memcpy(((char*) target) + size * target_index,
((char*) source) + size * left,
size * conti_times_left);
left += conti_times_left;
target_index += conti_times_left;
}
else if (conti_times_right) {
memcpy(((char*) target) + size * target_index,
((char*) source) + size * right,
size*conti_times_right);
right += conti_times_right;
target_index += conti_times_right;
}
```
>> TODO: 參考 [Timsort](https://en.wikipedia.org/wiki/Timsort)，思考混合兩種以上 stable sorting algorithm 的機會。可透過 C macro 來建立 sorting template，參見 [sort](https://github.com/swenson/sort) [name=jserv]
---
### call graph
- ori
testcase - random - 1000000
testcase - random - 10000000
- opt
testcase - random - 10000000
**inline cmp**
因為cmp感覺是必要叫這麼多次的函式，目前也想不出其他可以減少的方法，所以選擇把 merge 部份的 cmp 做 inline
```c
inline __attribute__((always_inline)) int cmp(const void *a,const void *b)
{
return *(int *)a - *(int *)b;
}
```
### 建立測試資料
為了分析不同的測試資料會如何影響結果，所以建立了四種資料的形式
1. random : 資料完全隨機
2. sorted : 資料已經排序好(由小到大)
3. nearly-sorted : 資料呈現部份排序的樣子
4. reversed : 資料已經排序好(由大到小)
### 到目前為止所做的優化結果：
在不同種資料類型的mergesort時間,以及優化結果
資料為一百萬, range 為 0 到 50000000
以下為在資料類型是 random number 時在不同 range 下取1000個資料所得到的時間以及優化結果
在相同資料數量下, random number 在不同 range 下排列的時間是相近的,優化過後的結果更是趨於穩定
### branch predictor
記得在[Modern Microprocessors](http://www.lighterra.com/papers/modernmicroprocessors/) 提到一個小概念：
```clike
if (a > 7) {
b = c;
} else {
b = d;
}
```
* 原始組語
```
cmp a, 7 ; a > 7 ?
ble L1
mov c, b ; b = c
br L2
L1: mov d, b ; b = d
L2: ...
```
* Eliminating Branches with Predication
```
cmp a, 7 ; a > 7 ?
mov c, b ; b = c
cmovle d, b ; if le, then b = d
```
==不先判斷，而是先做如果式子成立的結果，再來判斷如果式子不成立的話要 override== ---> 試試看能不能把這種概念應用上去以減少判斷
且這個概念在[第18題 branch predictor](https://embedded2016.hackpad.com/2016q1-Week-2--Sc7AmIvN7EN#:h=%E7%AC%AC18%E9%A1%8C:-branch-predictor)中有類似的
```c
while (i1 < n && i2 < n) {
int v1 = src1[i1];
int v2 = src2[i2];
int take1 = v1 < v2;
dest[id++] = take1 ? v1 : v2;
i1 += take1;
i2 += (1 - take1);
}
```
但是要注意的是，merge 這個動作總共有三種寫法：origin, teacher version, student version 而詳細的探討要考慮到 ==gcc 的最佳化==，以及 CPU 架構（不確定能不能探討到）
>> 實際的結果呢？如何設計實驗？ [name=jserv]
* gcc的參數
`-fno-guess-branch-probability`
Do not guess branch probabilities using heuristics.
`-O0`表示無最佳化
`-O1`到`-O3`表示最佳化不同的程度
[參考資料](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)
## 兩份參考文件的閱讀
### 第一篇 :**uop cache**
## [Segmented Sort and Locality Sort](https://nvlabs.github.io/moderngpu/segsort.html)
我們看了幾次這份文件，本來以為 segmented sort 跟 locality sort 也是兩種 sort 的方法 ( 演算法或處理方法 ) ，但是發覺其實不是，這兩個 sort 是表現我們在 programming 時可能會面臨到的問題。
Segmented sort 與 locality sort : merge sort 的變形，使用在分配不平均的亂數資料上( non-uniform random data )
### 1. Segmented Sort
Segmented sort allows us to sort many variable-length arrays in parallel. ->將 sort 分成幾個 block 同時排列
圖片顯示如果平均 segment (資料長度)在較小的情況，平行處理時有較高的效能，反之如果 segment 太大，反而比一般 mergesort 效能低
* Segmented sort pairs
一組陣列的值對應到一個位置, sort 之後的陣列數值排列,因此也要將一開始的位置對應到 sort 後的值
>> 這是我們對 code 的理解,但仍不知道會應用在哪
>>
>>source code 的連結掉了，不知道裡面是否有應用的內容,而且github提供fork的repository找不到對應的檔案
>>[name=魏孜昀]
* Algorithm 解釋 : 文章後半段有討論到 segmented 想呈現的問題
一個單一的問題很好解決，但是當我們面臨到的是小區小區破碎的問題時，就要有不同手段的考量，就像當我們在 partition 不同的資料區間之後，不可能不論長度大小都用要叫多個 thread 來 sort ，因為 CP 值不符。
==未完==
### 2. Locality Sort
~~Locality sort : 可以偵測接近排序的資料範圍~~
這是關於 locality sort 的 benchmark，如果 noise amplitude 愈大，效能表現會遇差
noise amplitude 指的其實是資料上的浮動範圍，浮動範圍愈大會使得 locality sort 的效能降低
如果我們觀察 demo 的 input keys
```clike=
for(int i = 0; i < N; ++i)
keysHost[i] = i + Rand(0, 25);
```
呈現的資料:
```
Input keys:
0: 15 26 16 9 26 27 12 16 16 28
10: 13 32 36 18 30 40 28 35 34 44
20: 34 40 38 28 38 34 44 32 41 50
30: 55 55 37 52 36 57 38 48 39 47
40: 50 62 53 57 53 48 65 52 64 61
50: 70 61 76 72 79 64 60 77 61 84
60: 78 83 64 84 77 74 79 68 90 94
70: 82 92 82 95 91 76 95 77 91 94
80: 89 100 85 99 99 102 92 111 89 95
90: 109 114 98 96 105 103 113 119 107 105
```
這個 demo 程式讓資料呈現"有點排序"的樣子，noise amplitude 就是因為 `Rand(0,25)` 所造成資料上的浮動程度。
### 第二篇 :[Mergesort](https://nvlabs.github.io/moderngpu/mergesort.html)
1. Mergesort is type of ==multi-pass vectorized merge==: the first iteration executes N / 2 merges with inputs of length 1; the second iteration executes N / 4 merges with inputs of length 2; etc.
2. The number of batched merge passes is log(N) and the work per pass is N. This O(N log N) work-efficiency hurts mergesort's scalability compared to radix sort.
**number of batched merge passes :** 呼叫下一層 function
**work per pass :** 要處理的 list 的長度
所以 O(N log N) 比較 radix sort 的 O(variable x (N + 10)) ，當 N 變大，也就是 scalability 變大，mergesort 所受的影響會比較大
>> 上面整段數學式要說明什麼? 請重新整理式子並說明[color=red][name=課程助教]
>> 這個是當初在討論的計算式，忘記刪掉
>> [mgpuhost.cuh:163](https://github.com/ankan-ban/perft_gpu/blob/master/moderngpu-master/include/mgpuhost.cuh) 有提到 The data is sorted in-place. 顛覆我們一般對 mergesort 的認知，不知道是怎麼辦到的，有待確認 [name=petermouse]
>>> in-place : 不代表完全不使用額外空間，而是使用小數量或常數數量的額外空間 ( i.e. 空間複雜度 O(1))
3. This mergesort implementation runs at about half the throughput on large arrays with 32-bit keys as the fastest GPU radix sorts.
4. If the sort is truly on a critical path, it may be worth pulling a radix sort from B40C/FastSortSm20.
>> 不太懂第三句的 as the fastest GPU radix sorts 意思
>> 以及第四句的 a critical path 是什麼
>> **CTA-level blocksorts** ?? [name=heathcliffYang]
>>> [NVIDIA Compute](https://www.doc.ic.ac.uk/~wl/teachlocal/arch2/papers/nvidia-PTX_ISA_1.0.pdf) 中的 CTA 解釋 : A Cooperative Thread Array, or CTA, is an array of threads that execute a kernel concurrently or in parallel.
>>> [CUDPP library](http://cudpp.github.io/cudpp/2.1/dir_c96d66e6db050f7b70147ba9bf6583c2.html) 有一系列的 sort 函式[name=petermouse]
>>> [blocksort kernel](https://github.com/OrangeOwlSolutions/CUB/blob/master/BlockSort.cu)
>>
### Advantages of mergesort against radix sort
1. mergesort is comparison sort-> 支援多種 type 的 input ,而 radix sort 限制 input 在整數和小數
2. 數字大小不會影響 mergesort 時間複雜度,而數字大小會影響 radix 有幾個 bit 的多寡
3. 提供 CTA-level blocksort
4. 面對已經排序好的 input , mergesort 可以只是單純地將 data 複製，而其只需執行 log(N) ，radix sort 沒有檢查的機制
---
#### 補充 : radix sort
1. 將數值切割進行排序 ---> 92 變成 9 , 2
2. 以 base 來決定要分幾類 ---> 例如要排序的數值是10進位，那 base 就是10，依照單個位數大小分成0~9個類
3. 有幾個位數分幾類，從左到右與從右到左又是兩種方法
4. 最想要分析的 O(位數 x (數列長度 + base ))
[參考資料](http://notepad.yehyeh.net/Content/Algorithm/Sort/Radix/Radix.php)
---
5. Mergesort on GPU runs best when written in two distinct stages:
a) A blocksort kernel sorts random inputs into tile-length sorted lists, communicating with low-latency, high-bandwidth shared memory .The CTA blocksort forms a convenient re-usable component for MGPU's customers.
b) Multiple launches of MGPU Merge iteratively merge sorted lists, starting with the output of 1, communicating between passes with high-latency, high-capacity DRAM.
兩個都用到以下方法：
Pairs of threads (or CTAs for the global merge passes) cooperatively merge two VT-length lists (or two NV-length lists) into one list.
---> 覺得很神奇的地方是，居然是 " 一對 " 甚至更多個 thread 來一起處理 merge two list into one 的事
範例:
圖中每個正方形的上邊與右邊代表的是**已經排序好的資料**
對第一個正方形來說:
thread 1 與 thread 2 使用這兩個排序好的資料，排序到下圖中第一個正方形的**上邊**
對第二個正方形來說:
thread 3 與 thread 4 使用這兩個排序好的資料，排序到下圖中第一個正方形的**右邊**
接著，下圖 thread 1、2、3、4 會將第一個正方形的上邊與右邊，排序到更大一個正方形的上邊
thread 之間並行排序的細節:
在之前的圖中可以看到有條綠色的線從左上向下或向右走到右下，這條線稱作 Merge Path ，它由這兩筆排序好的資料的互相比對來呈現這兩筆資料的關係
接著對這張圖使用斜線切割成相等的大小，可以對照到x軸與y軸間，每個斜線之間的資料數目是一樣的
因此根據這個斜線切到的位置 diag 以及 Merge path，可以決定斜線之間的資料應該要放到哪個位置，比如說，從這張圖來看， a0 與 a1 之間的資料，以及 b0 與 b1 之間的資料，可以確定他們將會排在 diag0 與 diag1 之間，同理 a3 與 a4 之間的資料與 b3 與 b4 之間的資料將會排到 diag3 與 diag4 之間，有了區間的資訊，每個 thread 可以取得獨立的資料，排在各自的位置上，使得 thread 可以平行處理 merge
實際上在找尋斜線之間的區間(就是那個 a b 區間)時，並不是真的建立 Merge Path (不然直接 merge 就好啦!)，而是用 binary search 來找尋這個區間
>不太清楚 binary search 實際要怎麼應用在這裡 [name=petermouse]
---
### sorting networks
blocksort 在每個 CTA 中使用 NT 個 threads ，每個 thread 有 VT 個值，
在 merge 的前幾個回合可以使用 in-register 的 sort 方式 ， 分為
1. Batcher's odd-even mergesort (O(n log^2^ n))
2. odd-even transposition sort (O(n^2))
兩種在平行處理時有很好的效果，因為沒有使用到 shared memory
odd-even transposition sort 雖然比較慢，但是具備穩定性，所以推薦使用
### blocksort
CTAMergesort : 是一個 block-level mergesort
1. OddEvenTransposeSort : 先做簡單的 local sort
2. DeviceThreadToShared : 將 local sorted keys 放入 shared memory
3. DeviceBlocksortLoop : 將 list merge
---
---
6. Note that ==the keys and indices are unioned in shared memory, so as to not waste resources==. This is an idiom used throughput Modern GPU, and an important one to follow if you goal is high throughput.
---> 這個概念在下一章節會被用到
7. Segmented and locality sorts exploit the sortedness of inputs to reduce processes.
---> 比起 mergesort 的較為粗略的 partitioning ，Segmented sort 跟 locality sort 可以找到 " 已排序好的 " input 去減少執行過程
8. Mergesort is a multi-pass, out-of-place algorithm. The blocksort reduces global memory traffic by **sorting blocks of NV elements locally, performing key exchange through low-latency shared memory.**
---> 兩個方法減少資料溝通所需（因為 global 需要被所有 thread 看見)
>> 但是不是很了解實際執行的方法
>> 這個要順便複習其他組別研究 prefetch 和計算機組織結構中的 cacheline 對效能的影響 [name=jserv]
>> 跟 hugikun999 討論到可以試著用 AVX + 配上用 clz 當 key 來增加一次可以判斷的數值數量，但還沒開始實作
## out of our dictionary
> vanilla mergesort：normal mergesort
### 參考資料
* [intel 64 software developer instruction set reference manual ](http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-software-developer-instruction-set-reference-manual-325383.pdf)
* [gcc X86-Built-in-Functions](https://gcc.gnu.org/onlinedocs/gcc-4.8.4/gcc/X86-Built-in-Functions.html)
* [gcc implementation of builtin clz](http://stackoverflow.com/questions/9353973/implementation-of-builtin-clz)
* [gcc 查組語對照](http://leavinel.blogspot.tw/2012/06/gcc.html)
* [gcc 內嵌組語](http://ccckmit.wikidot.com/as:inlinec)
* [how undefined are builtin ctz0](http://stackoverflow.com/questions/19527897/how-undefined-are-builtin-ctz0-or-builtin-clz0)
* [GNU 的連結工具 (ld, nm, objdump, ar)](http://sp1.wikidot.com/gnulinker)
* [objdump](https://sourceware.org/binutils/docs/binutils/objdump.html)
* [NVlabs Cub](https://nvlabs.github.io/cub/)
* [moderngpu intro](https://nvlabs.github.io/moderngpu/intro.html)
* [optimizing merge sort](http://jeffreystedfast.blogspot.tw/2011/04/optimizing-merge-sort.html)
* [An Experimental Study of Sorting and Branch
Prediction](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9149&rep=rep1&type=pdf)
Sign in with Wallet
Connect another wallet