# 2019q3 Homework4 (compute-pi) contributed by < `xl86305955` > ###### tags:`sysprog2019` `Homework4` `compute-pi` ## 目標 * 複習數值系統 * 觀察和使用 SIMD 指令 * 練習透過離散微積分求圓周率 * [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80) * [積分計算圓周率π](http://book.51cto.com/art/201506/480985.htm) * 著手透過 SIMD 指令作效能最佳化 ---- ## 開發環境 * 作業系統 * Ubuntu 18.04.3 LTS * CPU ``` $ lscpu Architecture: x86_64 CPU op-mode(s): 32-bit, 64-bit Byte Order: Little Endian CPU(s): 12 On-line CPU(s) list: 0-11 Thread(s) per core: 2 Core(s) per socket: 6 Socket(s): 1 NUMA node(s): 1 Vendor ID: GenuineIntel CPU family: 6 Model: 158 Model name: Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz Stepping: 10 CPU MHz: 799.113 CPU max MHz: 4600.0000 CPU min MHz: 800.0000 BogoMIPS: 6384.00 Virtualization: VT-x L1d cache: 32K L1i cache: 32K L2 cache: 256K L3 cache: 12288K NUMA node0 CPU(s): 0-11 ``` ---- ## 有關 OpenMP 與 AVX ### OpenMP `OpenMP` 是一種 API，使程式開發者能夠過這類的指令達成程式的平行執行，其中就是透過多執行緒來實現 當然，內建的 thread library 如 PThreads 也可以達成透過多執行緒達成程式的平行分解，但是若是僅要對一個 for 迴圈分解，使用 thread library 是否有點大材小用？除此之外，OpenMP 也據可跨平台的特性 OpenMP 語法分為三類 * Directives * Clauses * Functions 使用語法大致如下： ``` #pragma omp <directive> [clause[[,] clause] ...] ``` 使用上僅需 `#include <omp.h>` 並在欲平行化的程式前加入 `pragma omp <directive> [clause[[,] clause] ...]` 即可，並在編譯選項加入 `-fopenmp`，使用上相當簡單 以本次計算 baseline 計算 pi 之 OpenMP 版本為例： * 在 `line 6` 中，代表下面程式將被 `num_threads` 所指定的執行緒數目平行處理 * 在 `line 9` 中，將 x 變數設為各 threads 私有 （OpenMP default 設定為 public，這樣可能會有 data race 狀況發生），並對每條 threads 對 pi 加法的結果合併加入到主 threads 之中 ```clike= 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 AVX（Advanced Vector Extensions）是一種由 intel 所提出來的指令擴充，可藉此實現 SIMD compute_pi 所使用到的程式碼功能解析： * __m256d _mm256_set1_pd (double a) * 將一個長度 64 位元的浮點數複製四份儲存成長度 256 位元 * 行為： ``` FOR j := 0 to 3 i := j*64 dst[i+63:i] := a[63:0] ENDFOR dst[MAX:256] := 0 ``` * __m256d _mm256_set_pd (double e3, double e2, double e1, double e0) * pd 代表 packed double，將 4 個長度 64 位元的數視為一個 256 位元數值儲存 * 行為： ``` dst[63:0] := e0 dst[127:64] := e1 dst[191:128] := e2 dst[255:192] := e3 dst[MAX:256] := 0 ``` * __m256d _mm256_setzero_pd (void) * 將暫存器設為零 * 行為： ``` dst[MAX:0] := 0 ``` * _mm256_add_pd (__m256d a, __m256d b) * 將暫存器內的資料切割成 4 等份，對應相同位置相加 * 行為： ``` FOR j := 0 to 3 i := j*64 dst[i+63:i] := a[i+63:i] + b[i+63:i] ENDFOR dst[MAX:256] := 0 ``` * __m256d _mm256_mul_pd (__m256d a, __m256d b) * 對應位置相乘 * 行為： ``` FOR j := 0 to 3 i := j*64 dst[i+63:i] := a[i+63:i] * b[i+63:i] ENDFOR dst[MAX:256] := 0 ``` * __m256d _mm256_div_pd (__m256d a, __m256d b) * 對應位置相除 * 行為： ``` FOR j := 0 to 3 i := 64*j dst[i+63:i] := a[i+63:i] / b[i+63:i] ENDFOR dst[MAX:256] := 0 ``` * _mm256_store_pd (double * mem_addr, __m256d a) * 將資料儲存到指定記憶體位置 * 行為： ``` MEM[mem_addr+255:mem_addr] := a[255:0] ``` 參考:[Intel Intrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#techs=AVX&expand=4922,4973,4922,4973,5582,124,5016,3922,2150) ```clike= 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; } ``` ### AVX + Loop Unrolling ## 測試程式效能 ### baseline $公式：\frac{\pi}{4} =\int_0^1 \frac 1{1+x^2} \, dx$ ![](https://i.imgur.com/7fmEVBO.png) ```=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; } ``` 在 Makefile 中，可以根據 `METHOD ?= ` 選項，對相對應計算方法繪圖 畫圖之前，需先從 Makefile 選項執行`$ make gencsv`檔案生成 `result_clock_gettime.csv` 檔 其中測試從將 N 切割成 1008 份開始，每次增加 4000 直到 1000000 為止，從 .csv 檔中可以很清楚得到當 N 切割的愈細，所需執行時間愈久 ```= gencsv: default for i in `seq 1008 4000 1000000`; do \ printf "%d " $$i;\ ./benchmark_clock_gettime $$i; \ done > result_clock_gettime.csv ``` 接著，利用 `$ make plot` 透過 `gnuplot` 將圖匯出 gnuplot 會依據在 `/scripts` 下的 `runtime.gp` 檔，利用 gencsv.csv 欄位內容畫出圖表 ``` plot: gencsv gnuplot scripts/runtime.gp ``` ![](https://i.imgur.com/Zswud1Z.png) * Error rate ![](https://i.imgur.com/QHZN2Zm.png) ### LEIBNIZ $公式：\frac{\pi}4 = \sum_{k=0}^\infty\frac{(-1)^k}{2k+1}$ #### proof $tan\theta = 1$ $\theta = \frac{\pi}{4}= arctan(1)$ 其中 $\forall x\in [-1,1]\quad{arctan} (x)=\sum_{k=0}^{\infty} (-1)^k\frac{x^{2k+1}}{2k+1}= x - \frac13 x^3 + \frac15 x^5 - \frac17 x^7+ \cdots$ $arctan'(x) = \frac{1}{1+x^2}$ $\begin{split}\\ arctan(1) &= \int_0^1\frac{1}{1+x^2}dx\\ &= \int_0^1(\sum_{k=0}^{n}(-1)^kx^{2k}+\frac{(-1)^{n+1}x^{2n+2}}{1+x^2})dx\\ &= \int_0^1\sum_{k=0}^{n}(-1)^kx^{2k}dx+\int_0^1\frac{(-1)^{n+1}x^{2n+2}}{1+x^2}dx\\ &=\sum_{k=0}^{n}\frac{(-1)^k}{2k+1}+\int_0^1\frac{(-1)^{n+1}x^{2n+2}}{1+x^2}dx\\ &=\sum_{k=0}^{n}\frac{(-1)^k}{2k+1}+(-1)^{n+1}\int_0^1\frac{x^{2n+2}}{1+x^2}dx\\ \end{split}$ 考慮右式，第 n+1 項會等於零當 n 驅近於無限大 $0\leq\int_0^1\frac{x^{2n+2}}{1+x^2}dx\leq\int_0^1x^{2n+2} = \frac{1}{2n+3} \;\rightarrow 0 \text{ as } n \rightarrow \infty$ 因此，當 n 趨近於無限大時 $\frac{\pi}{4}=\sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}$ 得證 ![](https://i.imgur.com/cy8kgp8.png) * Error rate ![](https://i.imgur.com/3IsDRGK.png) ### Euler $公式：\frac{\pi^2}{6} = \sum_{k=1}^\infty\frac{1}{k^2}$ ![](https://i.imgur.com/lMLjQv0.png) * error rate ![](https://i.imgur.com/vuCQcnz.png) ## 其他計算方法 ### 阿基米德割圓法 閱讀：[圆周率是怎么计算的？祖冲之的缀术已经失传？李永乐老师5分钟带你了解割圆法](https://www.youtube.com/watch?v=lS1gnJxHPX4) 阿基米德利用正多邊形慢慢逼近出圓的周長，從正六邊形開始，一路逼近到正九十六邊型 若正 N 邊型的邊長為 $L_n$ 則正 2N 邊型的邊長為 $L_{2N}$ $L_{2N} = \sqrt{{2}-{\sqrt{{4}-{{L_N}^2}}}}$ 假設一圓半徑為一公分，則由半徑切割出來的正六邊形邊長 $L_6 = 1$ $L_{12} = \sqrt{{2}-{\sqrt{{4}-{{L_1}^2}}}} = \sqrt{{2}-{\sqrt{3}}}$ 依此可一直逼近到$6·2^k$邊型，當邊越多時，則誤差越小 ![](http://p1.pstatp.com/large/50ed000741bc5544cfee) :::danger 未完待續 :::