Try   HackMD

2016q3 Homework1(compute_pi)

目標

  • 學習透過離散微積分求圓周率
    • Leibniz formula for π
    • 積分計算圓周率π
    • Function
  • 著手透過 SIMD 指令作效能最佳化
  • 作業說明

Makefile

每寫一份功課又是學習Makefile的時候=w=

  • 編譯時有-DBASELINE等等，其中的-D指的是#define BASELINE 1
    • 這樣程式碼就可以塞在同一個.c .h中！只要在編譯的時候加上-D就好了！！
    • Makefile的篇幅可以大幅降低！
    • 我想phonebook也可以這樣改寫，不然不同的方法就還要在多生一組.c .h

實驗

  • make check結果：
























time ./time_test_baseline
N = 400000000 , pi = 3.141593
7.59user 0.00system 0:07.59elapsed 100%CPU (0avgtext+0avgdata 1560maxresident)k
8inputs+0outputs (0major+82minor)pagefaults 0swaps

time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
8.03user 0.00system 0:04.01elapsed 200%CPU (0avgtext+0avgdata 1572maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps

time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
8.32user 0.00system 0:02.16elapsed 384%CPU (0avgtext+0avgdata 1676maxresident)k
0inputs+0outputs (0major+89minor)pagefaults 0swaps

time ./time_test_avx
N = 400000000 , pi = 3.141593
3.30user 0.00system 0:03.30elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps

time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
1.73user 0.00system 0:01.73elapsed 100%CPU (0avgtext+0avgdata 1608maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps
  • make gencsv ==> 輸出至表單
    輸出結果與在Makefile中看到的程式碼只差在$$i
    以下為輸出：




for i in `seq 100 5000 25000`; do \
		printf "%d," $i;\
		./benchmark_clock_gettime $i; \
	done > result_clock_gettime.csv

原來除了寫進txt檔之外，還可以寫進libre office的表單裡，好酷！

  • 理解程式碼：
    • 上面的程式碼相當於：





for(i=100; i<25000; i+=5000)
{
	printf("%d",i);
	benchmake_clock_gettime(i);
}

然後把結果輸出至.csv檔

細節的語法待會查 先highlight

  • 在Makefile中，在執行各檔時前面多加了time

    • time: 查看指令的執行時間

  • 直接執行time ./time_test_baseline

N = 400000000 , pi = 3.141593

real	0m7.805s
user	0m7.808s
sys	0m0.000s

  • 顯示了三種時間：real,user和sys
    Real refers to actual elapsed time; User and Sys refer to CPU time used only by the process.

    • real: real time,a wall clock time,指令從開始到結束執行的時間，含其他process
    • user: user CPU time,指令在user mode的時間,其他process的時間不算
    • sys: system CPU time,指令在kernel mode的時間,其他process的時間不算

  • real = user + sys

    • 並非如此，如sleep(),不論在user mode或kernel mode沒佔用CPU時間，那就不會被算入，反觀real time仍會把sleep()的時間算進去

Baseline


為上述積分式的實作，最後回傳的值要*4，因為算出來是
π4

  • N愈大[0,1]就被切割成更多等份，dt也就愈小










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

  • 也是用跟baseline一樣的積分式求得π
  • 可以傳遞參數決定要開幾個threads
  • 在time_test.c中分別傳入2個和4個threads的實驗
  • 程式碼：
    • private(x)：
      讓每個thread都有自己的x,不然threads之間會把x當成共用的變數
    • reduction(+:pi)：讓各個thread針對指定的變數擁有一份有起始值的複本
      • 支援的運算元有 +, *, –, &, ^, |, &&, ||；而變數則必須要是 shared 的















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 SIMD

  • 甚麼是AVX？

    • AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集
    • AVX指令集借鑑了一些AMD SSE5的設計思路，進行擴展和加強，形成一套新一代的完整SIMD指令集規範。
    • SSE將向量處理能力從64位提升到128位後，SSE系列指令都只能使用128位XMM暫存器，這次AVX將所有16個128位XMM暫存器擴充為256位的YMM暫存器，從而支持256位的向量計算。

  • prototype:

    • _mm256 : single precision,32-bit
    • _mm256d : double precision, 64bit

配上註解與下方的連結終於稍微看懂了~

  • __attribute__ 機制，可以用來設置函數屬性（Function Attribute）、變數屬性（Variable Attribute）和類型屬性（Type Attribute）。
  • aligned 則規定變數或結構的最小對齊格式，以 Byte 為單位。

不知道為甚麼沒有highlight@@
























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; 
}
  • 程式碼解讀(可能有誤！！)
    • __m256d:存放double,256-bit可用，因此可以存放256/64 = 4 個double
    • 每次的運算是4個double 個別 加總
    • 最後使用tmp[4]把4個值存起來，最後在加總
      • (aligned(32)) 是因為double有8個byte,又宣告 double tmp[4] ==> 8*4= 32(byte)

Leibniz

  • 程式：










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

一直出現錯誤資訊：
應該是link的問題，但是即使加了-lm還是沒用@@

computepi.c:(.text+0xd80): 未定義參考到 pow
collect2: error: ld returned 1 exit status
make: *** [default] Error 1

computepi.c 有included <math.h>嗎？ 或許可以試試，我花很多時間處理這個…SarahYuHanCheng

  • 只好先改成



int tmp = (i%2) ? -1 : 1;    
x = (double)tmp / (2*i +1) ; //記得轉型
pi += x;

Leibniz AVX

看了一些資料在搭配原有的程式碼，很快就可以寫出來了~ 覺得開心=w=























double compute_pi_leibizavx(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); //for pow(-1,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); //i i+1 i+2 i+3
		ymm3 = _mm256_mul_pd(ymm1, ymm3);	//2*i 2*(i+1)...
		ymm3 = _mm256_add_pd(ymm3,ymm2);	//2*i+1 2*(i+1)+1 ...
		ymm3 = _mm256_div_pd(ymm0,ymm3);	//(-1)^i/(2*i+1) ...
		ymm4 = _mm256_add_pd(ymm3,ymm4);	//sum
		 
    }
    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;
}

其他計算圓周率的方法

  • 蒙地卡羅法：
    • 邊長為2的正方形，正方形內有一個內切圓，面積為
      1⋅1⋅π=π
    • 圓面積和正方形面積之比為π:4 = 正方形內均勻隨機丟石頭，落在圓形內的機率













double compute_pi_MonteCarlo(size_t N)
{
     int count = 0;
    srand(time(NULL));
	
    for (size_t i = 0; i < N; i++) {
        double x = (double)rand()/RAND_MAX;      
		double y = (double)rand()/RAND_MAX;
        if((x*x+y*y)<=1) count++;
    }

    return 4.0 * count / N ;
}

Plot

  • 複習語法
    • using [a：b] : 以第a column當x座標,第b column當y座標

製圖時產生的warning:

warning: Skipping data file with no valid points

查了一段時間，玄機就在副檔名.csv和gnuplot他吃數據的方式啊！！

csv : comma separated value
所以資料之間是以 逗號 區隔，因此要在gnuplot script裡加一行set datafile separator ","才讀的到資料

  • 不同方法比較圖

抖的超厲害！

取樣從N=1000~N=250000，間隔1000

隔一天跑出來的結果差好多@@
omp_4沒那麼劇烈了，還不知道為甚麼會這樣

新增Leibniz

  • 可以發現Leibniz與omp_2的時間差不多，但又比omp_2還穩定
  • LeibnizAVX也比avx的時間再更短
  • LeibnizAVX_unroll語omp_4幾乎重合

思考如何精準計算時間

  • walk clock time: 由xtime記錄，系統每次啟動時會先從設備上的 RTC 上讀入 xtime
    • xtime:從cmos電路中取得的時間，一般是從某一歷史時刻開始到現在的時間
    • Jeffie:紀錄系統開幾以來，經過多少個 tick，取決於系統的頻率，單位是Hz
  • system time:
    • 系統上的時間，開機時它會讀取RTC 來設定，可能過程中還會有時區換算等之類的設置。

信賴區間

參考資料

TempoJiJi hackmd

tags: system embedded HW1-3
Last changed by 
 
Vivian Lin
0
676
Published on HackMD