# 2017q1 Homework1 (compute-pi)
contributed by < `changyuanhua` >
### Reviewed by `baimao8437`
* commit message 開頭應為大寫。
* 可以對一開始`make check`的各 method 時間結果做一些解釋。
* 可嘗試95信賴區間解決資料抖動問題。
* Leibniz 的程式碼中 `int x = (i%2) ? -1 : 1;`這行判斷拖慢了執行速度，可考慮其他方法來解決變號，e.g.`int x = -1; x*=-1;`。
* 最後一張 error rate 的圖片其實看不到很多 method 的結果，若有重疊或接近 0 的現象應特別解釋。
* 由於此作業對於計算 $\pi$ 會使用多種 method ，使得每新增一種，就得在 makefile 中多寫很多行，可嘗試簡化 makefile 已提升測試效率，可參考[makefile修改](https://github.com/0140454/compute-pi/commit/89ab12a144d0ce6d2d3c498e50277fe81c42b8d4)。
---
## 開發環境
* 輸入指令 ` lscpu `
```
Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
每核心執行緒數：1
每通訊端核心數：4
Socket(s): 1
NUMA 節點： 1
供應商識別號： GenuineIntel
CPU 家族： 6
型號： 94
Model name: Intel(R) Core(TM) i5-6300HQ CPU @ 2.30GHz
製程： 3
CPU MHz： 799.890
CPU max MHz: 3200.0000
CPU min MHz: 800.0000
BogoMIPS: 4608.00
虛擬： VT-x
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 6144K
NUMA node0 CPU(s): 0-3
```
## 作業要求
* 在 GitHub 上 fork [compute-pi](https://github.com/sysprog21/compute-pi)，嘗試重現實驗，包含分析輸出
* 注意: 要增加圓周率計算過程中迭代的數量，否則時間太短就看不出差異
* 詳閱廖健富提供的[詳盡分析](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)，研究 error rate 的計算，以及為何思考我們該留意後者
* 用 [phonebok](/s/S1RVdgza) 提到的 gnuplot 作為 `$ make gencsv` 以外的輸出機制，預期執行 `$ make plot` 後，可透過 gnuplot 產生前述多種實做的效能分析比較圖表
* 可善用 [rayatracing](/s/B1W5AWza) 提到的 OpenMP, software pipelining, 以及 loop unrolling 一類的技巧來加速程式運作
* 詳閱前方 "Wall-clock time vs. CPU time"，思考如何精準計算時間
* 除了研究程式以外，請證明 [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)，並且找出其他計算圓周率的方法
* 請詳閱[許元杰的共筆](https://embedded2015.hackpad.com/-Ext1-Week-1-5rm31U3BOmh)
* 將你的觀察、分析，以及各式效能改善過程，並善用 gnuplot 製圖，紀錄於「[作業區](https://hackmd.io/s/Hy2rieDYe)」
## baseline 版本
```clike=
double computePi_v1(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;
}
```
* 因為積分區間為[0,1]，N為這之中的分段，所以理論上N值越大，精準度就越高
## 時間處理與 time 函式使用
* [相關資料](http://www.cnblogs.com/krythur/archive/2013/02/25/2932647.html)
* [clock() 和 time() 的回傳值資料](http://pydoing.blogspot.tw/2010/07/c-stdtime.html)
* [clock()](https://linux.die.net/man/3/clock)
```
clock_t clock(void);
```
* clock() 函数的返回值類型是 clock_t ,它除以 CLOCKS_PER_SEC 来得出时间
* clock()有3個問題：
1.超過一小時，會 overflow
2.沒有考慮 CPU 被子進程使用的情況
3.不能區分 user mode 和 kernel mode
* 回傳程式開始執行後所使用的 cpu 時間
* [clock_gettime()](https://linux.die.net/man/3/clock_gettime)
```
int clock_gettime(clockid_t clk_id, struct timespec *tp);
```
* 在計時的時候，第一個參數 clk_id 可填入(ex. CLOCK_REALTIME )，而 struct timespec* tp 則是回傳的結果
* 有一個時間結構體 timespec，計時的單位是奈秒（nanosecond）
```clike=
strace timespec {
time_t tv_sec; /* seconds */
long tv_nsec; /* nanoseconds */
}
```
* 是獲取 process 所經過的時間，不隨系统管理員修改了機器的時間而改變
* [time()](https://linux.die.net/man/2/time)
```
time_t time(time_t *t);
```
* 回傳日曆時間，也就是自 1970 年 1 月 1 日到現在所過的總秒數
* [gettimeofday()](https://linux.die.net/man/2/gettimeofday)
```
int gettimeofday ( struct timeval * tv , struct timezone * tz )
```
* 有一個時間結構體 timeval，計時的單位精確到微秒
```clike=
struct timeval {
time_t tv_sec;
suseconds_t tv_usec;
};
```
* 如果系统管理員修改了機器的時間(調整日期)，則所得的時間也會改變
* 四種函數比較
* clock() 的精確度是 10 毫秒 (ms)
* time() 的精確度是 10 毫秒 (ms)
* gettimofday() 的精確度是微秒 (μs)
* clock_gettime() 的計量單位是奈秒 (ns)
* 為什麼 clock_gettime() 結果飄忽不定？
* 會因爲不同的 process 在CPU之間切換而受到影響
* 為什麼 time() 和 gettimeofday() 不適合拿來作 benchmark ?
* 因為回傳的是系統時間，如果在執行期間改變的話，得到的值就不正確
## [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操作，使得平行處理無法帶來顯著的效能提昇 ( 或沒有提昇)。
## 操作
* 編譯指令 ```make check```
* 輸出
```
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
```
* Makefile 裡頭有` -fopenmp` `-mavx `，是作為指定開啟 OpenMP 和 AVX
```
CFLAGS = -O0 -std=gnu99 -Wall -fopenmp -mavx
```
* baseline
```
time ./time_test_baseline
N = 400000000 , pi = 3.141593
1.25user 0.00system 0:01.25elapsed 99%CPU (0avgtext+0avgdata 1776maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps
```
* openmp_2
```
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
1.41user 0.00system 0:00.70elapsed 199%CPU (0avgtext+0avgdata 1868maxresident)k
0inputs+0outputs (0major+87minor)pagefaults 0swaps
```
* openmp_4
```
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
1.55user 0.00system 0:00.41elapsed 373%CPU (0avgtext+0avgdata 1808maxresident)k
0inputs+0outputs (0major+89minor)pagefaults 0swaps
```
* avx
```
time ./time_test_avx
N = 400000000 , pi = 3.141593
0.62user 0.00system 0:00.62elapsed 99%CPU (0avgtext+0avgdata 1704maxresident)k
0inputs+0outputs (0major+81minor)pagefaults 0swaps
```
* avxnuroll
```
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
0.44user 0.00system 0:00.44elapsed 99%CPU (0avgtext+0avgdata 1856maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps
```
* 編譯指令 ```make gencsv```
* 輸出
```shell
for i in `seq 100 5000 25000`; do \
printf "%d," $i;\
./benchmark_clock_gettime $i; \
done > result_clock_gettime.csv
```
* 可視為
```cstyles=
for(i=100; i<25000; i+=5000)
{
printf("%d",i);
benchmark_clock_gettime(i);
}
```
* 編譯指令 ```time ./time_test_baseline```
* 輸出
```
N = 400000000 , pi = 3.141593
real 0m1.276s
user 0m1.272s
sys 0m0.000s
```
* 顯示了三種時間：real,user和sys
* real time : 也就是 Wall-clock time，當然若有其他程式同時在執行，必定會影響到。
* user time : 表示程式在 user mode 佔用所有的 CPU time 總和。
* sys time : 表示程式在 kernel mode 佔用所有的 CPU time 總和 ( 直接或間接的系統呼叫 )
* 我將 seq 100 5000 25000 改為 seq 500 500 200000 作圖後得
* 從圖中發現 Baseline>openmp_2>avx>avx_unroll>openmp_4
## 實驗一 : Leibniz formula for π
* **proof**
* Considering only the integral in the last line, we have:
* Therefore, by the squeeze theorem, as n → ∞ we are left with the Leibniz series:
* [Leibniz](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)
* **Leibniz**
```clike=
double compute_pi_leibniz(size_t N)
{
double pi = 0.0;
for(int i=0; i<N; i++) {
int x = (i%2) ? -1 : 1;
pi += (double)x / (2*i +1) ;
}
return pi * 4.0;
}
```
>請注意程式碼縮排
>[name=課程助教][color=red]
>>ok, 謝謝[name=changyuanhua][color=black]
* **Leibniz AVX SIMD**
```clike=
double compute_pi_leibnizavx(size_t N)
{
double pi = 0.0;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
ymm0 = _mm256_set_pd(1.0,-1.0,1.0,-1.0);
ymm1 = _mm256_set1_pd(2.0);
ymm2 = _mm256_set1_pd(1.0);
ymm4 = _mm256_setzero_pd();
for (int i = 0; i <= N - 4; i += 4) {
ymm3 = _mm256_set_pd(i, i+1.0, i+2.0, i+3.0);
ymm3 = _mm256_mul_pd(ymm3,ymm1);
ymm3 = _mm256_add_pd(ymm3,ymm2);
ymm3 = _mm256_div_pd(ymm0,ymm3);
ymm4 = _mm256_add_pd(ymm4,ymm3);
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4);
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}
```
* **Leibniz AVX SIMD + Unroll looping**
```clike=
double compute_pi_leibnizavxunroll(size_t N)
{
double pi = 0.0;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4,
ymm5, ymm6, ymm7, ymm8, ymm9,
ymm10;
ymm0 = _mm256_set_pd(1.0,-1.0,1.0,-1.0);
ymm1 = _mm256_set1_pd(2.0);
ymm2 = _mm256_set1_pd(1.0);
ymm7 = _mm256_setzero_pd();
ymm8 = _mm256_setzero_pd();
ymm9 = _mm256_setzero_pd();
ymm10 = _mm256_setzero_pd();
for (int i = 0; i <= N - 16; i += 16) {
ymm3 = _mm256_set_pd(i , i+1.0 , i+2.0, i+3.0);
ymm4 = _mm256_set_pd(i+4.0 , i+5.0 , i+6.0, i+7.0);
ymm5 = _mm256_set_pd(i+8.0 , i+9.0 , i+10.0, i+11.0);
ymm6 = _mm256_set_pd(i+12.0 , i+13.0 , i+14.0, i+15.0);
ymm3 = _mm256_mul_pd(ymm1, ymm3);
ymm4 = _mm256_mul_pd(ymm1, ymm4);
ymm5 = _mm256_mul_pd(ymm1, ymm5);
ymm6 = _mm256_mul_pd(ymm1, ymm6);
ymm3 = _mm256_add_pd(ymm2, ymm3);
ymm4 = _mm256_add_pd(ymm2, ymm4);
ymm5 = _mm256_add_pd(ymm2, ymm5);
ymm6 = _mm256_add_pd(ymm2, ymm6);
ymm3 = _mm256_div_pd(ymm0, ymm3);
ymm4 = _mm256_div_pd(ymm0, ymm4);
ymm5 = _mm256_div_pd(ymm0, ymm5);
ymm6 = _mm256_div_pd(ymm0, ymm6);
ymm7 = _mm256_add_pd(ymm7, ymm3);
ymm8 = _mm256_add_pd(ymm8, ymm4);
ymm9 = _mm256_add_pd(ymm9, ymm5);
ymm10 = _mm256_add_pd(ymm10, ymm6);
}
double tmp1[4] __attribute__((aligned(32)));
double tmp2[4] __attribute__((aligned(32)));
double tmp3[4] __attribute__((aligned(32)));
double tmp4[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp1, ymm7);
_mm256_store_pd(tmp2, ymm8);
_mm256_store_pd(tmp3, ymm9);
_mm256_store_pd(tmp4, ymm10);
pi += tmp1[0] + tmp1[1] + tmp1[2] + tmp1[3] +
tmp2[0] + tmp2[1] + tmp2[2] + tmp2[3] +
tmp3[0] + tmp3[1] + tmp3[2] + tmp3[3] +
tmp4[0] + tmp4[1] + tmp4[2] + tmp4[3];
return pi * 4.0;
}
```
* **使用 gnuplot 作圖分析**
* 從圖中發現 Baseline > leibniz > leibniz avx > openmp_2 > avx > leibniz avx unroll > avx_unroll > openmp_4
* leibniz avx 和 openmp_2 看起來相近
## 實驗二 : Euler for π
* **code**
```clike=
double compute_pi_euler(size_t N)
{
double pi = 0.0;
for (size_t i = 1; i <= N; i++)
pi += (1 / pow(i, 2.0));
pi *= 6;
return sqrt(pi);
}
```
* 之前程式碼因為 pow() 無法使用，所以就用最直接的方式表示之，但 sqrt() 不知道該怎呈現，所以找了資料，在 gcc 後加入 -lm 就可以解決就算加入 math.h 也不能使用 pow() 和 sqrt() 的問題
* **使用 gnuplot 作圖分析**
* 從圖中發現 euler > Baseline > leibniz > leibniz avx > openmp_2 > avx > leibniz avx unroll > avx_unroll > openmp_4
## 實驗三 : Nilakantha for π
* **code**
```clike=
double compute_pi_nilakantha(size_t N)
{
double pi = 0.0;
int x = 1;
for(size_t i=1; i <= N; i++){
pi += x / ((2.0 * (double)i)*(2.0*(double)i+1.0)*(2.0*(double)i+2.0));
x *= (-1);
}
pi *= 4.0;
pi += 3;
return pi;
}
```
* **使用 gnuplot 作圖分析**
* 從圖中發現 euler > nilakantha > Baseline > leibniz > leibniz avx > openmp_2 > avx > leibniz avx unroll > avx_unroll > openmp_4
## error rate
* **code**
```clike=
// Error rate calculation
#define M_PI acos(-1.0)
double pi = compute_pi(n);
double diff = pi - M_PI > 0 ? pi - M_PI : M_PI - pi;
double error = diff / M_PI;
```
* 觀看 [廖健富學長的筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)
## 參考資料
* [TempoJiJi](https://hackmd.io/KwBgpgxghhAcsFpYgGyICwEYXqZgJgEYL4DMmYAnAEz4Bmm41QA=?view)
* [shelly4132](https://hackmd.io/s/HJrLU0dp#)
* [compute pi](https://en.wikipedia.org/wiki/Pi)
* [廖健富學長](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)
* [clock gettime()資料](https://starforcefield.wordpress.com/2012/08/06/c%E8%AA%9E%E8%A8%80%EF%BC%9A%E4%BD%BF%E7%94%A8clock_gettime%E8%A8%88%E7%AE%97%E7%A8%8B%E5%BC%8F%E7%A2%BC%E7%9A%84%E6%99%82%E9%96%93%E9%9C%80%E6%B1%82/)
Sign in with Wallet
Connect another wallet