2017q1 Homework1(compute-pi)

contributed by < baimao8437 >

Reviewed by stanleytazi

• 信賴區間的公式目前看到的都是 standard error 下去算
  ，而非standard deviation
  confidence interval
• 從比較圖來看好像都有看到OpenMP with 4 threads有比較大的抖動，有點好奇為什麼
• default 設定跑25次來算，可以再試試看跑多次一點，看不同做法的時間會不會差的比較明顯

Reviewed by hunng

• 利用多種方式計算並使用 AVX 實作
• 中間提到的不確定有沒有信賴區間，可以多次執行儲存結果後分析看看
  • 而且有加上信賴區間比沒有加的抖動還要更大是不是出了什麼問題

花費時間 _{2/28 - 3/2}

• 讀書：9 hr
• 文件：4 hr
• coding: 3 hr
• 看著發呆：3 hr

目標設定

• 習慣 linux 系統基本操作
• 熟練 git 基本版本控管操作
• 繪圖工具 LibreOffice、gnuplot
• SIMD 實作練習
• 複習微積分？
• Hackmd 文件 LaTeX 數學式寫法
• 減少發呆時間…

前置作業

這個作業一開始我整個霧煞煞 沒了老溼的精美影片教學 頓失依靠
只好先讀神人們的共筆了解一些基本名詞與時間分析

開發環境

baimao@baimao-Aspire-V5-573G:~$ lscpu
Architecture:          x86_64
CPU 作業模式：    32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                4
On-line CPU(s) list:   0-3
每核心執行緒數：2
每通訊端核心數：2
Socket(s):             1
NUMA 節點：         1
供應商識別號：  GenuineIntel
CPU 家族：          6
型號：              69
Model name:            Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
製程：              1
CPU MHz：             1711.242
CPU max MHz:           2600.0000
CPU min MHz:           800.0000
BogoMIPS:              4589.38
虛擬：              VT-x
L1d 快取：          32K
L1i 快取：          32K
L2 快取：           256K
L3 快取：           3072K
NUMA node0 CPU(s):     0-3

SIMD指令集

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

效能分析

測試原版程式效能：

time ./time_test_baseline
N = 400000000 , pi = 3.141593
4.40user 0.00system 0:04.40elapsed 99%CPU (0avgtext+0avgdata 1808maxresident)k
0inputs+0outputs (0major+80minor)pagefaults 0swaps

time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
4.91user 0.00system 0:02.46elapsed 199%CPU (0avgtext+0avgdata 1772maxresident)k
0inputs+0outputs (0major+80minor)pagefaults 0swaps

time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
9.76user 0.00system 0:02.45elapsed 397%CPU (0avgtext+0avgdata 1792maxresident)k
0inputs+0outputs (0major+85minor)pagefaults 0swaps

time ./time_test_avx
N = 400000000 , pi = 3.141593
1.26user 0.00system 0:01.27elapsed 99%CPU (0avgtext+0avgdata 1816maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps

time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
1.22user 0.00system 0:01.22elapsed 99%CPU (0avgtext+0avgdata 1816maxresident)k

可以看到 user(user mode CPU time) 跟 system (kernel mode CPU time) 時間
也可以看到佔用了多少 CPU OpenMP 2 threads 最多佔用 200% CPU
OpenMP 4 threads 最多就佔用 400% CPU

由於我還沒做好算數學的心理準備 就先來畫畫好了~

想要建立出作業範例以 LibreOffice 建立的圖表（如下圖）：

先到 Makefile 的 gencsv 看到 seq 100 5000 25000
等於是 for(i=100; i<25000; i+=5000)
我調成 seq 1000 200 400000 跑跑看將近兩千筆資料
然後把出來的 csv 以 LibreOffice Calc 開啟 選擇以逗號分隔
再把全部的資料選取好後 按下圖表 再做一些設定：

1. 圖表類型：設定為 線條圖 僅限線條 類型：平滑
2. 資料範圍：勾選 已欄表示的資料序列 + 第一欄作為標籤
3. 資料序列：把每條線名子改為對應 method
4. 圖表元素：x 軸為 N, y 軸為 Time(sec)

這份報告之後可以拿來做研究所推甄、工作實習，甚至可作為論文的雛形，用與請精確 –jserv

Yes, Sir!baimao8437

畫出來後:

預料之中 跟原圖差距甚遠
連一點指數上升的感覺都沒有 無法理解
而且原圖是 Baseline > OpenMP_2 > AVX > OpenMP_4 > AVX_unrolling
我的結果是 Baseline > OpenMP_4 > OpenMP_2 > AVX > AVX_unrolling
好像差了有點多 不過參考了別人的圖發現有人跟我一樣結果 就安心多了

補上 gnuplot perfomance + error rate:

信賴區間

可看出圖表顯示有一些極端值造成線條的抖動
因此參考學長的筆記 試圖建立95信賴區間來排除極端值

機統沒學好 信賴區間忘光光： 簡單來說就這三句話

計算平均值
計算標準差(standard deviation)
平均值的正負標準差乘以 2 就是 95% 信賴區間

標準差長這樣:

$$s=\sqrt{\frac{1}{n-1}\sum _{i=1}^n(x_i-\overline{x}{)}^2}$$
然後 makefile 中加入 -lm 一樣讓他跑 2000 筆資料
結果:

OpenMP 4 threads 還是抖了一下到讓我懷疑到底有沒有信賴區間…

補上 gnuplot:

可以看到 一樣的演算法 但是用不同的優化結果
AVX_unrolling 的優化會造成再某些 N 計算結果誤差比較大
其他的算法 error rate 皆相同

Leibniz 公式證明

證明：

$$\frac{\pi }{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots$$
首先考慮

$$1-x^2+x^4-x^6+x^8-\cdots =\frac{1}{1+x^2},|x|<1$$
積分後

$$x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}-\cdots =ta{n}^{-1}x,|x|<1$$
接著證明收斂

$$\frac{1}{1+x^2}=1-x^2+x^4-x^6+x^8-\cdots +(-1)x^{2n}+\frac{(-1{)}^{n+1}{x}^{2n+2}}{1+x^2}$$
兩邊積分後等於

$$\frac{\pi }{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots +\frac{-1^n}{2n+1}+(-1{)}^{n+1}{\int }_0^1\frac{x^{2n+2}}{1+x^2}dx$$
當

$n\to \mathrm{\infty }$時, 除積分項以外的項收斂到萊布尼茲級數, 同時積分項收斂至0

$${\int }_0^1\frac{x^{2n+2}}{1+x^2}dx<{\int }_0^1{x}^{2n+2}dx=\frac{1}{2n+3}\to 0,whenn\to \mathrm{\infty }$$
得證!

$$\pi =4\sum _{n=1}^{\mathrm{\infty }}\frac{(-1{)}^n}{2n+1}$$
接著寫成 code:

double compute_pi_leibniz(size_t N)
{
    double pi = 0.0;
    int tmp = 1;
    for (size_t i = 0; i < N; i++) {
        pi += ((double)tmp / (2.0 * (double)i + 1.0));
        tmp *= -1;
    }
    return pi * 4.0;
}

第一次寫用了 pow(-1,i)去變換正負
但太花費時間 所以用 tmp 去變換正負
然後 Leibniz 優化挑了 AVX 來實作:

double compute_pi_leibniz_avx(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);//represent tmp
    ymm1 = _mm256_set1_pd(1);
    ymm2 = _mm256_set1_pd(2);
    ymm4 = _mm256_setzero_pd();             // sum of pi

    for (int i = 0; i <= N - 4; i += 4) {
        ymm3 = _mm256_set_pd((double)i, (double)i + 1.0, (double)i + 2.0, (double)i + 3.0); // represent n
        ymm3 = _mm256_mul_pd(ymm2, ymm3);   // [2*n]
        ymm3 = _mm256_add_pd(ymm1, ymm3);   // 2*n[+1]
        ymm3 = _mm256_div_pd(ymm0, ymm3);   // [tmp/](2*n+1)
        ymm4 = _mm256_add_pd(ymm4, ymm3);   // pi += tmp/(2*n+1)
    }
    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 但有了上面的範例就好懂多了 也比較容易寫出來

• gnuplot 圖表
  
  有點緊密 不過可以看到 AVX_unrolling 還是比 Leibniz_avx 快一點點
  
  但是 error rate 卻明顯 Leibniz_avx 優於 AVX_unrolling
  從輸出 csv 來看這邊 因為 Leibniz & Leibniz_avx 的 error rate 又全部跟其他重疊了 所以只有一條線

其他算

$\pi$ 的方法 : Nilakantha method

$$\pi =3+\frac{4}{2\times 3\times 4}-\frac{4}{4\times 5\times 6}+\frac{4}{6\times 7\times 8}-\frac{4}{8\times 9\times 10}+\cdots$$
寫成 code 也不難 大概長這樣

double compute_pi_nilakantha(size_t N)
{
    double pi = 0.0;
    int tmp = 1;
    double start = 0.0;
    for (size_t i = 1; i < N; i++) {
        start = (double)i * 2.0;
        pi += (double)tmp / ((start) * (start + 1.0) * (start + 2.0));
        tmp *= -1;
    }
    pi *= 4;
    pi += 3;
    return pi;
}

但是如果要用 AVX 優化的話 就要稍微細心一點了 改到頭昏腦脹

double compute_pi_nilakantha_avx(size_t N)
{
    double pi = 0.0;
    register __m256d ymm0, ymm1, ymm2, ymm3, ymm4, ymm5, ymm6;
    ymm0 = _mm256_set_pd(1.0, -1.0, 1.0, -1.0);
    ymm1 = _mm256_set1_pd(1);
    ymm2 = _mm256_set1_pd(2);
    ymm3 = _mm256_set1_pd(4);
    ymm4 = _mm256_setzero_pd();//_mm256_set1_pd(4);

    for (int i = 0; i < N - 4; i += 4) {
        ymm6 = _mm256_set_pd((double)i + 1.0, (double)i + 2.0, (double)i + 3.0, (double)i + 4.0); //n
        ymm6 = _mm256_mul_pd(ymm2, ymm6); //start = i*2
        ymm5 = _mm256_add_pd(ymm6, ymm1); //(start+1)
        ymm6 = _mm256_mul_pd(ymm6, ymm5); //start*(start+1)
        ymm5 = _mm256_add_pd(ymm5, ymm1); //(start+2)
        ymm6 = _mm256_mul_pd(ymm6, ymm5); //start*(start+1)*(start+2)
        ymm6 = _mm256_div_pd(ymm3, ymm6); // 4/start*(start+1)*(start+2)
        ymm6 = _mm256_mul_pd(ymm6, ymm0);
        ymm4 = _mm256_add_pd(ymm4, ymm6);
    }
    double tmp[4] __attribute__((aligned(32)));
    _mm256_store_pd(tmp, ymm4);
    pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
    pi += 3.0;
    return pi;
}

我相信應該有更簡潔的寫法 但是至少這個結果是對的
來看圖表

效能來看最好的三個是 :
(時間少) AVX > AVX_unrolling > Leibniz_avx (時間多)

這個圖表看起來跟上面的 error 圖看起來一樣
因為 Nilakantha 的兩個算法 錯誤率都非常接近 0
錯誤率低的嚇人 我把 N 調小 單獨看 Nilakantha

 N,error rate       N,error rate
 5,0.000608        55,0.000000
10,0.000079        60,0.000000
15,0.000023        65,0.000000
20,0.000010        70,0.000000
25,0.000005        75,0.000000
30,0.000003        80,0.000000
35,0.000002        85,0.000000
40,0.000001        90,0.000000
45,0.000001        95,0.000000
50,0.000001       100,0.000000

只要算到 50 幾項 就幾乎等於

$\pi$ 了

參考資料：

• HankLo 的筆記
• HaoTse 的筆記
• 廖建富學長筆記
• 萊布尼茲公式
• Nilakantha 公式
• 勃興大神的 makefile