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:

  ​​int clock_gettime(clockid_t clk_id, struct timespec *tp);
  

  精準度可以到達奈秒等級。
  根據函數的第一個參數，可以決定要取得的來源，若是傳入 CLOCK_REALTIME 則可能出現非單調遞增的情況，此時要用 CLOCK_MONOTONIC 或 CLOCK_MONOTONIC_RAW 才對。

  ​​struct timespec {
  ​​    time_t   tv_sec;        /* seconds */
  ​​    long     tv_nsec;       /* nanoseconds */
  ​​};
  

• time()
  回傳自 1970年01月01日 以來所經過的秒數。

• gettimeofday()
  取得目前的時間，精確度可到微秒等級。但跟 time() 一樣不一定是單調遞增，所以也不常用來量測時間。

  ​​struct timeval {
  ​​    time_t      tv_sec;     /* seconds */
  ​​    suseconds_t tv_usec;    /* microseconds */
  ​​};
  

Baseline 效能

執行 make check 得到以下結果。

另外，在 Makefile 中使用的 time 是 /usr/bin/time，而不是 bash 內建的 time。所以在輸出格式上才會有所差異。

$ 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 畫圖

每次都還要開 LibreOffice 才能生成圖片，有點慢。依照老師說的，用 gnuplot 來完成這一切，只要下個 make plot 就好。

Pi 的計算方式

Baseline 版本

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=ta{n}^{-1}1-ta{n}^{-1}0=\frac{\pi }{4}$，所以算出來要乘以4倍。

Leibniz 版本

公式

$$\frac{\pi }{4}=\sum _{n=0}^{\mathrm{\infty }}\frac{(-1{)}^n}{2n+1}$$

證明

直接相除可得

$$\frac{1}{x^2+1}=1-x^2+x^4-\cdots +(-1{)}^n{x}^{2n}+\frac{(-1{)}^{n+1}{x}^{2n+2}}{x^2+1}$$

接著對等號兩邊做定積分得

$$\frac{\pi }{4}=1-\frac{1}{3}+\frac{1}{5}-\cdots +\frac{(-1{)}^n}{2n+1}+(-1{)}^{n+1}{\int }_0^1\frac{x^{2n+2}}{x^2+1}dx$$

最後，當

$n\to \mathrm{\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}^{\mathrm{\infty }}\frac{(-1{)}^n}{2n+1}$

執行結果

$ 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

讓我出乎意料的是 AVX 並沒有像 Baseline 版本那樣，比 OpenMP 快，反而變慢了…

Euler 版本

公式

$$\frac{{\pi }^2}{6}=\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+\dots$$

執行結果

$ 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

AVX 的部份也是比 OpenMP 還要慢。

Error Rate

比較不同版本的 Error Rate

Baseline 的 error rate 與 Leibniz 的重疊了，而從圖中可以看到 Euler 的方法在大部份的時候，精準度都是比其他兩個好一點點。

同版本下 AVX Error Rate 與眾不同的問題

看到老師提及賴劭芊同學的HackMD以及王紹華同學的HackMD才發現AVX的問題，只要將 Makefile 中 for 迴圈的資料點改為可被16整除即可修正。

將

for i in `seq 1000 5000 1000000`; do

改成

for i in `seq 1008 4000 1000000`; do

原先是

修改後變成

參考資料

• Leibniz formula for π - Wikipedia
• 圓周率 - 維基百科
• ettimeofday(2) - Linux manual page
• clock_getres(2) - Linux manual page
• 廖健富學長的Hackpad
• TempoJiJi的HackMD