2016q3 Homework1 (compute-pi)

contributed by <HaoTse>

---

SIMD指令集

• AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集。
• 我的電腦 CPU 為 intel i5-4200U
  
• 使用時需#include <immintrin.h>，並在 Makefile 中的編譯選項加入-mavx
• Intel Intrinsics Guide

---

OpenMP

在 raytracing 中使用過的 OpenMP 也可以應用在 compute-pi 當中。

不過在 compute-pi 中每次迭代都會讀取 pi 值並更新 pi 值，不同次的迴圈中並不相互獨立，有兩個方法可以避免：

1. 可以使用 atmoic、critical，防止變數同時被多個執行緒修改，但速度就會慢上許多，因為一次只允許一個執行續進行存取。
2. 使用 reduction，語法為 reduction(<op>:<variable>)。

   作業採用方法2 鄭皓澤

---

Wall-clock time vs. CPU time

---

Wall-clock time

• 顧名思義就是掛鐘時間，也就是現實世界中實際經過的時間，是由 kernel 裡的 xtime 來紀錄，系統每次啟動時會先從設備上的 RTC 上讀入 xtime。
  • xtime 精確度到微秒
  • xtime vs. jiffies
• System time：指系統上的時間，開機時它會讀取 RTC 來設定，可能過程中還會有時區換算等之類的設置。一般說來 System time 就是我們執行 date 命令看到的時間。
• Wall-clock time 可能會和 NTP 伺服器、時區、日光節約時間同步或使用者自己調整。所以「通常」不建議拿來量測程式時間，因為它不是一個穩定的時間，用專業一點的用語講，Wall-clock time 不一定是單調遞增 (monotonic)。

---

CPU time

• 程序在 CPU 上面運行消耗 (佔用) 的時間，clock() 就是很典型用來計算 CPU time 的時間函式。
  • 注意. 如果有多條 threads，那 CPU time 算的是「每條 thread 的使用時間加總」，所以 CPU time 可能比 Wall-clock time 還大。
• 執行時間除了 CPU time 以外，可能還包括 I/O time、 communication channel delay、synchronization overhead…等等。

---

time 指令

今日的 UNIX 系統內事實上存在兩個 time 指令，一個是系統所提供的 /usr/bin/time，這個 time 指令是 System V 版本所提供的工具，這是為了 Bourne shell 使用者所寫的工具；而另一個則是 C shell 內建指令 time。

因此 $ make check 輸出格式與直接 $ time ./程式名稱 不同。鄭皓澤

1. real time : 也就是 Wall-clock time，當然若有其他程式同時在執行，必定會影響到。
2. user time : 表示程式在 user mode 佔用所有的 CPU time 總和。
3. sys time : 表示程式在 kernel mode 佔用所有的 CPU time 總和 ( 直接或間接的系統呼叫 )。

---

timer 比較

clock()

• #include <time.h>
• 回傳程式執行後佔用的 CPU 時間，單位為 tick 。

time()

• #include <time.h>
• 回傳日曆時間，也就是自 1970 年 1 月 1 日後經過的總秒數。

clock_gettime()

• #include <time.h>

• 編譯時連結 librt.so
  gcc編譯時加入 -Wl,-lrt 選項，或連結時加入 -lrt 選項。

• 函式原型

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

  struct timespec* tp 是函式回傳結果，struct timespec 宣告如下

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

• 精確到奈秒等級。

gettimeofday()

• #include<sys/time.h>

• 函式原型

  ​int gettimeofday(struct  timeval*tv,struct  timezone *tz )
  

---

分析

首先直接將 Makefile gencsv 中的 for loop 改為 seq 1000 5000 1000000，得到下圖

這張圖長的有點醜，這裡研究了 Libre Office 很久，最後直接使用 gnuplot 畫出來，如下圖

之後才發現 Libre Office 中的 Chart Type 要選 XY Scatter。鄭皓澤

注意. gnuplot 讀入 .csv檔時，","後面要加上空白或將他改為空白。鄭皓澤

---

信賴區間

在 trace benchmark_clock_gettime.c 時，發現裡面有許多 for loop 去跑計算 pi 的 function，回想起觀看廖健富學長筆記時老師提到的信賴區間，將這裡的 for loop 視為取樣，並將信賴區間紀錄起來。

計算信賴區間函式如下

double compute_ci(double *min, double *max, double data[SAMPLE_SIZE])
{
    double mean = 0.0;
    double std_dev = 0.0; //standard deviation
    double std_err;
    int i = 0;

    //calculate mean value
    for(i = 0; i < SAMPLE_SIZE; i++)
    {
        mean += data[i];
    }
    mean /= (double)SAMPLE_SIZE;

    //calculate standard deviation
    for(i = 0; i < SAMPLE_SIZE; i++)
    {
        std_dev += (data[i] - mean) * (data[i] - mean);
    }
    std_dev = sqrt(std_dev / (double)SAMPLE_SIZE);

    //calculate standard deviation
    std_err = std_dev / sqrt((double)SAMPLE_SIZE);

    *min = mean - (1.96 * std_err);
    *max = mean + (1.96 * std_err);

    return mean;
}

如何使用信賴區間使得數據更好看還在研究了不少時間，當初機統應該好好學的QQ，下圖為使用信賴區間跑的成果

不過用 thread 的兩個版本還是有凸出，原因還在探討中。

---

error rate

目前只能看出花費時間的比較，無法比較計算出來的 pi 的正確性，因此決定再畫一張圖做比較。首先定義 pi 在 double 型態下的正確值，在<math.h>中有定義了M_PI可以用做比較。加入以下 macro

#define compute_error(pi) ( (pi > M_PI) ? pi - M_PI : M_PI - pi)

• 輸出結果
  
  這的結果與我預期的不一樣，我原本認為用同一個算式的 error rate 應該都相同，但在這張圖中我卻發現使用 AVX + unroll looping 的版本在 1000、11000、21000…中的數字差異甚大，原本我認為是我更動到了程式碼，一直找不到原因下，我將原本的程式碼 clone 下來，將 time_test.c 中的 N 改為 1000，$ make check 後發現 AVX + unroll looping 版本仍然輸出差異比別人，至於為什麼會這樣還在思考中。

---

Leibniz

原本的 baseline 算法為

$$\pi =4{\int }_0^1\frac{1}{1+x^2}dx$$

Leibniz 級數

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

• 證明
  

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

  等號兩邊做定積分

  $\frac{\pi }{4}={\int }_0^1\frac{1}{1+x^2}dx=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}=1-\frac{1}{3}+\frac{1}{5}-\cdots +\frac{(-1{)}^n}{2n+1}=\sum _{i=1}^{\mathrm{\infty }}\frac{(-1{)}^n}{2n+1}$

  故得證

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

---

結果分析

將 Leibniz 公式轉成程式碼如下

double compute_pi_leibniz(size_t N)
{
    double pi = 0.0;
    for(size_t i = 0; i < N; i++)
    {
        int sign = (i % 2) ? -1 : 1;
        pi += (sign / (2.0 * (double)i +1.0));
    }

    return pi * 4.0;
}

得到結果

• runtime
  
  這裡看出 Leibniz 的速度比 baseline 還快，我猜測是因為 Leibniz 的乘除運算較少。
• error_rate
  

---

參考資料

• SIMD wiki
• time 指令輸出格式
• timer
  • c timer 整理
  • clock_gettime
  • gettimeofday
• gnuplot
• standard error
• LaTex
• 信賴區間的迷思
• 廖健富學長筆記

tags: HaoTse sysprog21 compute-pi SIMD AVX timer Leibniz