2016q3 Homework01 (compute-pi)

contrubuted by <hugikun999>
Reviewed by <heathcliffYang>

Reviewed by `heathcliffYang`

除了比較優化之餘，也可以嘗試其他更多pi的演算法去比較。

inplement

cycle finding
對於一個f(x)，給定一初始值，請找到第一次開使重覆的地方。（即Xi = Xj）一旦出現這種情況就會週期性出現。

memset(void *str,int c,size_t z)

說明：將這個指標所指向的前z個字節設為"c";

回傳值：pointer。

Memoization
一種最佳化的技巧，紀錄function call的結果，當有同樣的輸入發生時直接回傳cache的東西。
相關pi的演算法
時間函數
（1）UTC：世界標準時間。
（2）calender time:日曆時間，即相對時間，對不同的編譯系統會有不同標準時間點。
（3）clock tick:計時單元，其時間長短由CPU決定。
時間精準度：
（1）clock_gettime():取得目前時間。（1ns）
（2）Windows:QueryPerformanceCounter():取得rdtsc的counter數。( < 1ns)
(3)clock():精確度只到ms。
xtime：
（1）microsecond:微秒
（2）ms:毫秒
(3) 從1970/1/1開始計時，所謂的實時時鐘，精確度是us。
jiffies:
(1)每一個tick是以系統頻率決定，若是200MHz則可以精確到5ns。

AVX

延伸SSE指令集
這裡有歷代指令集的介紹。最大的特點是將XMM 128bit -> YMM 256bit，提高CPU在處理特定向量和浮點數運算的效率。
4個operands
通過VEX prefix實現4個operands，表示為dest、src1、src2、src3。雖然有4個operand，但速度不會到4x。VEX又可分為128和256版，分別集成escape prefix、SIMD prefix、REX prefix。
指令編碼格式

VEX的部份
（1）提供寄存器寻址。
（2）提供 escape prefix 与 SIMD prefix
（3）指示 operand 是 128 位还是 256 位
（4）提供 REX.RXB 的作用
實現不改變destination
AVX512
再將256bit擴大到512bit。
AVX相關語法
這裡是基本介紹每個字的意思。
這裡是AVX還有其它SSE的版本一些相關用法，裡面有每個function的做法和對應的assembly code。

compute_pi

在開使之前先用time測試了raytracing，結果發現由main.c的clock_gettime()會小於由time所測出來的real。

# Rendering scene
Done!
Execution time of raytracing() : 1.059423 sec

real	0m1.064s
user	0m8.364s
sys	0m0.028s

Makefile

(一)
我在寫Makefile的時候發現我的tab後面加的東西沒有變成指令，查了一下覺得可能是因為我好像把tab變成四個空白。
參考資料
（二）
在Makefile中的gencsv裡有寫到100、500、25000，經過測試知道100是N的初始值，500是每次增加的大小，25000是N的上限。

gencsv: default
    for i in `seq 100 500 25000`; do \
        printf "%d," $$i;\
        ./benchmark_clock_gettime $$i; \
    done > result_clock_gettime.csv

gnuplot

這次使用的是.csv檔和上一次不一樣，按照上次的去做runtime.gp會一直出現"runtime.gp", line 11: warning: Skipping data file with no valid points這個錯勿，參考Vivian Lin在其中多加一行set datafile separator ","就可以成功吃到資料了。

set xlabel 'N'
set ylabel 'Time(sec)'
set autoscale fix
set style fill solid
set datafile separator ","
set title 'runtime'
set term png enhanced font 'Verdana,10'
set output 'time.png'

plot 'result_clock_gettime.csv' using 1:2 title 'baseline' with linespoints, \
     'result_clock_gettime.csv' using 1:3 title 'openmp_2' with linespoints, \
     'result_clock_gettime.csv' using 1:4 title 'openmp_4' with linespoints, \
     'result_clock_gettime.csv' using 1:5 title 'AVX' with linespoints, \
     'result_clock_gettime.csv' using 1:6 title 'AVX+unroll' with linespoints, \

file

openMP
在編譯時記得要在gcc後加上參數-fopenmp，檔案內也要記得#include <omp.h>，不然就算加上一些openmp的語法也不會有作用。這裡有語法的介紹和一些範例。

#pragma omp parallel  num_threads(threads)

這個意思是指定threads，在（）內填入想指定的數目，通常thread越高就會跑越快，但是會有一個極限。

#pragma omp for private(x) reduction(+:pi)

for是將底下的for迴圈拆開來分給每個thread去執行。
private(x)是讓每個thread都有一份自己的x而不去影響到其它thread所存的東西，如果改到別的thread存的東西會造成錯誤。
reduction(+:pi)很類似private()，讓大家都有一份複本各自運算，不同的是最後再將所有的pi相加起來。

    register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;

宣告出4個256bit的register。__m256d是指4個64bit的double。

ymm0 = _mm256_set1_pd(1.0);

將4個double分別紀錄1.0進去。

ymm2 = _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0);

將4個double分別照順序記錄dt3、dt2 、dt*1、dt。

ymm4 = _mm256_setzero_pd();

全部bit記錄為0。

ymm3 = _mm256_add_pd(ymm3, ymm2);

將ymm3加上ymm2存到ymm3中。add可以替換sub、mul、div之類的算術方法，也可換為store之類的功能。

double tmp[4] __attribute__((aligned(32)));

(1)__attribute__可以設置函數屬性（Function Attribute）、變量屬性（Variable Attribute）和類型屬性（Type Attribute）。
(2)
aligned(32)
該聲明將強制編譯器確保（盡它所能)分配空間時採用32字節對齊方式。取32byte因為一個double是8byte，須要4個即為32byte。

未優化

(一)
(1)make check

time ./time_test_baseline
N = 400000000 , pi = 3.141593
3.49user 0.00system 0:03.49elapsed 100%CPU (0avgtext+0avgdata 1780maxresident)k
0inputs+0outputs (0major+85minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
3.58user 0.00system 0:01.79elapsed 199%CPU (0avgtext+0avgdata 1788maxresident)k
0inputs+0outputs (0major+86minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
3.80user 0.00system 0:00.95elapsed 398%CPU (0avgtext+0avgdata 1784maxresident)k
0inputs+0outputs (0major+92minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
1.01user 0.00system 0:01.01elapsed 100%CPU (0avgtext+0avgdata 1756maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
0.98user 0.00system 0:00.98elapsed 100%CPU (0avgtext+0avgdata 1792maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps

（2）由100～50000的圖

（二）
我更改了Makefile裡面的N，發現了N從0～20的區間都會跑很久，甚至是懷疑N=1就卡住了，從perf top 監看，發現一直compute_pi_avx處於熱點，目前還不確定是什麼原因，等看完AVX相關的文件應該就會比較了解。

優化

信賴區間
(一)
從圖中可以看到有些超出範圍太多，可以視為例外並沒有太多參考價值，故對這些數據取2個標準差可以使其有更精準的量測結果。
我在benchmark_clock_gettime.c中多寫了一個用來計算標準差的function，當中有多用到math.h這個library，所以在Makefile中必須多加上-lm才能執行。
（1）Makefile

default: computepi.o
    $(CC) $(CFLAGS) computepi.o time_test.c -DBASELINE -o time_test_baseline
    $(CC) $(CFLAGS) computepi.o time_test.c -DOPENMP_2 -o time_test_openmp_2
    $(CC) $(CFLAGS) computepi.o time_test.c -DOPENMP_4 -o time_test_openmp_4
    $(CC) $(CFLAGS) computepi.o time_test.c -DAVX -o time_test_avx
    $(CC) $(CFLAGS) computepi.o time_test.c -DAVXUNROLL -o time_test_avxunroll
    $(CC) $(CFLAGS) computepi.o benchmark_clock_gettime.c -o benchmark_clock_gettime -lm

注意-lm必須加在gcc的最後面，如果直接加在CFLAGS後面會錯。
（2）圖（sample = 10）

雖然次數的間隔比之前大，但是從圖中可以看到趨勢線已經比較少有極端值的存在了。

(3)圖（sample = 25）

這次取25個sample震盪幅度又變更小了，可是我的openmp4很奇怪，會有一些超乎常理的暴增。目前覺得是因為我的電腦只有2個threads，所以宣告成四個threads時會有互搶的問題，所以才有些會超高。另一個是有些openmp4比AVX還要低，AVX一次處理4個register，openmp4也像是一次跑4個的概念，但是卻比AVX還要快。

(4)圖（sample = 25）

驗證在（三）中的想法，我將AVX中加上了openmp去嘗試，結果發現圖非常的醜，應該是不能用private去讓每個thread有自己的256bit register，可是也有一些低的AVX+openmp可以低到跟openmp4一樣了，可見最低值還是有下降。
另外我增跑了openmp8，雖然沒有比openmp4更低卻有比openmp4更穩定。

Error check
（一）
除了花費時間外，從下面的圖中可以看出來由五種方法所算出來的正確性，剛開始N很小時所做出來的正確性很低，當N變大時正確性也提高。

參考自heathcliffYang

（二）
AVX和AVX+unroll會有震盪是因為每次加4和每次加16，如果N不是4或著16的倍數就會造成少算的結果，所以會產生震盪。將N和每次的差都改為16的倍數就不會和其它演算法有差。