2016q3 Homework 1 (compute-pi)

# 2016q3 Homework 1 (compute-pi) contributed by <`jeff60907`> ##### tags: `jeff60907` ### Reviewed by <`andy19950`> - 有詳細介紹 user/sys/real time 的差異，I/O-bound CPU-bond 也有用簡單的中文解釋讓我們容易了解，但希望能夠把這些知識加入到效能分析內，比如 clock_gettime() 傳不同的參數分別拿到的是什麼樣的時間，用哪種時間算會比較精準。 - 有充分了解各種方法的執行成果，也可以嘗試看看用這些方法來證明看看別的演算法，像是萊布尼茲或蒙地卡羅之類的。 - 可以在 Makefile 修改迭代的次數，讓 gnuplot 畫出來的結果看起來比較好看，如果能取 95% 信賴區間更好。 - 最後是下面 compute_ci 的程式碼，在 hackMD 上看格式好像是錯的，我有嘗試修改但好像沒用，可能是在複製貼上的時候有貼到特殊字元吧。 --- ## 開發環境 ``` 作業系統 Lubuntu 16.04 CPU:Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz L1d cache: 32K L1i cache: 32K L2 cache: 256K L3 cache: 8192K ``` * 複習 [微積分概論](https://www.junyiacademy.org/junyi-pre-univ/root/junyi-math/mu-cal-ch1) * 閱讀圓周率π [學長詳盡分析共筆](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I) * [閱讀Function](http://www.csie.ntnu.edu.tw/~u91029/Function2.html) 其他演算法近期全部看一輪，惡補自己的未修過演算法的不足 * [2016q3 Homework 1 (compute-pi)](https://hackmd.io/JwdgZmCGBsDGCmBaAzPWAGRAWWtqOABNkAORESAVnUqzACMxkx0g?both) --- ## 進行實驗之前 ### 時間處理與 time 函式使用 [What do 'real', 'user' and 'sys' mean ](http://stackoverflow.com/questions/556405/what-do-real-user-and-sys-mean-in-the-output-of-time1) [linux time 的指令](http://yuanfarn.blogspot.tw/2012/08/linux-time.html) * Real is wall clock time time from start to finish of the call. This is all elapsed time including time slices used by other processes and time the process spends blocked (for example if it is waiting for I/O to complete). > 從一開始到結束所花費的時間，包含了 processes跟blocked 及 I/O的等待 * User is the amount of CPU time spent in user-mode code (outside the kernel) within the process. This is only actual CPU time used in executing the process. Other processes and time the process spends blocked do not count towards this figure. > 在user mode code 執行的CPU時間，只有CPU在程序的時間，其他processes跟blocked不計算再內。 * Sys is the amount of CPU time spent in the kernel within the process. This means executing CPU time spent in system calls within the kernel, as opposed to library code, which is still running in user-space. Like 'user', this is only CPU time used by the process. See below for a brief description of kernel mode (also known as 'supervisor' mode) and the system call mechanism. > process 在kernel進行的CPU時間，像是 systme call * real time = sys + user ? >1.下指令 sleep（） real time > sys + user ，因為sleep 指令在 user mode 及 kernel mode 都沒有佔用 CPU >2.使用多核心計算可能會造成 real time < sys + user ，uesr是各個核心計算的總合 * 時間單位 * 秒用s表現，毫秒用ms,微秒用μs表示，納秒用ns表示，皮秒用ps表示 * **1\. Wall-clock time** 顧名思義就是掛鐘時間，也就是現實世界中實際經過的時間，是由 kernel 裡的 xtime 來紀錄，系統每次啟動時會先從設備上的 [RTC](https://zh.wikipedia.org/wiki/%E5%AF%A6%E6%99%82%E6%99%82%E9%90%98) 上讀入 xtime。 * **2\. CPU time** 指的是程序在 CPU 上面運行消耗 (佔用) 的時間，如果多核心運算 threads，CPU time 是「每條 thread 的使用時間加總」，單純計算 CPU time 是不夠的，因為執行時間可能還包括 I/O time、 communication channel delay、synchronization overhead....等等。 [CPU bound](https://en.wikipedia.org/wiki/CPU-bound) 與 [I/O bound](https://en.wikipedia.org/wiki/I/O_bound) 所謂 CPU bound，指的是程序需要大量的 CPU computation (對 CPU time 需求量大)，相比之下僅需少量的 I/O operation，效能取決於 CPU 速度。所謂 I/O bound，需要大量 I/O operation (I/O request 的頻率很高或一次 I/O request 成本很高)，但僅需少量的 CPU computation，效能取決於 I/O device 速度。很多時候我們可以使用 time 判斷程式是 CPU bound 還是 I/O bound。 1. `real < user`: 表示程式是 CPU bound，使用平行處理有得到好處，效能有提昇。 2. `real ≒ user`: 表示程式是 CPU bound，即使使用平行處理，效能也沒有明顯提昇，常見的原因有可能是其實計算量沒有很大，生成、處理、結束多執行緒的 overhead 吃掉所有平行處理的好處，或是程式相依性太高，不適合拿來作平行處理。 3. `real > user`: 表示程式是 I/O bound，成本大部分為 I/O操作，使得平行處理無法帶來顯著的效能提昇 ( 或沒有提昇)。 ## 嘗試實驗 ### 探討Makefile 及各檔案 [閱讀Louei Lu](https://hackmd.io/EYTg7MAmBmULQDZgGMCGcAsAmADAUzgA5IBWA1EsDZBEcZQvIA==#) ``` check: default time ./time_test_baseline time ./time_test_openmp_2 time ./time_test_openmp_4 time ./time_test_avx time ./time_test_avxunroll ``` > $ make check 會執行default在各自執行這些指令 > 如果想單獨看某一檔案執行時間 > $ time ./time_test_baseline ``` gencsv: default for i in `seq 100 5000 25000`; do \ printf "%d," $$i;\ ./benchmark_clock_gettime $$i; \ done > result_clock_gettime.csv ``` > 可以看成for 初始100 每次 +5000 到 25000 > 然後把資料放進 .csv excel可讀檔 * 增加 make plot 到Makefile，先執行gencsv 在運作gnuplot ``` plot: gencsv gnuplot runtime.gp ``` > 要先把runtime.gp加進資料夾內，可參考phonebook ### 測試效能 ``` $ make check ``` ``` time ./time_test_baseline N = 400000000 , pi = 3.141593 2.93user 0.00system 0:02.93elapsed 99%CPU (0avgtext+0avgdata 1788maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_openmp_2 N = 400000000 , pi = 3.141593 2.93user 0.00system 0:01.47elapsed 198%CPU (0avgtext+0avgdata 1724maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_openmp_4 N = 400000000 , pi = 3.141593 3.71user 0.00system 0:01.10elapsed 335%CPU (0avgtext+0avgdata 1604maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_avx N = 400000000 , pi = 3.141593 1.34user 0.00system 0:01.34elapsed 99%CPU (0avgtext+0avgdata 1696maxresident)k 0inputs+0outputs (0major+84minor)pagefaults 0swaps time ./time_test_avxunroll N = 400000000 , pi = 3.141593 1.18user 0.00system 0:01.18elapsed 99%CPU (0avgtext+0avgdata 1796maxresident)k 0inputs+0outputs (0major+84minor)pagefaults 0swaps ``` ### 各自探討效能 [TempoJiJi hackmd](https://hackmd.io/KwBgpgxghhAcsFpYgGyICwEYXqZgJgEYL4DMmYAnAEz4Bmm41QA=?view) ### baseline 透過積分式子的實作 ```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; } ``` ``` $ time ./time_test_baseline ``` ``` N = 400000000 , pi = 3.141593 real 0m3.011s user 0m3.008s sys 0m0.000s ``` > real ≒ user 表示程式是 CPU bound，即使使用平行處理，效能也不會明顯提昇。 ### openmp2 ``` $ time ./time_test_openmp(2/4) ``` ```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; } ``` ``` N = 400000000 , pi = 3.141593 real 0m1.688s user 0m3.368s sys 0m0.000s ``` ### openmp4 ``` N = 400000000 , pi = 3.141593 real 0m1.562s user 0m4.748s sys 0m0.000s ``` > real < user: 表示程式是 CPU bound，使用平行處理有得到好處，效能有提昇。 ### avx AVX 是一組 SIMD(Single Instruction Multiple Data)的指令集，擴展了原有的 MMX 及 Intel Streaming SIMD Extensions(Intel SSE)，增加了下列的 features * 原來 128-bit 的 SIMD 暫存器擴展至 256-bit，未來可能會支援 512-bit 甚至 1024-bit * Three-operand，原來的 two-operand 運算是 A=A+B，將結果寫到來源暫存器，而新的是 A=B+C，舊的仍然存在 * 有一些指令可以支援四個暫存器的運算元，使得 code 更小，並且減少一些不必要的指令，提升速度 * 有一個新的 extension coding scheme(VEX) ``` $ time ./time_test_avx ``` ``` N = 400000000 , pi = 3.141593 real 0m1.359s user 0m1.360s sys 0m0.000s ``` > real ≒ user 表示程式是 CPU bound，即使使用平行處理，效能也不會明顯提昇。 ### avxunroll ``` $ time ./time_test_avxunroll ``` ``` N = 400000000 , pi = 3.141593 real 0m1.246s user 0m1.244s sys 0m0.000s ``` > real ≒ user 表示程式是 CPU bound，即使使用平行處理，效能也不會明顯提昇。 * 設定gnuplot格式進行圖表 ``` reset set key width 2 set key left set style fill solid set yrange[0:] set xrange[0:61000] set grid set title 'perfomance comparison' set ylabel 'time(sec)' set xlabel 'N' set datafile separator "," set term png enhanced font 'Verdana,10' set output 'runtime.png' plot "result_clock_gettime.csv" using 1:2 with lines title 'baseline', \ '' using 1:3 with lines title 'omp_2', \ '' using 1:4 with lines title 'omp_4', \ '' using 1:5 with lines title 'avx', \ '' using 1:6 with lines title 'avxunroll' ``` > runtime .gp 格式 xyrange可自行設定 > 圖一直編譯失敗，才發現CSV檔是以","做格子區分 > set datafile separator "," [參考Vivan Lin ](https://hackmd.io/MYZgpgTAJmCsUFpYDNiwQFnMhBOAhmAGx5gYCMY+IUuRy5QA?view#) > clokc_getime loop = 25 是25次取樣後才計算時間 ![](https://i.imgur.com/Itr02ei.png) > 之後將取樣loop 改為 1 ![](https://i.imgur.com/rTQ2s51.png) > 這張圖發現 omp4 一開始抖動幅度很誇張，但是後面幾乎趨於平穩不在抖動 ![](https://i.imgur.com/ODlgsh5.png) > 把疊帶次數縮短來看 omp4在500 >發現把疊代數變多才看的出抖動的情形 >這幾張baseline抖動的效果最大，omp4執行效能卻是最高的但是抖動程度也是很高 >為什麼會造成這樣的抖動呢？ omp4其他同學效果都不是最好的 why? ### 實作學長分析共筆的error rate * 先探討benchmark_clock_gettime.c內使用的 clock_gettime函式的原型為： ``` int clock_gettime(clockid_t clk_id, struct timespec *tp); ``` 應用在計時的時候，第一個參數clk_id可填入CLOCK_REALTIME，而struct timespec* tp則是本函式回傳的結果。struct timespec的宣告如下： ``` struct timespec { time_t tv_sec; /* seconds */ long tv_nsec; /* nanoseconds */ }; ``` 使用這種時間函數的話，可能可以精確到奈秒等級。如果想要進一步確定可容許的精確度，則可以使用clock_getres，其原型如下： ``` int clock_getres(clockid_t clk_id, struct timespec *res); ``` * time() returns the wall-clock time from the OS, with precision in seconds. * clock() seems to return the sum of user and system time. * clock_gettime(CLOCK_MONOTONIC, ...) provides nanosecond resolution. * gettimeofday() returns the wall-clock time with (nominally) µs precision [參考C語言：使用clock_gettime計算程式碼的時間需求](https://starforcefield.wordpress.com/2012/08/06/c語言：使用clock_gettime計算程式碼的時間需求/) [Measure time in Linux - time vs clock vs getrusage vs clock_gettime vs gettimeofday vs timespec_get?](http://stackoverflow.com/questions/12392278/measure-time-in-linux-time-vs-clock-vs-getrusage-vs-clock-gettime-vs-gettimeof) 參考學長共筆的程式碼 ```c= double compute_error(double Pi){ Pi=fabs(Pi-M_PI)/M_PI; return Pi; } ``` #### error rate 分析 ![](https://i.imgur.com/ODUCUpD.png) > 看到這張圖的時候以為自己是不是做錯了，怎麼avxunroll呈現鋸齒狀 > 然後9/29早上老師PO了一篇文原來也有同學發現N在1000結尾的時候會發生問題 ![](https://i.imgur.com/MycL0Kp.png) > 把間距拉小看，可以發現AVX再約疊代一千次後幅度空間會變小 > avxunroll 可以清楚的看到每約1000會有一個折反點，是什麼原因造成呢？ ![](https://i.imgur.com/ALQGdLW.png) > [參考王紹華同學](https://hackmd.io/KYQwjARgZgTArAZgLQgJwAYpICyoMYAmSAHJiEggGzzqQEIDsqYQA===#) 他提到avxunroll是因為avxunroll的 loop code執行一次等於 baseline 執行16次，所以會造成誤差，AVX 的版本是因為每個 for loop 只等同於 baseline 版本執行了4次，剛好可以整除資料點，才沒有算錯 > ex：N=20，上述迴圈只會執行一次，而導致資料少了20 - 16 = 4筆 ### 實作看看N疊代差異 > seq 100 4800 500000 ![](https://i.imgur.com/q9BaIzp.png) > 當設定間距為4800 / 16 的整數倍可以讓 avxunroll 不會因為誤差抖動，產生平滑曲線 > seq 112 4800 500000 ![](https://i.imgur.com/RcFPSC8.png) > 當起使點設定 112 /16 的整數倍卻讓avx 誤差完全消失，疑惑為什麼當起始點為16倍數反而造成這個現象？ ### 信賴區間(Confidence Interval) #### Steps: * 計算平均值 * 計算標準差(standard deviation) * 平均值的正負標準差乘以 2 就是 95% 信賴區間 #### 標準差算法 ![](https://i.imgur.com/ecyC4vF.png) #### 加上信賴區間實作測 > 取樣 loop = 25 ![](https://i.imgur.com/WxtI5Ua.png) > 發現抖動情形沒有改變很多，在想程式碼是不是打錯 > [參考王紹華同學](https://hackmd.io/KYQwjARgZgTArAZgLQgJwAYpICyoMYAmSAHJiEggGzzqQEIDsqYQA===#) 程式碼做修改 > 把取樣數 25個放進來計算 95% ```c= double compute_ci(double data[]) { double mean = 0; double SD = 0; // Standard Deviation double SDerror = 0; double low,top; double result = 0; int count = 0; //compute mean for(int i = 0; i <loop; i++){ mean += data[i]; } mean = mean / loop; //compute Standard Deviation for(int i = 0; i < loop; i++) { SD += pow((data[i] - mean),2); } //Calculate standard error SDerror = sqrt(SD / (double)loop); //standard range low = mean - 2 * SDerror; top = mean + 2 * SDerror; // in 95% for(int i = 0; i < loop; i++) { if(data[i] >= low && data[i] <= top) { count++; result+= data[i]; } } result = result / count; return result; } ``` ### 驗證 Leibniz formula for π ![](https://i.imgur.com/LUKVYwO.png) ```c= double compute_pi_leibniz(size_t N) { double pi = 0.0; for (size_t i = 0; i < N; i++) { pi += pow(-1.0,i) / (2*i + 1); } return pi * 4.0; } ``` > 測試了加上 -lm 讓math.h運作也無法，原來也有同學遇到這個問題 > [參考Vivan Lin](https://hackmd.io/MYZgpgTAJmCsUFpYDNiwQFnMhBOAhmAGx5gYCMY+IUuRy5QA?view#) > 把平方拆成直接用 ![](https://i.imgur.com/wPmNblD.png) > 看起來leibniz 抖動效果看起來不怎麼明顯 * 未完其他之後補上

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.