# 2016q3 Homework1(compute_pi)
## 目標
* 學習透過離散微積分求圓周率
* Leibniz formula for π
* 積分計算圓周率π
* Function
* 著手透過 SIMD 指令作效能最佳化
* [作業說明](https://hackmd.io/JwdgZmCGBsDGCmBaAzPWAGRAWWtqOABNkAORESAVnUqzACMxkx0g#a03-compute-pi)
## Makefile
每寫一份功課又是學習Makefile的時候=w=
* 編譯時有`-DBASELINE`等等，其中的`-D`指的是`#define BASELINE 1`
* 這樣程式碼就可以塞在同一個.c .h中！只要在編譯的時候加上`-D`就好了！！
* Makefile的篇幅可以大幅降低！
* 我想phonebook也可以這樣改寫，不然不同的方法就還要在多生一組.c .h
## 實驗
* `make check`結果：
```shell==
time ./time_test_baseline
N = 400000000 , pi = 3.141593
7.59user 0.00system 0:07.59elapsed 100%CPU (0avgtext+0avgdata 1560maxresident)k
8inputs+0outputs (0major+82minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
8.03user 0.00system 0:04.01elapsed 200%CPU (0avgtext+0avgdata 1572maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
8.32user 0.00system 0:02.16elapsed 384%CPU (0avgtext+0avgdata 1676maxresident)k
0inputs+0outputs (0major+89minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
3.30user 0.00system 0:03.30elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
1.73user 0.00system 0:01.73elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps
```
* `make gencsv` ==> 輸出至表單
輸出結果與在Makefile中看到的程式碼只差在`$$i`
以下為輸出：
```shell==
for i in `seq 100 5000 25000`; do \
printf "%d," $i;\
./benchmark_clock_gettime $i; \
done > result_clock_gettime.csv
```
原來除了寫進txt檔之外，還可以寫進libre office的表單裡，好酷！
* 理解程式碼：
* 上面的程式碼相當於：
```C==
for(i=100; i<25000; i+=5000)
{
printf("%d",i);
benchmake_clock_gettime(i);
}
```
然後把結果輸出至.csv檔
:::warning
細節的語法待會查 先highlight
:::
* 在Makefile中，在執行各檔時前面多加了`time`
* time: 查看指令的執行時間
* 直接執行`time ./time_test_baseline`
```
N = 400000000 , pi = 3.141593
real 0m7.805s
user 0m7.808s
sys 0m0.000s
```
* 顯示了三種時間：real,user和sys
Real refers to actual elapsed time; User and Sys refer to CPU time used only by the process.
* **real**: real time,a wall clock time,指令從開始到結束執行的時間，含其他process
* **user**: user CPU time,指令在user mode的時間,其他process的時間不算
* **sys**: system CPU time,指令在kernel mode的時間,其他process的時間不算
* `real = user + sys`
* ==**並非如此**==，如sleep(),不論在user mode或kernel mode沒佔用CPU時間，那就不會被算入，反觀real time仍會把sleep()的時間算進去
### Baseline
為上述積分式的實作，最後回傳的值要*4，因為算出來是 $π \over4$
* N愈大[0,1]就被切割成更多等份，dt也就愈小
```C==
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;
}
```
### OpenMP
* 也是用跟baseline一樣的積分式求得π
* 可以傳遞參數決定要開幾個threads
* 在`time_test.c`中分別傳入2個和4個threads的實驗
* 程式碼：
* private(x)：
讓每個thread都有自己的x,不然threads之間會把x當成共用的變數
* reduction(+:pi)：讓各個thread針對指定的變數擁有一份有起始值的複本
* 支援的運算元有 +, *, –, &, ^, |, &&, ||；而變數則必須要是 shared 的
```C==
double compute_pi_openmp(size_t N, int threads)
{
double pi = 0.0;
double dt = 1.0 / N;
double x;
#pragma omp parallel num_threads(threads)
{
#pragma omp for private(x) reduction(+:pi)
for (size_t i = 0; i < N; i++) {
x = (double) i / N;
pi += dt / (1.0 + x * x);
}
}
return pi * 4.0;
}
```
### AVX SIMD
* 甚麼是AVX？
* AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集
* AVX指令集借鑑了一些AMD SSE5的設計思路，進行擴展和加強，形成一套新一代的完整SIMD指令集規範。
* SSE將向量處理能力從64位提升到128位後，SSE系列指令都只能使用128位XMM暫存器，這次AVX將所有16個128位XMM暫存器擴充為256位的YMM暫存器，從而支持256位的向量計算。
* prototype:
* _mm256 : single precision,32-bit
* _mm256d : double precision, 64bit
配上註解與下方的連結終於稍微看懂了~
* `__attribute__ `機制，可以用來設置函數屬性（Function Attribute）、變數屬性（Variable Attribute）和類型屬性（Type Attribute）。
* `aligned` 則規定變數或結構的最小對齊格式，以 Byte 為單位。
>不知道為甚麼沒有highlight@@
```C==
double compute_pi_avx(size_t N)
{
double pi = 0.0;
double dt = 1.0 / N;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
ymm0 = _mm256_set1_pd(1.0);
ymm1 = _mm256_set1_pd(dt);
ymm2 = _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0);
ymm4 = _mm256_setzero_pd(); // sum of pi
for (int i = 0; i <= N - 4; i += 4) {
ymm3 = _mm256_set1_pd(i * dt); // i*dt, i*dt, i*dt, i*dt
ymm3 = _mm256_add_pd(ymm3, ymm2); // x = i*dt+3*dt, i*dt+2*dt, i*dt+dt, i*dt+0.0
ymm3 = _mm256_mul_pd(ymm3, ymm3); // x^2 = (i*dt+3*dt)^2, (i*dt+2*dt)^2, ...
ymm3 = _mm256_add_pd(ymm0, ymm3); // 1+x^2 = 1+(i*dt+3*dt)^2, 1+(i*dt+2*dt)^2, ...
ymm3 = _mm256_div_pd(ymm1, ymm3); // dt/(1+x^2)
ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += dt/(1+x^2)
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4); // move packed float64 values to 256-bit aligned memory location
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}
```
* 程式碼解讀(可能有誤！！)
* `__m256d`:存放double,256-bit可用，因此可以存放256/64 = 4 個double
* 每次的運算是4個double **個別** 加總
* 最後使用tmp[4]把4個值存起來，最後在加總
* `(aligned(32))` 是因為double有8個byte,又宣告 double tmp[4] ==> 8*4= 32(byte)
### Leibniz
* 程式：
```C==
double compute_pi_Leibniz(size_t N)
{
double pi = 0.0;
double x;
for (size_t i = 0; i < N; i++) {
x = pow(-1,i) / (2*i +1);
pi += x;
}
return pi * 4.0;
}
```
一直出現錯誤資訊：
應該是link的問題，但是即使加了`-lm`還是沒用@@
:::warning
computepi.c:(.text+0xd80): 未定義參考到 pow
collect2: error: ld returned 1 exit status
make: *** [default] Error 1
:::
> computepi.c 有included <math.h>嗎？ 或許可以試試，我花很多時間處理這個...[name=SarahYuHanCheng]
* 只好先改成
```C==
int tmp = (i%2) ? -1 : 1;
x = (double)tmp / (2*i +1) ; //記得轉型
pi += x;
```
### Leibniz AVX
看了一些資料在搭配原有的程式碼，很快就可以寫出來了~ 覺得開心=w=
```C==
double compute_pi_leibizavx(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); //for pow(-1,i)
ymm1 = _mm256_set1_pd(2.0);
ymm2 = _mm256_set1_pd(1.0);
ymm4 = _mm256_setzero_pd(); // sum of pi
for (int i = 0; i <= N - 4; i += 4) {
ymm3 = _mm256_set_pd(i, i+1.0, i+2.0, i+3.0); //i i+1 i+2 i+3
ymm3 = _mm256_mul_pd(ymm1, ymm3); //2*i 2*(i+1)...
ymm3 = _mm256_add_pd(ymm3,ymm2); //2*i+1 2*(i+1)+1 ...
ymm3 = _mm256_div_pd(ymm0,ymm3); //(-1)^i/(2*i+1) ...
ymm4 = _mm256_add_pd(ymm3,ymm4); //sum
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4); // move packed float64 values to 256-bit aligned memory location
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}
```
### 其他計算圓周率的方法
* [蒙地卡羅法](https://zh.wikipedia.org/wiki/%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85%E6%96%B9%E6%B3%95)：
* 邊長為2的正方形，正方形內有一個內切圓，面積為 $1 \cdot 1 \cdotπ = π$
* 圓面積和正方形面積之比為π:4 = 正方形內均勻隨機丟石頭，落在圓形內的機率
```C==
double compute_pi_MonteCarlo(size_t N)
{
int count = 0;
srand(time(NULL));
for (size_t i = 0; i < N; i++) {
double x = (double)rand()/RAND_MAX;
double y = (double)rand()/RAND_MAX;
if((x*x+y*y)<=1) count++;
}
return 4.0 * count / N ;
}
```
### Plot
* 複習語法
* `using [a：b]` : 以第a column當x座標,第b column當y座標
製圖時產生的warning:
:::warning
warning: Skipping data file with no valid points
:::
查了一段時間，玄機就在副檔名`.csv`和gnuplot他吃數據的方式啊！！
**csv : comma separated value**
所以資料之間是以 ==**逗號**== 區隔，因此要在gnuplot script裡加一行`set datafile separator ","`才讀的到資料
* 不同方法比較圖
:::warning
抖的超厲害！
:::
取樣從N=1000~N=250000，間隔1000
隔一天跑出來的結果差好多@@
omp_4沒那麼劇烈了，還不知道為甚麼會這樣
**新增Leibniz**
* 可以發現Leibniz與omp_2的時間差不多，但又比omp_2還穩定
* LeibnizAVX也比avx的時間再更短
* LeibnizAVX_unroll語omp_4幾乎重合
### 思考如何精準計算時間
* walk clock time: 由xtime記錄，系統每次啟動時會先從設備上的 [RTC](https://zh.wikipedia.org/wiki/%E5%AF%A6%E6%99%82%E6%99%82%E9%90%98) 上讀入 xtime
* xtime:從cmos電路中取得的時間，一般是從某一歷史時刻開始到現在的時間
* Jeffie:紀錄系統開幾以來，經過多少個 tick，取決於系統的頻率，單位是Hz
* system time:
* 系統上的時間，開機時它會讀取RTC 來設定，可能過程中還會有時區換算等之類的設置。
### 信賴區間
## 參考資料
[TempoJiJi hackmd](https://hackmd.io/KwBgpgxghhAcsFpYgGyICwEYXqZgJgEYL4DMmYAnAEz4Bmm41QA=?view)
* Time
* [time_ref1](http://yuanfarn.blogspot.tw/2012/08/linux-time.html)
* [time_ref2](http://stackoverflow.com/questions/556405/what-do-real-user-and-sys-mean-in-the-output-of-time1)
* [容易混淆LINUX時鐘的xtime和jiffies](http://www.coctec.com/docs/linux/show-post-186188.html)
* OpenMP
* [openmp_private](http://viml.nchc.org.tw/blog/paper_info.php?CLASS_ID=1&SUB_ID=1&PAPER_ID=17)
* [openmp_reduction](https://kheresy.wordpress.com/2006/09/22/%E7%B0%A1%E6%98%93%E7%9A%84%E7%A8%8B%E5%BC%8F%E5%B9%B3%E8%A1%8C%E5%8C%96%EF%BC%8Dopenmp%EF%BC%88%E4%BA%94%EF%BC%89-%E8%AE%8A%E6%95%B8%E7%9A%84%E5%B9%B3%E8%A1%8C%E5%8C%96/)
* AVX
* [introduction](http://www.expreview.com/13236.html)
* [manual](http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-optimization-manual.pdf)(chap11)
* [functions](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=579,2050,100,4629,4578,112,100&techs=AVX)
###### tags: `system embedded` `HW1-3`
Sign in with Wallet
Connect another wallet