--- tags: ncku-course --- # 2016q3 Homework1 (compute-pi) ## 背景知識 在開始之前，先來釐清一些名詞。 ### Wall-clock time vs. CPU time * Wall-clock time 基本上是現實世界中實際經過的時間。 但有可能因為 NTP 或時區相關因素，導致不是單調遞增。 * CPU time 程式在 CPU 實際所花費的時間。 ### CPU bound vs. IO bound * CPU bound 主要的操作都是 CPU 運算，僅有少量的 IO 操作。這一類型的程式 ==有機會== 可以透過平行化的方式來加速。 * IO bound 與 CPU bound 相反，主要都是 IO 操作，無法受益於平行化處理。 ### Bourne shell 中的 time 為系統內建的工具之一，會計算程式執行的時間並回報給使用者。可以透過設定環境變數 `TIMEFORMAT` 來指定輸出格式。 預設而言，會輸出3個時間，分別是 `user`、`sys`、`real`。 * user time 在 user space 中花費的 CPU time。 * sys time 在 kernel space 中花費的 CPU time。 與 user time 相加即為 total CPU time。 * real time 也就是上方提及的 wall-clock time。不僅包含 CPU time，還包含 I/O 等等的時間。 ## 時間測量函數 * clock() 返回自程式開始執行到現在為止所使用的 CPU time，type 是 clock_t。若要換算成時間，則必須要除以 `CLOCKS_PER_SEC`。 除此之外，`clock()` 並沒有考慮到多執行緒的情況，因此會有誤差。 * clock_gettime() Prototype: ```clike= int clock_gettime(clockid_t clk_id, struct timespec *tp); ``` 精準度可以到達奈秒等級。 根據函數的第一個參數，可以決定要取得的來源，若是傳入 `CLOCK_REALTIME` 則可能出現非單調遞增的情況，此時要用 `CLOCK_MONOTONIC` 或 `CLOCK_MONOTONIC_RAW` 才對。 ```clike= struct timespec { time_t tv_sec; /* seconds */ long tv_nsec; /* nanoseconds */ }; ``` * time() 回傳自 1970年01月01日 以來所經過的秒數。 * gettimeofday() 取得目前的時間，精確度可到微秒等級。但跟 `time()` 一樣不一定是單調遞增，所以也不常用來量測時間。 ```clike= struct timeval { time_t tv_sec; /* seconds */ suseconds_t tv_usec; /* microseconds */ }; ``` ## Baseline 效能 執行 `make check` 得到以下結果。 另外，在 Makefile 中使用的 `time` 是 `/usr/bin/time`，而不是 bash 內建的 time。所以在輸出格式上才會有所差異。 ```shell= $ make check time ./time_test_baseline N = 400000000 , pi = 3.141593 3.32user 0.04system 0:03.37elapsed 99%CPU (0avgtext+0avgdata 1752maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_openmp_2 N = 400000000 , pi = 3.141593 3.72user 0.00system 0:01.86elapsed 199%CPU (0avgtext+0avgdata 1724maxresident)k 0inputs+0outputs (0major+89minor)pagefaults 0swaps time ./time_test_openmp_4 N = 400000000 , pi = 3.141593 6.63user 0.01system 0:01.71elapsed 388%CPU (0avgtext+0avgdata 1820maxresident)k 0inputs+0outputs (0major+92minor)pagefaults 0swaps time ./time_test_avx N = 400000000 , pi = 3.141593 0.97user 0.01system 0:00.98elapsed 100%CPU (0avgtext+0avgdata 1712maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps time ./time_test_avxunroll N = 400000000 , pi = 3.141593 0.94user 0.01system 0:00.95elapsed 100%CPU (0avgtext+0avgdata 1704maxresident)k 0inputs+0outputs (0major+86minor)pagefaults 0swaps ``` 觀察 `time` 輸出的時間可以知道 * 除了 OpenMP 以外的時間都是 user time ≒ real time。 表示他們屬於 CPU bound，但無法受益於平行化處理。 * 而 OpenMP 則是 user time > real time。 表示他們是屬於可以受益於平行化處理的 CPU bound。 * OpenMP 的部份有 CPU utilization 超過 100% 的情況。 CPU utilization 100% 的時候表示充份地利用了一顆核心，以我的筆電為例，有四顆核心，最高只會出現 400%。 再來是利用 LibreOffice 畫圖 ![](https://i.imgur.com/KXJ5h0G.png) 每次都還要開 LibreOffice 才能生成圖片，有點慢。依照老師說的，用 gnuplot 來完成這一切，只要下個 `make plot` 就好。 ![](https://i.imgur.com/QbGusNZ.png) ## Pi 的計算方式 ### Baseline 版本 ```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; } ``` 因為 $\int_{0}^{1} \frac{1}{x^2 + 1} dt = tan^{-1}1 - tan^{-1}0 = \frac{\pi}{4}$，所以算出來要乘以4倍。 ### Leibniz 版本 #### 公式 $$ \frac{\pi}{4} = \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} $$ #### 證明 直接相除可得 $$ \frac{1}{x^2 + 1} = 1 - x^2 + x^4 - \dots + (-1)^{n}x^{2n} + \frac{(-1)^{n+1}x^{2n+2}}{x^{2}+1} $$ 接著對等號兩邊做定積分得 $$ \frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \dots + \frac{(-1)^{n}}{2n+1} + (-1)^{n+1} \int_{0}^{1} \frac{x^{2n+2}}{x^{2}+1} dx $$ 最後，當 $n \to \infty$ 時，可以得知最後一項的積分會趨近於 0 $$ \int_{0}^{1} \frac{x^{2n+2}}{x^{2}+1} dx < \int_{0}^{1} x^{2n+2} dx = \frac{1}{2n+3} \to 0 $$ 至此，即可得證 $\frac{\pi}{4} = \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1}$ #### 執行結果 ```shell= $ make check time ./time_test_leibniz N = 400000000 , pi = 3.141593 1.68user 0.00system 0:01.68elapsed 100%CPU (0avgtext+0avgdata 1796maxresident)k 0inputs+0outputs (0major+83minor)pagefaults 0swaps time ./time_test_leibniz_openmp_2 N = 400000000 , pi = 3.141593 1.70user 0.00system 0:00.85elapsed 199%CPU (0avgtext+0avgdata 1788maxresident)k 0inputs+0outputs (0major+90minor)pagefaults 0swaps time ./time_test_leibniz_openmp_4 N = 400000000 , pi = 3.141593 3.36user 0.00system 0:00.85elapsed 392%CPU (0avgtext+0avgdata 1740maxresident)k 0inputs+0outputs (0major+91minor)pagefaults 0swaps time ./time_test_leibniz_avx N = 400000000 , pi = 3.141593 1.02user 0.00system 0:01.02elapsed 99%CPU (0avgtext+0avgdata 1936maxresident)k 0inputs+0outputs (0major+87minor)pagefaults 0swaps time ./time_test_leibniz_avxunroll N = 400000000 , pi = 3.141593 1.29user 0.00system 0:01.29elapsed 99%CPU (0avgtext+0avgdata 1784maxresident)k 0inputs+0outputs (0major+87minor)pagefaults 0swaps ``` ![](https://i.imgur.com/lDhgjj5.png) 讓我出乎意料的是 AVX 並沒有像 Baseline 版本那樣，比 OpenMP 快，反而變慢了... ### Euler 版本 #### 公式 $$ \frac{\pi^{2}}{6} = \frac{1}{1^{2}} + \frac{1}{2^{2}} + \frac{1}{3^{2}} + \dots $$ #### 執行結果 ```shell= $ make check time ./time_test_euler N = 400000000 , pi = 3.141593 1.66user 0.00system 0:01.67elapsed 99%CPU (0avgtext+0avgdata 1896maxresident)k 0inputs+0outputs (0major+92minor)pagefaults 0swaps time ./time_test_euler_openmp_2 N = 400000000 , pi = 3.141593 1.70user 0.00system 0:00.85elapsed 199%CPU (0avgtext+0avgdata 2144maxresident)k 0inputs+0outputs (0major+95minor)pagefaults 0swaps time ./time_test_euler_openmp_4 N = 400000000 , pi = 3.141593 3.38user 0.00system 0:00.85elapsed 396%CPU (0avgtext+0avgdata 2008maxresident)k 0inputs+0outputs (0major+96minor)pagefaults 0swaps time ./time_test_euler_avx N = 400000000 , pi = 3.141593 1.02user 0.00system 0:01.02elapsed 99%CPU (0avgtext+0avgdata 2092maxresident)k 0inputs+0outputs (0major+95minor)pagefaults 0swaps time ./time_test_euler_avxunroll N = 400000000 , pi = 3.141593 1.29user 0.00system 0:01.29elapsed 99%CPU (0avgtext+0avgdata 2104maxresident)k 0inputs+0outputs (0major+94minor)pagefaults 0swaps ``` ![](https://i.imgur.com/YSaYuUb.png) AVX 的部份也是比 OpenMP 還要慢。 ## Error Rate ### 比較不同版本的 Error Rate ![](https://i.imgur.com/i1Nm2EM.png) Baseline 的 error rate 與 Leibniz 的重疊了，而從圖中可以看到 Euler 的方法在大部份的時候，精準度都是比其他兩個好一點點。 ### 同版本下 AVX Error Rate 與眾不同的問題 看到老師提及[賴劭芊同學的HackMD](https://hackmd.io/s/HJrLU0dp)以及[王紹華同學的HackMD](https://hackmd.io/KYQwjARgZgTArAZgLQgJwAYpICyoMYAmSAHJiEggGzzqQEIDsqYQA===)才發現AVX的問題，只要將 Makefile 中 for 迴圈的資料點改為可被16整除即可修正。 將 ```shell= for i in `seq 1000 5000 1000000`; do ``` 改成 ```shell= for i in `seq 1008 4000 1000000`; do ``` 原先是 ![](https://i.imgur.com/oAmhhY6.png) 修改後變成 ![](https://i.imgur.com/cbc6vDZ.png) ## 參考資料 * [Leibniz formula for π - Wikipedia](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80) * [圓周率 - 維基百科](https://zh.wikipedia.org/wiki/%E5%9C%93%E5%91%A8%E7%8E%87#.E6.95.B8.E5.AD.B8.E5.88.86.E6.9E.90) * [ettimeofday(2) - Linux manual page](http://man7.org/linux/man-pages/man2/gettimeofday.2.html) * [clock_getres(2) - Linux manual page](http://man7.org/linux/man-pages/man2/clock_gettime.2.html) * [廖健富學長的Hackpad](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I) * [TempoJiJi的HackMD](https://hackmd.io/s/Hy7tssw6)