2017 q1 Homework1 (compute-pi)

contributed by <zmke>

Reviewed by <`rayleigh0407`>

建議可以證明一下 Leibniz formula

Reviewed by <`twzjwang`>

執行時同時使用其他程式，為什麼 omp_4 受到的影響最大?
可以探討作業要求中時間函式部分提到的問題、延伸思考

開發環境

OS: Ubuntu 16.04 LTS
Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
Model name: Intel® Core™ i5-4200H CPU @ 2.80GHz
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 3072K

AVX

Advanced Vector Extensions
Intel 開發的 Single Instruction Multiple Data 指令集，在一個 cycle 內可以同時對多個數據執行同樣的指令
ymm0~ymm15 共 16 個暫存器，每個暫存器內可以放256 bits
- e.g. 8個 32-bit float; 4個 64-bit double

時間處理與 time 函式使用

NTP (Network Time Protocol)
- 透過網路封包交換，將電腦之間的時鐘校正和同步
clock()
- return CPU time
- 回傳值是 clock_t type ，要除以 CLOCKS_PER_SEC 才能取得以秒為單位的時間
- On a 32-bit system where CLOCKS_PER_SEC equals 1000000 this function will return the same value approximately every 72 minutes. （超過72分鐘會 wrap around ）
- 在 muti threads 的情況下， CPU time 是每條 thread 的加總
clock_gettime()
- int clock_gettime(clockid_t clk_id, struct timespec *tp);
- retrieve and set 的類型由 clk_id 決定
  - CLOCK_REALTIME : system-wide clok ，會回傳 real time 也就是 wall-clock ，如果使用者任意更動系統時間，會影響測量結果
  - CLOCK_MONOTONIC : 不受 discontinuous jumps 影響
- 精度為 nanoseconds
time()
- wall-clock
- return 1970-01-01 00:00:00 +0000 (UTC) 到目前為止的時間，精度為 seconds

以前寫作業的時候曾經遇過 multi thread 時間測量的問題，一開始是用 clock() ，發現測出來的時間跟預期相差很大，現在終於知道原因了

原始版本

修改 Makefile








gencsv: default
	for i in `seq 500 500 250000`; do \
		printf "%d," $$i;\
		./benchmark_clock_gettime $$i; \
	done > result_clock_gettime.csv	
    
plot: gencsv
        gnuplot script.gp

從 N=500 開始，間隔 500 測試到 N=250000
$ make plot

起伏有點誇張，尤其是 omp_4，發現可能的原因是我邊等輸出邊用 youtube 聽音樂，把網頁關起來再試一次

這次在等輸出的時候沒有對電腦進行任何操作，明顯比上一次穩定很多

Leibniz






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

omp 版本










double compute_pi_leibniz_openmp(size_t N, int threads){
	double pi = 0.0;
	#pragma omp parallel num_threads(threads)
    {
        #pragma omp for reduction(+:pi)
        for(size_t i = 0; i < N; i++)
			pi += pow(-1.0, i) / (2*i + 1);
    }
    return pi * 4.0;
}

$ make plot

發現 Leibniz 求圓周率明顯比其他方法慢，猜測是 pow(-1.0, i) 太浪費時間，修改版本如下







double compute_pi_leibniz(size_t N){
	double pi = 0.0;
	for(size_t i = 0; i < N; i++){		
		pi += (i % 2) ? (-1.0 / (2*i + 1)) : (1.0 / (2*i + 1));
	}
	return pi * 4.0;
}

omp










double compute_pi_leibniz_openmp(size_t N, int threads){
	double pi = 0.0;
	#pragma omp parallel num_threads(threads)
    {
        #pragma omp for reduction(+:pi)
        for(size_t i = 0; i < N; i++)
			pi += (i % 2) ? (-1.0 / (2*i + 1)) : (1.0 / (2*i + 1));
    }
    return pi * 4.0;
}

模仿原本的 avx，練習把 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);		//pow(-1.0, 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);
		ymm3 = _mm256_mul_pd(ymm3, ymm1);			//2*i
		ymm3 = _mm256_add_pd(ymm3, ymm2);			//2*i+1
		ymm3 = _mm256_div_pd(ymm0, ymm3);			//pow(-1.0, i) / (2*i+1)
		ymm4 = _mm256_add_pd(ymm4, ymm3);			//pi += pow(-1.0, i) / (2*i+1)
	}
	double tmp[4] __attribute__((aligned(32)));
	_mm256_store_pd(tmp, ymm4);
	pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
	return pi * 4.0;
}

$ make plot

發現把 pow(-1.0, i) 改掉之後， leibinz_omp4 變成執行速度最快的版本

Compute Error

以 arccos(-1) 為基準計算誤差









    // Baseline
	result = compute_pi_baseline(N);
	error = (result - PI > 0) ? result - PI : PI - result;
	printf("%lf,", error);

	// Leibniz
	result = compute_pi_leibniz(N);
	error = (result - PI > 0) ? result - PI : PI - result;
	printf("%lf\n", error);

$ make plot_error

兩條線幾乎重疊，表示在同樣的 N 精度差距不大

把判斷正負的判斷式移除後再畫一次，發現兩個求圓周率方法分別是從上下逼近
$ make plot_error

Reference

AVX是什麼？AVX指令集技術與應用解析
 TempoJiJi 筆記
 clock()
clock_gettime()
NTP
Leibniz formula for pi