OpenMP并行實(shí)驗(yàn)報(bào)告

1、本文格式為Word版，下載可任意編輯OpenMP并行實(shí)驗(yàn)報(bào)告 并行試驗(yàn)報(bào)告 一、 積分計(jì)算圓周率 1.1 積分計(jì)算圓周率的向量優(yōu)化 1.1.1 串行版本的設(shè)計(jì) 任務(wù)：理解積分求圓周率的方法，將其用 C 代碼實(shí)現(xiàn)。 留意：理論上，dx 越小，求得的圓周率越精確；在計(jì)算機(jī)中由于表示的數(shù)據(jù)是有精度范圍的，假如 dx 太小，積分次數(shù)過多，誤差積累導(dǎo)致結(jié)果不精確。 以下為串行代碼： #includestdio.h #includetime.h #define N 10000000 double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int

2、 i; for(i=0; idt; i+) double x=(double)i/dt; pi+=delta/(1.0+x*x); return pi*4; int main() int dx; double pai; double start,finish; dx=N; start=clock(); pai=get_pi(dx); finish=clock(); printf(%.8lfn,pai); printf(%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); return 0; 時(shí)間運(yùn)行如下： 第一次：time=0.02674000S 其次

3、次：time=0.02446500S 第三次：time=0.02402800S 三次平均為：0.02508S 1.1.2 SSE 向量優(yōu)化版本設(shè)計(jì) 任務(wù)：此部分需要給出單精度和雙精度兩個(gè)優(yōu)化版本。 留意： （1）測(cè)試均在劃分度為 10 的 7 次方下完成。 以下是 SSE 雙精度的代碼： #includestdio.h #includex86intrin.h #includetime.h #define N 10000000 double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i; for(i=0; idt; i+) d

4、ouble x=(double)i/dt; pi+=delta/(1.0+x*x); return pi*4; double get_pi_sse(size_t dt) double pi=0.0; double delta =1.0/dt; _m128d xmm0,xmm1,xmm2,xmm3,xmm4; xmm0=_mm_set1_pd(1.0); xmm1=_mm_set1_pd(delta); xmm2=_mm_set_pd(delta,0.0); xmm4=_mm_setzero_pd(); for(long int i=0; i=dt-2; i+=2) xmm3= _mm_set1

5、_pd(double)i*delta); xmm3= _mm_add_pd(xmm3,xmm2); xmm3= _mm_mul_pd(xmm3,xmm3); xmm3= _mm_add_pd(xmm0,xmm3); xmm3= _mm_div_pd(xmm1,xmm3); xmm4= _mm_add_pd(xmm4,xmm3); double tmp2 _attribute_(aligned(16); _mm_store_pd(tmp,xmm4); pi+=tmp0+tmp1/*+tmp2+tmp3*/; return pi*4.0; int main() int dx; double pai

6、; double start,finish; dx=N; start=clock(); pai=get_pi_sse(dx); finish=clock(); printf(%.8lfn,pai); printf(%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); return 0; 時(shí)間運(yùn)行如下： 第一次：time=0.00837500S 其次次：time=0.00741100S 第三次：time=0.00772021S 三次平均為：0.00783S 以下是 SSE 單 精度的代碼： #includestdio.h #includex86intri

7、n.h #includetime.h #define N 10000000 float get_pi_sse(size_t dt) float pi=0.0; float delta =1.0/dt; _m128 xmm0,xmm1,xmm2,xmm3,xmm4; xmm0=_mm_set1_ps(1.0); xmm1=_mm_set1_ps(delta); xmm2=_mm_set_ps(delta*3,delta*2,delta,0.0); xmm4=_mm_setzero_ps(); for(long int i=0; i=dt-4; i+=4) xmm3= _mm_set1_ps(fl

8、oat)i*delta); xmm3= _mm_add_ps(xmm3,xmm2); xmm3= _mm_mul_ps(xmm3,xmm3); xmm3= _mm_add_ps(xmm0,xmm3); xmm3= _mm_div_ps(xmm1,xmm3); xmm4= _mm_add_ps(xmm4,xmm3); float tmp4 _attribute_(aligned(16); _mm_store_ps(tmp,xmm4); pi+=tmp0+tmp1+tmp2+tmp3; return pi*4.0; int main() int dx; float pai; double star

9、t,finish; dx=N; start=clock(); pai=get_pi_sse(dx); finish=clock(); printf(%.8fn,pai); printf(%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); return 0; 時(shí)間運(yùn)行如下： 第一次：time=0.00406100S 其次次：time=0.00426400S 第三次：time=0.00437600S 三次平均為：0.00423S 1.1.3 A 向量優(yōu)化版本設(shè)計(jì) 任務(wù)：此部分需要給出單精度和雙精度兩個(gè)優(yōu)化版本 留意： （1）測(cè)試均在劃分度為 10 的 7

10、 次方下完成。 （2）在編譯時(shí)需要加-ma 編譯選項(xiàng)，才能啟用 A 指令集，否則默認(rèn) SSE 指令集 （3）理論上，向量版本對(duì)比 SSE 版本和串行版本有明顯加速，單精度版本速度明顯優(yōu)于雙精度，速度接近雙精度的兩倍。 以下是 A 雙精度的代碼： #includestdio.h #includex86intrin.h #includetime.h #define N 10000000 /*double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i; for(i=0; idt; i+) double x=(double)i/dt

11、; pi+=delta/(1.0+x*x); return pi*4; */ double get_pi_a(size_t dt) double pi=0.0; double delta =1.0/dt; _m256d ymm0,ymm1,ymm2,ymm3,ymm4; ymm0=_mm256_set1_pd(1.0); ymm1=_mm256_set1_pd(delta); ymm2=_mm256_set_pd(delta*3,delta*2,delta,0.0); ymm4=_mm256_setzero_pd(); for(long int i=0; i=dt-4; i+=4) ymm3=

12、 _mm256_set1_pd(double)i*delta); ymm3= _mm256_add_pd(ymm3,ymm2); ymm3= _mm256_mul_pd(ymm3,ymm3); ymm3= _mm256_add_pd(ymm0,ymm3); ymm3= _mm256_div_pd(ymm1,ymm3); ymm4= _mm256_add_pd(ymm4,ymm3); double tmp4 _attribute_(aligned(32); _mm256_store_pd(tmp,ymm4); pi+=tmp0+tmp1+tmp2+tmp3; return pi*4.0; int

13、 main() int dx; double pai; double start,finish; dx=N; start=clock(); pai=get_pi_a(dx); finish=clock(); printf(%.8lfn,pai); printf(%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); return 0; 時(shí)間運(yùn)行如下： 第一次：time=0.00720210S 其次次：time=0.00659800S 第三次：time=0.00670600S 三次平均為：0.00683S 以下是 A 單 精度的代碼： 時(shí)間運(yùn)行如下： 第一

14、次：time=0.00234200S 其次次：time=0.00234200S 第三次：time=0.00230000S 三次平均為：0.002328S 由以上試驗(yàn)統(tǒng)計(jì)得出結(jié)論： A-float=0.002328 S A-double=0.00683 S SSE-float=0.00423 S SSE-double= 0.00783 S 基本符合規(guī)律：（以下為速度比較） A-float A-double SSE-float SSE-double serial 1.2 積分計(jì)算圓周率的 OpenMP 優(yōu)化 1.2.1 OpenMP 并行化 任務(wù)：在串行代碼的基礎(chǔ)上進(jìn)行 OpenMP 并行優(yōu)化 留

15、意： ：測(cè)試在劃分度為 10 的 9 次方下完成。 參考代碼： #includestdio.h #includeomp.h #define N 1000000000 double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i; #pragma omp parallel for reduction(+:pi) for(i=0; idt; i+) double x=(double)i/dt; pi+=delta/(1.0+x*x); return pi*4; int main() int dx; double pai; /dou

16、ble start,finish; dx=N; double start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtime(); printf(%.8lfn,pai); printf(%lfn,finish-start); return 0; 運(yùn)行結(jié)果如下圖： 串行結(jié)果如下： 提速非常明顯。 1.2.2 OpenMP 并行化+SIMD 向量化 任務(wù)：實(shí)現(xiàn) OpenMP 線程級(jí)和 SIMD 兩級(jí)并行 自動(dòng)向量化代碼如下： #includestdio.h #includeomp.h #define N 1000000000 d

17、ouble get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i; #pragma omp parallel for simd reduction(+:pi) for(i=0; idt; i+) double x=(double)i/dt; pi+=delta/(1.0+x*x); return pi*4; int main() int dx; double pai; dx=N; double start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtime();

18、printf(%.8lfn,pai); printf(%lfn,finish-start); return 0; 留意：自動(dòng)向量化語句為#pragma omp parallel for simd for . 使用編譯語句為：gcc -fopenmp -ma -O3 . 運(yùn)行結(jié)果如下圖： 從結(jié)果看出：有很明顯的提速。 手動(dòng)向量化代碼如下： #includestdio.h #includex86intrin.h #includeomp.h #define N 1000000000 double get_pi(int dt) double pi=0.0; double delta =1.0/dt;

19、double tmp4 _attribute_(aligned(32); _m256d ymm0,ymm1,ymm2,ymm3,ymm4; ymm0=_mm256_set1_pd(1.0); ymm1=_mm256_set1_pd(delta); ymm2=_mm256_set_pd(delta*3,delta*2,delta,0.0); ymm4=_mm256_setzero_pd(); int i; #pragma omp parallel shared(ymm0,ymm1,ymm2) private(i,ymm3,tmp) #pragma omp for reduction(+:pi)

20、for(long int i=0; i=dt-4; i+=4) ymm3= _mm256_set1_pd(double)i*delta); ymm3= _mm256_add_pd(ymm3,ymm2); ymm3= _mm256_mul_pd(ymm3,ymm3); ymm3= _mm256_add_pd(ymm0,ymm3); ymm3= _mm256_div_pd(ymm1,ymm3); /ymm4= _mm256_add_pd(ymm4,ymm3); _mm256_store_pd(tmp,ymm3); pi+=tmp0+tmp1+tmp2+tmp3; /double tmp4 _att

21、ribute_(aligned(32); /_mm256_store_pd(tmp,ymm4); /pi+=tmp0+tmp1+tmp2+tmp3; return pi*4.0; int main() int dx; double pai; dx=N; double start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtime(); printf(%.8lfn,pai); printf(%lfn,finish-start); return 0; 通過對(duì)向量化代碼的分析，各個(gè)向量間的運(yùn)算是沒有任何依靠關(guān)系的，可以直接分線程并

22、行運(yùn)算，但需要留意最終要把各個(gè)線程的運(yùn)算結(jié)果累加。而線程的定義 openmp 的函數(shù) reduction 是沒有方法直接使用（+：）進(jìn)行累加，需要手動(dòng)完成。 引入數(shù)組 tmp 用于將 ymm3 向量分割存放，并累加到 pi 變量，使用openmp 函數(shù) reduction(+:pi)對(duì) pi 變量進(jìn)行累加（詳見代碼） 解決的問題： 并行塊中如何私有化一個(gè)數(shù)組：直接將數(shù)組名稱寫入 private（）函數(shù)中。 曾經(jīng)嘗試將數(shù)組各項(xiàng)都放入 private()函數(shù)中，錯(cuò)誤如下： 多次嘗試后，正確做法如下： 以 tmp4數(shù)組舉例：私有化描述如下： private(tmp); 運(yùn)行結(jié)果如下圖： 手動(dòng)化結(jié)果明

23、顯優(yōu)于自動(dòng)化結(jié)果。手動(dòng)化的修改更符合編寫的程序本身。 二、 矩陣- - 矩陣相乘的 p openmp 優(yōu)化 2.1 個(gè) 編寫一個(gè) " 矩陣- 向量 '或"矩陣-' 矩陣' 相乘的 OpenMP 并 并行程序 ，或其他矩陣運(yùn)算相關(guān)程序。 在 矩陣的驗(yàn)證均在 1024*1024 規(guī)模下完成 矩陣- 矩陣相乘的 openmp 代碼和串行代碼如下： #includestdio.h #includeomp.h #includetime.h #define N 1024 #define n 4 int aNN; int bNN; int cNN; int dNN; int main() int i,j,