2016q3 Homework1(compute_pi)

# 2016q3 Homework1(compute_pi) ## 目標 * 學習透過離散微積分求圓周率 * Leibniz formula for π * 積分計算圓周率π * Function * 著手透過 SIMD 指令作效能最佳化 * [作業說明](https://hackmd.io/JwdgZmCGBsDGCmBaAzPWAGRAWWtqOABNkAORESAVnUqzACMxkx0g#a03-compute-pi) ## Makefile 每寫一份功課又是學習Makefile的時候=w= * 編譯時有`-DBASELINE`等等，其中的`-D`指的是`#define BASELINE 1` * 這樣程式碼就可以塞在同一個.c .h中！只要在編譯的時候加上`-D`就好了！！ * Makefile的篇幅可以大幅降低！ * 我想phonebook也可以這樣改寫，不然不同的方法就還要在多生一組.c .h ## 實驗 * `make check`結果： ```shell== time ./time_test_baseline N = 400000000 , pi = 3.141593 7.59user 0.00system 0:07.59elapsed 100%CPU (0avgtext+0avgdata 1560maxresident)k 8inputs+0outputs (0major+82minor)pagefaults 0swaps time ./time_test_openmp_2 N = 400000000 , pi = 3.141593 8.03user 0.00system 0:04.01elapsed 200%CPU (0avgtext+0avgdata 1572maxresident)k 0inputs+0outputs (0major+84minor)pagefaults 0swaps time ./time_test_openmp_4 N = 400000000 , pi = 3.141593 8.32user 0.00system 0:02.16elapsed 384%CPU (0avgtext+0avgdata 1676maxresident)k 0inputs+0outputs (0major+89minor)pagefaults 0swaps time ./time_test_avx N = 400000000 , pi = 3.141593 3.30user 0.00system 0:03.30elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k 0inputs+0outputs (0major+83minor)pagefaults 0swaps time ./time_test_avxunroll N = 400000000 , pi = 3.141593 1.73user 0.00system 0:01.73elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k 0inputs+0outputs (0major+82minor)pagefaults 0swaps ``` * `make gencsv` ==> 輸出至表單輸出結果與在Makefile中看到的程式碼只差在`$$i` 以下為輸出： ```shell== for i in `seq 100 5000 25000`; do \ printf "%d," $i;\ ./benchmark_clock_gettime $i; \ done > result_clock_gettime.csv ``` 原來除了寫進txt檔之外，還可以寫進libre office的表單裡，好酷！ * 理解程式碼： * 上面的程式碼相當於： ```C== for(i=100; i<25000; i+=5000) { printf("%d",i); benchmake_clock_gettime(i); } ``` 然後把結果輸出至.csv檔 :::warning 細節的語法待會查先highlight ::: * 在Makefile中，在執行各檔時前面多加了`time` * time: 查看指令的執行時間 * 直接執行`time ./time_test_baseline` ``` N = 400000000 , pi = 3.141593 real 0m7.805s user 0m7.808s sys 0m0.000s ``` * 顯示了三種時間：real,user和sys Real refers to actual elapsed time; User and Sys refer to CPU time used only by the process. * **real**: real time,a wall clock time,指令從開始到結束執行的時間，含其他process * **user**: user CPU time,指令在user mode的時間,其他process的時間不算 * **sys**: system CPU time,指令在kernel mode的時間,其他process的時間不算 * `real = user + sys` * ==**並非如此**==，如sleep(),不論在user mode或kernel mode沒佔用CPU時間，那就不會被算入，反觀real time仍會把sleep()的時間算進去 ### Baseline ![](https://i.imgur.com/H0I0Cig.png) 為上述積分式的實作，最後回傳的值要*4，因為算出來是 $π \over4$ * N愈大[0,1]就被切割成更多等份，dt也就愈小 ```C== double compute_pi_baseline(size_t N) { double pi = 0.0; double dt = 1.0 / N; // dt = (b-a)/N, b = 1, a = 0 for (size_t i = 0; i < N; i++) { double x = (double) i / N; // x = ti = a+(b-a)*i/N = i/N pi += dt / (1.0 + x * x); // integrate 1/(1+x^2), i = 0....N } return pi * 4.0; } ``` ### OpenMP * 也是用跟baseline一樣的積分式求得π * 可以傳遞參數決定要開幾個threads * 在`time_test.c`中分別傳入2個和4個threads的實驗 * 程式碼： * private(x)：讓每個thread都有自己的x,不然threads之間會把x當成共用的變數 * reduction(+:pi)：讓各個thread針對指定的變數擁有一份有起始值的複本 * 支援的運算元有 +, *, –, &, ^, |, &&, ||；而變數則必須要是 shared 的 ```C== double compute_pi_openmp(size_t N, int threads) { double pi = 0.0; double dt = 1.0 / N; double x; #pragma omp parallel num_threads(threads) { #pragma omp for private(x) reduction(+:pi) for (size_t i = 0; i < N; i++) { x = (double) i / N; pi += dt / (1.0 + x * x); } } return pi * 4.0; } ``` ### AVX SIMD * 甚麼是AVX？ * AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集 * AVX指令集借鑑了一些AMD SSE5的設計思路，進行擴展和加強，形成一套新一代的完整SIMD指令集規範。 * SSE將向量處理能力從64位提升到128位後，SSE系列指令都只能使用128位XMM暫存器，這次AVX將所有16個128位XMM暫存器擴充為256位的YMM暫存器，從而支持256位的向量計算。 * prototype: * _mm256 : single precision,32-bit * _mm256d : double precision, 64bit 配上註解與下方的連結終於稍微看懂了~ * `__attribute__ `機制，可以用來設置函數屬性（Function Attribute）、變數屬性（Variable Attribute）和類型屬性（Type Attribute）。 * `aligned` 則規定變數或結構的最小對齊格式，以 Byte 為單位。 >不知道為甚麼沒有highlight@@ ```C== double compute_pi_avx(size_t N) { double pi = 0.0; double dt = 1.0 / N; register __m256d ymm0, ymm1, ymm2, ymm3, ymm4; ymm0 = _mm256_set1_pd(1.0); ymm1 = _mm256_set1_pd(dt); ymm2 = _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0); ymm4 = _mm256_setzero_pd(); // sum of pi for (int i = 0; i <= N - 4; i += 4) { ymm3 = _mm256_set1_pd(i * dt); // i*dt, i*dt, i*dt, i*dt ymm3 = _mm256_add_pd(ymm3, ymm2); // x = i*dt+3*dt, i*dt+2*dt, i*dt+dt, i*dt+0.0 ymm3 = _mm256_mul_pd(ymm3, ymm3); // x^2 = (i*dt+3*dt)^2, (i*dt+2*dt)^2, ... ymm3 = _mm256_add_pd(ymm0, ymm3); // 1+x^2 = 1+(i*dt+3*dt)^2, 1+(i*dt+2*dt)^2, ... ymm3 = _mm256_div_pd(ymm1, ymm3); // dt/(1+x^2) ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += dt/(1+x^2) } double tmp[4] __attribute__((aligned(32))); _mm256_store_pd(tmp, ymm4); // move packed float64 values to 256-bit aligned memory location pi += tmp[0] + tmp[1] + tmp[2] + tmp[3]; return pi * 4.0; } ``` * 程式碼解讀(可能有誤！！) * `__m256d`:存放double,256-bit可用，因此可以存放256/64 = 4 個double * 每次的運算是4個double **個別** 加總 * 最後使用tmp[4]把4個值存起來，最後在加總 * `(aligned(32))` 是因為double有8個byte,又宣告 double tmp[4] ==> 8*4= 32(byte) ### Leibniz ![](https://i.imgur.com/gdKlini.png) * 程式： ```C== double compute_pi_Leibniz(size_t N) { double pi = 0.0; double x; for (size_t i = 0; i < N; i++) { x = pow(-1,i) / (2*i +1); pi += x; } return pi * 4.0; } ``` 一直出現錯誤資訊：應該是link的問題，但是即使加了`-lm`還是沒用@@ :::warning computepi.c:(.text+0xd80): 未定義參考到 pow collect2: error: ld returned 1 exit status make: *** [default] Error 1 ::: > computepi.c 有included <math.h>嗎？或許可以試試，我花很多時間處理這個...[name=SarahYuHanCheng] * 只好先改成 ```C== int tmp = (i%2) ? -1 : 1; x = (double)tmp / (2*i +1) ; //記得轉型 pi += x; ``` ### Leibniz AVX 看了一些資料在搭配原有的程式碼，很快就可以寫出來了~ 覺得開心=w= ```C== double compute_pi_leibizavx(size_t N) { double pi = 0.0; register __m256d ymm0, ymm1, ymm2, ymm3, ymm4; ymm0 = _mm256_set_pd(1.0,-1.0,1.0,-1.0); //for pow(-1,i) ymm1 = _mm256_set1_pd(2.0); ymm2 = _mm256_set1_pd(1.0); ymm4 = _mm256_setzero_pd(); // sum of pi for (int i = 0; i <= N - 4; i += 4) { ymm3 = _mm256_set_pd(i, i+1.0, i+2.0, i+3.0); //i i+1 i+2 i+3 ymm3 = _mm256_mul_pd(ymm1, ymm3); //2*i 2*(i+1)... ymm3 = _mm256_add_pd(ymm3,ymm2); //2*i+1 2*(i+1)+1 ... ymm3 = _mm256_div_pd(ymm0,ymm3); //(-1)^i/(2*i+1) ... ymm4 = _mm256_add_pd(ymm3,ymm4); //sum } double tmp[4] __attribute__((aligned(32))); _mm256_store_pd(tmp, ymm4); // move packed float64 values to 256-bit aligned memory location pi += tmp[0] + tmp[1] + tmp[2] + tmp[3]; return pi * 4.0; } ``` ### 其他計算圓周率的方法 * [蒙地卡羅法](https://zh.wikipedia.org/wiki/%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85%E6%96%B9%E6%B3%95)： * 邊長為2的正方形，正方形內有一個內切圓，面積為 $1 \cdot 1 \cdotπ = π$ * 圓面積和正方形面積之比為π:4 = 正方形內均勻隨機丟石頭，落在圓形內的機率 ```C== double compute_pi_MonteCarlo(size_t N) { int count = 0; srand(time(NULL)); for (size_t i = 0; i < N; i++) { double x = (double)rand()/RAND_MAX; double y = (double)rand()/RAND_MAX; if((x*x+y*y)<=1) count++; } return 4.0 * count / N ; } ``` ### Plot * 複習語法 * `using [a：b]` : 以第a column當x座標,第b column當y座標製圖時產生的warning: :::warning warning: Skipping data file with no valid points ::: 查了一段時間，玄機就在副檔名`.csv`和gnuplot他吃數據的方式啊！！ **csv : comma separated value** 所以資料之間是以 ==**逗號**== 區隔，因此要在gnuplot script裡加一行`set datafile separator ","`才讀的到資料 * 不同方法比較圖 :::warning 抖的超厲害！ ::: 取樣從N=1000~N=250000，間隔1000 ![](https://i.imgur.com/A5GYP7c.png) 隔一天跑出來的結果差好多@@ omp_4沒那麼劇烈了，還不知道為甚麼會這樣 ![](https://i.imgur.com/GimWEAW.png) **新增Leibniz** * 可以發現Leibniz與omp_2的時間差不多，但又比omp_2還穩定 * LeibnizAVX也比avx的時間再更短 * LeibnizAVX_unroll語omp_4幾乎重合 ![](https://i.imgur.com/BI10Fxo.png) ### 思考如何精準計算時間 * walk clock time: 由xtime記錄，系統每次啟動時會先從設備上的 [RTC](https://zh.wikipedia.org/wiki/%E5%AF%A6%E6%99%82%E6%99%82%E9%90%98) 上讀入 xtime * xtime:從cmos電路中取得的時間，一般是從某一歷史時刻開始到現在的時間 * Jeffie:紀錄系統開幾以來，經過多少個 tick，取決於系統的頻率，單位是Hz * system time: * 系統上的時間，開機時它會讀取RTC 來設定，可能過程中還會有時區換算等之類的設置。 ### 信賴區間 ## 參考資料 [TempoJiJi hackmd](https://hackmd.io/KwBgpgxghhAcsFpYgGyICwEYXqZgJgEYL4DMmYAnAEz4Bmm41QA=?view) * Time * [time_ref1](http://yuanfarn.blogspot.tw/2012/08/linux-time.html) * [time_ref2](http://stackoverflow.com/questions/556405/what-do-real-user-and-sys-mean-in-the-output-of-time1) * [容易混淆LINUX時鐘的xtime和jiffies](http://www.coctec.com/docs/linux/show-post-186188.html) * OpenMP * [openmp_private](http://viml.nchc.org.tw/blog/paper_info.php?CLASS_ID=1&SUB_ID=1&PAPER_ID=17) * [openmp_reduction](https://kheresy.wordpress.com/2006/09/22/%E7%B0%A1%E6%98%93%E7%9A%84%E7%A8%8B%E5%BC%8F%E5%B9%B3%E8%A1%8C%E5%8C%96%EF%BC%8Dopenmp%EF%BC%88%E4%BA%94%EF%BC%89-%E8%AE%8A%E6%95%B8%E7%9A%84%E5%B9%B3%E8%A1%8C%E5%8C%96/) * AVX * [introduction](http://www.expreview.com/13236.html) * [manual](http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-optimization-manual.pdf)(chap11) * [functions](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=579,2050,100,4629,4578,112,100&techs=AVX) ###### tags: `system embedded` `HW1-3`

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.