Try   HackMD

2017q1 Homework3 (compute-pi)

contributed by <xdennisx>

Reviewed by baimao8437

  • 每次新增 method 都有新的分析圖形,但不能只丟圖上來就好,應稍微解釋圖片中的現象,如: Nilakantha vs Leibniz & Euler 的 error rate 圖形看起來一模一樣,沒有其他文字解釋的話看不出差異。
  • 對於資料的抖動的情況,可以嘗試使用95%信賴區間對資料進行篩選。
  • 由於此作業對於計算
    π
    會使用多種 method ,使得每新增一種,就得在 makefile 中多寫很多行,可嘗試簡化 makefile 已提升測試效率,可參考神之修改

原始圖形

將於本 result_clock_gettime.csv 的圖形用 gunplot 輸出

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

OpneMP(4) 一開始都比大家高?

error rate

baseline

原始公式:

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; }

OpenMP

指定要用幾個thread來執行

reduction

支援的運算元有 +, , –, &, ^, |, &&, ||,它讓各個執行緒針對指定的變數擁有一份有起始值的複本(起始值是運算元而定,像 +, – 的話就是 0, 就是 1),平行化計算時,以各自複本做運算,等到最後再以指定的運算元,將各執行緒的複本合在一起。

AVX SIMD

AVX 是一組 SIMD(Single Instruction Multiple Data)的指令集,擴展了原有的 MMX 及 Intel Streaming SIMD Extensions(Intel SSE)

支援 AVX 的硬體上會有 16 個 256-bit 的 YMM(YMM0-YMM15) 暫存器,還有一個 32-bit 的控制暫存器 MXCSR,舊有的 XMM 暫存器可以視為取 YMM 的 lower half,如下圖。

  • _mm256d : 它並不是一種暫存器,是指可以用來載入到 AVX 暫存器的 “Data type”,double precision, 64bit

  • _mm256_set1_pd(1.0):將參數浮點數值放到 _mm256 變數的所有位置。

  • _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0):將dt * 3, dt * 2, dt * 1, 0.0這些參數依序放入_mm256 變數,參數順序和放進去的次序相反。

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; }

時間處理與 time 函式使用

clock()

clock_gettime()

Wall-clock time vs. CPU time

Wall-clock time

顧名思義就是掛鐘時間,也就是現實世界中實際經過的時間,是由 kernel 裡的 xtime 來紀錄,系統每次啟動時會先從設備上的 RTC 上讀入 xtime。

struct timespec xtime __attribute__ ((aligned (16)));
struct timespec {
    __kernel_time_t tv_sec; /* seconds */
    long tv_nsec;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; /* nanoseconds */
};

這個值是自 1970-01-01 起經歷的秒數,另外 Linux 核心也會決定每秒要發生幾次中斷 (透過常數 HZ) ,每次 timer interrupt,就會去更新 xtime

另外注意,有時系統要與網絡中某個節點時間同步,那這個 Wall-clock time 可能會和 NTP 伺服器時區日光節約時間同步或使用者自己調整。所以「通常」不建議拿來量測程式時間,因為它不是一個穩定的時間,用專業一點的用語講,Wall-clock time 不一定是單調遞增 (monotonic)。

CPU time

指的是程序在 CPU 上面運行消耗 (佔用) 的時間,clock() 就是很典型用來計算 CPU time 的時間函式,但要注意的是,如果有多條 threads,那 CPU time 算的是「每條 thread 的使用時間加總」,所以如果我們發現 CPU time 竟然比 Wall-clock time 還大!這可能是個很重要的原因。

另外很多時候只計算 CPU time 是不夠的,因為執行時間可能還包括 I/O time、 communication channel delay、synchronization overhead等等。

  • real time : 也就是 Wall-clock time,當然若有其他程式同時在執行,必定會影響到。
  • user time : 表示程式在 user mode 佔用所有的 CPU time 總和。
  • sys time : 表示程式在 kernel mode 佔用所有的 CPU time 總和 ( 直接或間接的系統呼叫 )。

其他計算 pi 的方式

Leibniz formula for π

考慮以下分解:

image alt

對兩邊從0到1去做積分

時,除積分項以外的項收斂到萊布尼茨級數。同時,積分項收斂到0:

所以這便證明了萊布尼茨公式。

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

Euler


通式
π26
=
i=1n1n2

double compute_pi_euler(size_t N) { double pi = 0.0; for (size_t i = 1; i <= N; i++) { pi += (1 / pow(i, 2.0)); } pi *= 6; return sqrt(pi); }


Nilakantha 級數

π = 3 + 4/(234) - 4/(456) + 4/(678) - 4/(8910) + 4/(101112) - (4/(121314)

通式大概長這樣
π = 3+

i=2n4(1)n(2n)(2n1)(2n2)

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

某個資料說這個比萊布尼茨快,結果好像沒有,可能是環境或是其他因素


Reference