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

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 →

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

\frac{π}{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 \cdot 1 \cdot π = π$
- 圓面積和正方形面積之比為π: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

Time
OpenMP
- openmp_private
- openmp_reduction
AVX
- introduction
- manual(chap11)
- functions