2019q3 Homework4 (compute-pi)

contributed by < xl86305955 >

tags:sysprog2019 Homework4 compute-pi

目標

• 複習數值系統
• 觀察和使用 SIMD 指令
• 練習透過離散微積分求圓周率
  • Leibniz formula for π
  • 積分計算圓周率π
• 著手透過 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 之中

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

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

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

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

• Error rate
  
  

LEIBNIZ

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

proof

$tan\theta =1$

$\theta =\frac{\pi }{4}=arctan(1)$

其中

$\mathrm{\forall }x\in [-1,1]arctan(x)=\sum _{k=0}^{\mathrm{\infty }}(-1{)}^k\frac{x^{2k+1}}{2k+1}=x-\frac{1}{3}{x}^3+\frac{1}{5}{x}^5-\frac{1}{7}{x}^7+\cdots$

$arcta{n}^{\prime }(x)=\frac{1}{1+x^2}$

$\begin{array}{r}\\ arctan(1)& ={\int }_0^1\frac{1}{1+x^2}dx\\ & ={\int }_0^1(\sum _{k=0}^n(-1{)}^k{x}^{2k}+\frac{(-1{)}^{n+1}{x}^{2n+2}}{1+x^2})dx\\ & ={\int }_0^1\sum _{k=0}^n(-1{)}^k{x}^{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{array}$

考慮右式，第 n+1 項會等於零當 n 驅近於無限大

$0\le {\int }_0^1\frac{x^{2n+2}}{1+x^2}dx\le {\int }_0^1{x}^{2n+2}=\frac{1}{2n+3}\to 0\text{ as }n\to \mathrm{\infty }$

因此，當 n 趨近於無限大時

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

得證

• Error rate
  
  

Euler

$\mathrm{公}\mathrm{式}：\frac{{\pi }^2}{6}=\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^2}$

• error rate
  
  

其他計算方法

阿基米德割圓法

閱讀：圆周率是怎么计算的？祖冲之的缀术已经失传？李永乐老师5分钟带你了解割圆法

阿基米德利用正多邊形慢慢逼近出圓的周長，從正六邊形開始，一路逼近到正九十六邊型

若正 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$邊型，當邊越多時，則誤差越小

未完待續