2016q3 Homework1 (compute-pi)

contributed by <`andy19950`>

執行 make check

/*----- baseline -----*/
16.23 user	0.00 system	0:16.23 elapsed	99% CPU

/*----- OpenMP 2 threads -----*/
16.21 user	0.00 system	0:08.11 elapsed	199% CPU

/*----- OpenMP 4 threads -----*/
17.70 user	0.00 system	0:04.75 elapsed	372% CPU

/*----- AVX SIMD -----*/
7.05 user	0.00 system	0:07.05 elapsed	99% CPU

/*----- AVX SIMD + loop unrolling-----*/
6.44 user	0.00 system	0:06.45 elapsed	99% CPU

可以發現使用 AVX SIMD 256bits暫存器來平行處理可以在執行時間上大幅縮短

使用 gnuplot 繪製執行結果

首先需要新增.gp檔，我從第一個作業那邊複製過來做修改。
執行 make gencsv 會產生 result_clock_gettime.csv 裡面有執行各個funcion的時間
因為 gnuplot 不能用逗號 "，" 來做分隔所以必須修改 benchmark_clock_gettime.c
Makefile也要稍做修改，之後就可以使用 using 每筆資料都把第一列當作 x軸，輸出當 y軸。
如此執行就可以畫出每筆資料的迭代次數(X) 以及執行時間(Y) 的比較圖。
logscale 可以對某個軸取 log ，我對 x 軸取 log10 讓圖看起來比較好看。







set logscale x 10

plot "result_clock_gettime.csv" using 1:2 with lines title 'baseline', \
'' using 1:3 with lines title 'omp_2', \
'' using 1:4 with lines title 'omp_4', \
'' using 1:5 with lines title 'avx', \
'' using 1:6 with lines title 'avxunroll'

實做 Leibniz formula

其實用上面的式子可以很簡單的實做出來











double compute_pi_Leibniz(size_t N)
{
    double pi = 0.0;
    double x;
    for (size_t n = 0; n < N; n++) {
        int temp = (n % 2) ? -1 : 1;
        x = (double) temp / (2 * n + 1);
        pi += x;
    }
    return pi * 4.0;
}

嘗試使用 AVX SIMD 來優化 Leibniz formula

AVX 256 bits 暫存器為 ymm0~ymm15。
AVX 指令為 _mm256_op_suffix
op 我們有使用到 set, set1, add, mul, div, store
suffix 是參數的型態，因為我們用doulbe來計算所以使用 packed double(pd)
因為使用 double 的關係 set 可以傳四個 doulbe (64bits) 參數把整個暫存器塞滿。
而 set1 則是傳一個 double 參數，暫存器的四個部份都儲存同一個參數。





















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);
    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(ymm1,ymm3);
        ymm3 = _mm256_add_pd(ymm3,ymm2);
        ymm3 = _mm256_div_pd(ymm0,ymm3);
        ymm4 = _mm256_add_pd(ymm3,ymm4);
    }
    double tmp[4] __attribute__((aligned(32)));
    _mm256_store_pd(tmp, ymm4);
    pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
    return pi * 4.0;
}

因為要實現 (-1)ⁿ 所以把 1 跟 -1 交錯存入ymm0
2*n + 1的部份用 ymm1, ymm2, ymm3來實做
ymm4則當作答案的暫存器，最後使用 store 把 ymm4 存入記憶體 tmp 內
再把四個計算好的部份全部加起來 *4 就是答案 π !!

再使用 make check 跑跑看

/*----- Leibniz -----*/
N = 1000000000 , pi = 3.1415926546
7.61 user	0.00 system	0:07.61 elapsed	99% CPU 

/*---- Leibniz + AVX SIMD -----*/
N = 1000000000 , pi = 3.1415926546
3.32 user	0.00 system	0:03.32 elapsed	100% CPU

結果證明萊布尼茲(Leibniz)的結論是對的。
而使用 AVX 來加速一樣有顯著的效果。

使用 gnuplot 把所有結果畫成圖表

其他估計 π 的方法

Nilakantha 級數 ––圖片來源：wikiHow

據說使用電腦運算分析後， Nilakantha 級數會比 Leibniz formula 算的更快更準！
有時間我會寫程式來印證看看他的說法。

結論

又把快遺忘的微積分拿出來複習一下，也發現很多方法可以來計算 π。
gnuplot 又更會使用了，圖也畫的比較好看了。
使用電腦來做數學計算並驗證公式真的很有趣，也是我從來沒有想過的電腦使用方式。
本來前一個作業要用懂的 AVX SIMD 在這個作業終於搞懂他了。

2016q3 Homework1 (compute-pi)

contributed by <andy19950>

執行 make check

使用 gnuplot 繪製執行結果

實做 Leibniz formula

嘗試使用 AVX SIMD 來優化 Leibniz formula

再使用 make check 跑跑看

使用 gnuplot 把所有結果畫成圖表

其他估計 π 的方法

Nilakantha 級數 ––圖片來源：wikiHow

結論

Reference

contributed by <`andy19950`>