[compute-pi](https://hackmd.io/s/rJARexQT)

# [compute-pi](https://hackmd.io/s/rJARexQT) [github](https://github.com/diana0651/compute-pi) contributed by <`Diana Ho`> ###### tags: `d0651` `sys` ## 案例分析 ### 學習目標 - [ ]學習透過離散微積分求圓周率 * [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80) * [積分計算圓周率π](http://book.51cto.com/art/201506/480985.htm) * [Function](http://www.csie.ntnu.edu.tw/~u91029/Function2.html) - [ ]著手透過 SIMD 指令作效能最佳化 ### 開發環境 Ubuntu 14.04 LTS --- ## Wall-clock time vs. CPU time * clock(), clock_gettime() 的使用 * clock() returns an approximation of processor time used by the program. ```clike= clock_t begin = clock(); compute_pi; clock_t end = clock(); double time = (double)(end - begin) / CLOCKS_PER_SEC; ``` * clock_gettime() retrieves the time of the specified clock clk_id. ```clike= struct timespec { time_t tv_sec; // seconds long tv_nsec; //nanoseconds } struct timespec start,end; clock_gettime(CLOCK_MONOTONIC_RAW,&start); compute_pi; clock_gettime(CLOCK_MONOTONIC_RAW,&end); double time = (double)((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/1000000000); ``` * time(), gettimeofday() 的使用 * time() returns the wall-clock time from the OS, with precision in seconds. * gettimeofday() returns the wall-clock time with nominally µs precision. * 為什麼 clock_gettime() 結果飄忽不定？ * function的執行時間單位以nanoseconds來決定，執行時間很短時並無法排除Linux背景處理其他事情，導致額外的cache miss，所以需要很多針對某一筆參數作累積再平均，才能避免飄忽不定 * 為什麼 time() 和 gettimeofday() 不適合拿來作 benchmark ? * clock_gettime()是monotonic，而time() 和 gettimeofday()的精準度較低，所以不適合 >>[參考概念](https://embedded2015.hackpad.com/HW1-Mj4fZIhTKVH) --- ## 效能分析比較圖表預期執行 `$ make plot` 後，可透過 gnuplot 產生效能分析比較圖表 - 舊電腦: 每當執行到 avxunroll 時, 就會發生不合法命令的錯誤 ```clike= $ make check make: *** [check] Error 132 $ make gencsv Illegal instruction (core dumped) Illegal instruction (core dumped) Illegal instruction (core dumped) Illegal instruction (core dumped) Illegal instruction (core dumped) make: *** [gencsv] Error 132 $ time ./time_test_avxunroll 不合法的命令 (core dumped) real 0m1.219s user 0m0.000s sys 0m0.000s ``` :::warning > [參考](https://hackmd.io/KwUwLCAmBsYIwFoAcwCGBjBZICYDMyOkSCADEqgOwBmSIAnAEak5hA==) 如果在跑make check時發現下面的回應Illegal instruction (core dumped)代表電腦不支援AVX SIMD 不知道是什麼原因造成無法順利執行舊電腦的開發環境： - `$ cat /etc/issue` Ubuntu 14.04.5 LTS \n \l - `$ cat /proc/version` Linux version 4.4.0-53-generic (buildd@lgw01-54) (gcc version 4.8.4 (Ubuntu 4.8.4-2ubuntu1~14.04.3) ) ::: - 新電腦: ```clike= gcc -c -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.c -o computepi.o gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DBASELINE -o time_test_baseline gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DOPENMP_2 -o time_test_openmp_2 gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DOPENMP_4 -o time_test_openmp_4 gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DAVX -o time_test_avx gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DAVXUNROLL -o time_test_avxunroll gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o benchmark_clock_gettime.c -o benchmark_clock_gettime time ./time_test_baseline N = 400000000 , pi = 3.141593 4.90user 0.00system 0:04.90elapsed 99%CPU (0avgtext+0avgdata 1760maxresident)k 0inputs+0outputs (0major+85minor)pagefaults 0swaps time ./time_test_openmp_2 N = 400000000 , pi = 3.141593 5.04user 0.00system 0:02.52elapsed 199%CPU (0avgtext+0avgdata 1780maxresident)k 0inputs+0outputs (0major+87minor)pagefaults 0swaps time ./time_test_openmp_4 N = 400000000 , pi = 3.141593 5.54user 0.00system 0:01.49elapsed 370%CPU (0avgtext+0avgdata 1784maxresident)k 0inputs+0outputs (0major+95minor)pagefaults 0swaps time ./time_test_avx N = 400000000 , pi = 3.141593 1.57user 0.00system 0:01.58elapsed 99%CPU (0avgtext+0avgdata 1768maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_avxunroll N = 400000000 , pi = 3.141593 1.70user 0.00system 0:01.70elapsed 99%CPU (0avgtext+0avgdata 1688maxresident)k 0inputs+0outputs (0major+85minor)pagefaults 0swaps ``` ```clike= gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DBASELINE -o time_test_baseline gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DOPENMP_2 -o time_test_openmp_2 gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DOPENMP_4 -o time_test_openmp_4 gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DAVX -o time_test_avx gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o time_test.c -DAVXUNROLL -o time_test_avxunroll gcc -O0 -std=gnu99 -Wall -fopenmp -mavx computepi.o benchmark_clock_gettime.c -o benchmark_clock_gettime for i in `seq 100 5000 25000`; do \ printf "%d," $i;\ ./benchmark_clock_gettime $i; \ done > result_clock_gettime.csv ``` :::info #### AVX [Introduction to Intel® Advanced Vector Extensions](https://software.intel.com/en-us/articles/introduction-to-intel-advanced-vector-extensions) The hardware supporting Intel® AVX (and FMA) consists of the 16 256-bit YMM registers YMM0-YMM15 and a 32-bit control/status register called MXCSR. ::: ### 使用 gnuplot 製圖 - 在 `runtime.gp` 中編輯繪圖語法 - 在 `Makefile` 中定義製圖的 rule ```clike= plot: default gnuplot runtime.gp ``` - 執行 `$ make plot` 產生`runtime.png` - 執行 `$ eog runtime.png` 看到圖像 ### 測量時間的函式 > [實作參考](https://hackmd.io/KwUwLCAmBsYIwFoAcwCGBjBZICYDMyOkSCADEqgOwBmSIAnAEak5hA==) 原來程式碼輸出的圖形(25倍時間)，明顯看出資料分佈抖動的很厲害 clock 疊代了25次，所以應將程式改為取25次運算時間之平均 ```clike= clock_gettime(CLOCK_ID, &start); for(i = 0; i < loop; i++){ compute_pi_baseline(N); } clock_gettime(CLOCK_ID, &end); printf("%lf",(double)((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/ONE_SEC))/loop); ``` #### 解法: 用信賴區間消去極端值再將X軸取logscale :::info #### 信賴區間(Confidence interval) 由樣本資料定義一段數值區間，宣稱有多少信心以估計母體的參數包含於此區間內。該數值區間上、下限稱為信賴界限(confidence limit)。用以估計的信心程度稱為信賴(心)水準(confidence level)。 ##### 標準差 ![](https://i.imgur.com/tLEWvzY.jpg) * 一般常以 95% 或 99% 為信賴水準指標；相對應的 Z分數(相差幾個標準差) 分別為 1.96 與 2.58。 * 可表示為：有 95% 信心估計母群體平均數，在樣本平均數 ± 1.96 * (母群體標準差 / 樣本數 n 的平方根) 的範圍內。而 99% 信心估計母群體平均數，則在樣本平均數 ± 2.58 * (母群體標準差 / 樣本數 n 的平方根) 的範圍內。科學符號表示方式： μ 之 95% CI = X ± 1.96 * ( σ / √n ) μ 之 99% CI = X ± 2.58 * ( σ / √n ) ::: - `clock()` ![](https://i.imgur.com/PcyQU11.png) - `clock_gettime()` ![](https://i.imgur.com/hjMtt82.png) - data 飄動, 嘗試增加 Makefile 定義 data 的資料數量: ```clike= for i in `seq 100 100 25000`; do \ printf "%d " $$i;\ ./benchmark_clock_gettime $$i; \ done > result_clock_gettime.csv ``` ![](https://i.imgur.com/KhFVewD.png) ### Error Rate :::info The higher the error rate, the less reliable the connection or data transfer will be. ::: * references: * [廖健富學長的筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I) - error rate的計算: 利用定義在`<math.h>`中的M_PI與computepi中的計算結果對比 --- ## baseline 版本 * pi 的公式 ![](https://i.imgur.com/3mUF7g8.png) * 黎曼積分原理 ![](https://hackpad-attachments.s3.amazonaws.com/charles620016.hackpad.com_M6oHPzPNlO2_p.431773_1442555734600_integral.png) * 程式碼: ```clike= 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; } ``` ## Leibniz 版本 ![](https://hackpad-attachments.s3.amazonaws.com/embedded2015.hackpad.com_Mj4fZIhTKVH_p.334595_1443717119937_pi.PNG) ### 實作 ```clike= double compute_pi_leibniz(size_t N) { double pi = 0.0; for(int i=0; i<N; i++) { int temp = (i%2) ? -1 : 1; pi += (double) temp / (2*i+1); } return pi * 4.0; } ``` ### 測試程式效能 ![](https://i.imgur.com/bRBSPOy.png) #### 測量 error rate ### 優化方向 ... - OpenMP - AVX - AVX + loop unrolling >>[實作參考](https://hackmd.io/GYZgRgnALAhgHAJgLRgCYIkqBWMMlwgwhIAMMA7AGzbBSkKpQRA=) ## Euler 版本 ![](https://i.imgur.com/32qQwAJ.png) ### 實作 ```clike= double compute_pi_euler(size_t N) { double pi = 0.0; for(int i=0; i<N; i++) pi += (double)(1 / pow(i, 2)); pi *= 6; return sqrt(pi); } ``` >>[實作參考](https://hackmd.io/GYZgRgnALAhgHAJgLRgCYIkqBWMMlwgwhIAMMA7AGzbBSkKpQRA=) ### 測試程式效能 ![](https://i.imgur.com/UfYL7Po.png) 若呼叫 `#include <math.h>` 的函式, 需要花的時間更長 #### 測量 error rate ### 優化方向 ... - OpenMP - AVX - AVX + loop unrolling ## Monte Carlo 版本 ... ### 實作 ```clike= int main() { double x,y; int count=0; /* # of points in the 1st quadrant of unit circle */ double z; double pi; count=0; for(size_t i=0; i<N; i++) { x = (double)rand()/RAND_MAX; y = (double)rand()/RAND_MAX; z = x*x+y*y; if (z<=1) count++; } printf("pi = %lf \n",(double)count/n*4); return 0; } ``` >>[實作參考](https://hackmd.io/MYFgHAJgDAjBCcBaAbAIwKZkeESzoHZhEBDYE9ECAJnjHAGYg===) ... --- ## SIMD ### 閱讀資料 [SIMD](https://docs.google.com/presentation/d/1LeVe7EAmZvqD3KN7p4Wynbd36xPOk9biBCFrdmsqzKs/edit#slide=id.p3) SIMD = Single Instruction Multiple Data，主要把原本一個指令中的數值分成多個部份，分配給多個執行緒執行完後再統整程最後結果，加速執行速度 #### SIMD 優化方式 - Auto/Semi-Auto method - 讓compiler做vectorization處理（ex: `#pragma`） - 把原本serial的執行變成平行化處理的方式 - [IR Optimization](https://web.stanford.edu/class/archive/cs/cs143/cs143.1128/lectures/14/Slides14.pdf): 對中間語言(intermediate representation)的優化 - 在從高階語言轉換產生IR後，時常會有一些冗餘的計算，應該優化 - 使用者會使用for loop來做一些型式類似並且多次的運算，減少所寫的code數量，所以compiler做loop unrolling達到優化 - Compiler Intrinsics - 利用intrinsics讓compiler編譯時，會被直接編譯成組合語言 - Specific Framework/Infrastructure - 為了平行化 - 直接利用data Parallel Framework，或是透過SIMD優化過的library - Coding in Assembly - 最極致的最佳化的方式 #### SIMD 困難的地方 - Finding Parallelism in Algorithm - 找出原本程式中可以做平行運算的部份 - Portability between different intrinsics - 不同架構中的轉換 - sse -> neon - neon -> sse - 在不同平台上就要使用不同的intrinscis，可攜性是一個問題 - Boundary handling - 資料邊界的對齊與填充 - ==Divergence== - Register Spilling - 電腦中register數量有限，所以變數會對應到register或memory上 - 編譯產生組合語言時，所使用的variable數量 > register的數量的情況發生 - ==Non-Regular Access/Processing Pattern/Dependency== - Unsupported Operations - 比較高階的函數 - Floating-Point - 浮點數運算跨裝置上的相容性問題 >>[參考概念](https://hackmd.io/BwVgxmCMBMYJwFoCmYAmBDBAWExgOA0XSQDYkRUR1TUB2MIA?both) >>[參考概念](https://hackmd.io/KYQwbMBMkAxgtAVgOwE5LwCwmAY3qogIwBm8kywAzIhCCDLrkA==?both) ### implement 效能最佳化 --- <style> h2.part {color: red;} h3.part {color: green;} h4.part {color: blue;} h5.part {color: black;} </style>

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.