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

1、并行實(shí)驗(yàn)報(bào)告一、積分計(jì)算圓周率1.1 積分計(jì)算圓周率的向量?jī)?yōu)化1.1.1 串行版本的設(shè)計(jì)任務(wù)：理解積分求圓周率的方法，將其用C代碼實(shí)現(xiàn)。件1不=4jgdx注意：理論上，dx越小，求得的圓周率越準(zhǔn)確；在計(jì)算機(jī)中由于表示的數(shù)據(jù)是有精度范圍的，如果dx太小，積分次數(shù)過多，誤差積累導(dǎo)致結(jié)果不準(zhǔn)確。以下為串行代碼：#include<stdio.h>#include<time.h>#defineN10000000doubleget_pi(intdt)doublepi=0.0;doubledelta=1.0/dt;inti;for(i=0;i<dt;i+)doublex=(do

2、uble)i/dt;pi+=delta/(1.0+x*x);returnpi*4;intmain()intdx;doublepai;doublestart,finish;dx=N;start=clock();pai=get_pi(dx);finish=clock();printf("%.8lfn",pai);printf("%.8lfSn",(double)(finish-start)/CLOCKS_PER_SEC);return0;時(shí)間運(yùn)行如下：adninladnlni-virtual-machtne:"/Des-ktop$gcc-03tex

3、t1,cadninladniinl-vtrtual-maclitfte/Desktops./a,out3.14159275fl*G267400OSadninladniinl-vtrtua'l-nacliir»e:*-/De5ktop$#/a*out3.141592750*824465695adninl®adniinl-vtrtual-nachifte:-/Desktops*/a*out3.141592750*024628605adninladninl-vtrtual-riacliir*e:/Desktop$I第一次：time=0.02674000s第二次：time=

4、0.02446500s第三次：time=0.02402800s三次平均為：0.02508S1.1.2 SSE向量?jī)?yōu)化版本設(shè)計(jì)任務(wù)：此部分需要給出單精度和雙精度兩個(gè)優(yōu)化版本住思：(1)測(cè)試均在劃分度為10的7次方下完成。以下是SSE雙精度的代碼：#include<stdio.h>#include<x86intrin.h>#include<time.h>#defineN10000000doubleget_pi(intdt)doublepi=0.0;doubledelta=1.0/dt;inti;for(i=0;i<dt;i+)doublex=(double

5、)i/dt;pi+=delta/(1.0+x*x);returnpi*4;doubleget_pi_sse(size_tdt)doublepi=0.0;doubledelta=1.0/dt;_m128dxmm0,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(longinti=0;i<=dt-2;i+=2)xmm3=_mm_set1_pd(double)i*delta);xmm3=_mm_add_pd

6、(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);doubletmp2_attribute_(aligned(16);_mm_store_pd(tmp,xmm4);pi+=tmp0+tmp1/*+tmp2+tmp3*/;returnpi*4.0;intmain()intdx;doublepai;doublestartfinish;dx=N;start=clock();pai=get_pi_sse(dx);f

7、inish=clock();printf("%.8lfn",pai);printf("%.8lfSn",(double)(finish-start)/CLOCKS_PER_SEC);return0;時(shí)間運(yùn)行如下：admtnlgadmtnl-virtual-machtne:«/Desktop$pcc-03testl-SSE.cadminladminl-virtual-machine：w/Desktop$_/a.out3.141592750.000375065adninladmlnl-virtuaL-machtne:«/Desktop$*

8、/a,out3*141592750.007411008adminlQadmlnl-virtual-machtne:"/Desktop$./a.out3.14159275adHinladninil-vtrtual-riachlne:-/Desktop$,/a.out3*141592750.007720005第一次：time=0.00837500s第二次：time=0.00741100s第三次：time=0.00772000s三次平均為：0.00783S以下是SSE單精度的代碼：#include<stdio.h>#include<x86intrin.h>#incl

9、ude<time.h>#defineN10000000floatget_pi_sse(size_tdt)floatpi=0.0;floatdelta=1.0/dt;_m128xmm0,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(longinti=0;i<=dt-4;i+=4)xmm3=_mm_set1_ps(float)i*delta);xmm3=_mm

10、_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);floattmp4_attribute_(aligned(16);_mm_store_ps(tmp,xmm4);pi+=tmp0+tmp1+tmp2+tmp3;returnpi*4.0;intmain()intdx;floatpai;doublestartfinish;dx=N;start=clock();pai=get_pi_sse(dx);

11、finish=clock();printf("%.8fn",pai);printf("%.8lfSn",(double)(finish-start)/CLOCKS_PER_SEC);return0;時(shí)間運(yùn)行如下:adminladntnl-virtual-machines/Desktop$gcc-03testl-float-SSEicadminladninl-virtual-mathtne:/Desktop$./a.out3,145040510.OQ4O6100Sadninladminl-vi.rtuaL-mach1.ne:"/Desktop5o

12、utL14504051kDO4Z64OO5adminl®admini-virtual-machtne:"/Desktop$*/a.out3.14504051O.O6437600S第一次：time=0.00406100s第二次：time=0.00426400s第三次：time=0.00437600s三次平均為：0.00423s1.1.3 AVX向量?jī)?yōu)化版本設(shè)計(jì)任務(wù)：此部分需要給出單精度和雙精度兩個(gè)優(yōu)化版本、/、4'1、廠注息：(1)測(cè)試均在劃分度為10的7次方下完成。(2)在編譯時(shí)需要加-mavx編譯選項(xiàng)，才能啟用AVX指令集,否則默認(rèn)SSE指令集(3)理論上，向量版

13、本對(duì)比SSE版本和串行版本有明顯加速，單精度版本速度明顯優(yōu)于雙精度，速度接近雙精度的兩倍。以下是AVX雙精度的代碼：#include<stdio.h>#include<x86intrin.h>#include<time.h>#defineN10000000/*doubleget_pi(intdt)doublepi=0.0;doubledelta=1.0/dt;inti;for(i=0;i<dt;i+)doublex=(double)i/dt;pi+=delta/(1.0+x*x);returnpi*4;*/doubleget_pi_avx(size_t

14、dt)doublepi=0.0;doubledelta=1.0/dt;_m256dymm0,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(longinti=0;i<=dt-4;i+=4)ymm3=_mm256_set1_pd(double)i*delta);ymm3=_mm256_add_pd(ymm3,ymm2);ymm3=_mm256_mu

15、l_pd(ymm3,ymm3);ymm3=_mm256_add_pd(ymm0,ymm3);ymm3=_mm256_div_pd(ymm1,ymm3);ymm4=_mm256_add_pd(ymm4,ymm3);doubletmp4_attribute_(aligned(32);_mm256_store_pd(tmp,ymm4);pi+=tmp0+tmp1+tmp2+tmp3;returnpi*4.0;intmain()intdx;doublepai;doublestart,finish;dx=N;start=clock();pai=get_pi_avx(dx);finish=clock();

16、printf("%.8lfn",pai);printf("%.8lfSn",(double)(finish-start)/CLOCKS_PER_SEC);return0;時(shí)間運(yùn)行如下：adr)tnl(aadMtftl-virtual-nachtne:-/Desktop$gcc-navx-03testl-AVX.cadnitnl(3admtrkl-vtrtuaL-nachtne:-/Desktops./a.out3.141592750.007202005adnlnlgadnt>l-virtual-nachine;-/Desktop$./a.out3.

17、141592750.0065980GSadninlQadminl-virtual-Fiachtne:"/Desktop$./a.out3.141592750.OO67O6QOS第一次：time=0.00720200s第二次：time=0.00659800s第三次：time=0.00670600s三次平均為：0.00683S以下是AVX單精度的代碼:時(shí)間運(yùn)行如下:admtftK3admtnl-vtrtLal>ma£htne:Desktopsgcc-03testl-float*AVX.cadriirkldritnl-virtual-machine:Desktops*/a-

18、out3.141733170.662342G05admtnladmtnl-virtual-machine:/Desktop$./a.outB.14173317p.ee2342eosadminlfiadmi.nl-virtual-machine:Desktops'/3.out14173317O.OOZ300eOS第一次：time=0.00234200s第二次：time=0.00234200s第三次：time=0.00230000s三次平均為：0.002328s由以上實(shí)驗(yàn)統(tǒng)計(jì)得出結(jié)論：AVX-float=0.002328sAVX-double=0.00683sssE-float=0.004

19、23sssE-double=0.00783s基本符合規(guī)律：（以下為速度比較）AVX-float>AVX-doublessE-float>ssE-double>serial1.2 積分計(jì)算圓周率的OpenMP優(yōu)化1.2.1 OpenMP并行化任務(wù)：在串行代碼的基礎(chǔ)上進(jìn)行OpenMP并行優(yōu)化注意：測(cè)試在劃分度為10的9次方下完成。參考代碼：#include<stdio.h>#include<omp.h>#defineN1000000000doubleget_pi(intdt)doublepi=0.0;doubledelta=1.0/dt;inti;#pra

20、gmaompparallelforreduction(+:pi)for(i=0;i<dt;i+)doublex=(double)i/dt;pi+=delta/(1.0+x*x);returnpi*4;intmain()intdx;doublepai;/doublestartfinish;dx=N;doublestart=omp_get_wtime();pai=get_pi(dx);doublefinish=omp_get_wtime();printf("%.8lfn",pai);printf("%lfn",finish-start);return0

21、;運(yùn)行結(jié)果如下圖:adntnugadfttrl-virtualmachine:*/Desktopsgcc-fopennp-03test2*cadntnladntn1-virtual-riachtne:"/Desktops*/a.out3,141592051.Z54093odninlOodninl-virtuol-Hachtne:-/Desktops.，曰.out3.141592651.260fi4fiadntnl(3adntfil-v1.rtual-Fiachtne:«/Desktop$./a.out3.141592651.244996串行結(jié)果如下:dHinidadmtnl

22、-vtrtuaL-nachlne:-/De&ktapSgcc=03textileadninlQadninl-virtual-macMne:Desktop$-out3.141592652.3297720OSadntnlQadmtnl-vtrtuaL-machtnei/Desktops./a.out3.141592652.334373685ddntnladmtnl-vtrtudl>machtne:/Desktops./a.out3.141的26£2,3347Q1O0S提速十分明顯。1.2.2 OpenMP并行化+SIMD向量化任務(wù)：實(shí)現(xiàn)OpenMP線程級(jí)和SIMD兩級(jí)并行自

23、動(dòng)向量化代碼如下：#include<stdio.h>#include<omp.h>#defineN1000000000doubleget_pi(intdt)doublepi=0.0;doubledelta=1.0/dt;inti;#pragmaompparallelforsimdreduction(+:pi)for(i=0;i<dt;i+)doublex=(double)i/dt;pi+=delta/(1.0+x*x);returnpi*4;intmain()intdx;doublepai;dx=N;doublestart=omp_get_wtime();pai=

24、get_pi(dx);doublefinish=omp_get_wtime();printf("%.8lfn",pai);printf("%lfn",finish-start);return0;注意：自動(dòng)向量化語句為#pragmaompparallelforsimdfor.使用編譯語句為：gcc-fopenmp-mavx-O3.運(yùn)行結(jié)果如下圖：adntni«6dntni-vtrtual-nachtne:*/Deskto|j5gcc-fopenmp-mavx-03testZ-stmdicadritnladninl-virtual-nachint:

25、-*-/E>esktop5*/d*out3.141S926S6.663866adminladninl-virtual-machire：*/Desktop?./w.out3*14159265。"刃444adninlOadnitnl-virtueI-naGhine;-/Desktop5.out3.141592658612從結(jié)果看出：有很明顯的提速手動(dòng)向量化代碼如下：#include<stdio.h>#include<x86intrin.h>#include<omp.h>#defineN1000000000doubleget_pi(intdt)do

26、ublepi=0.0;doubledelta=1.0/dt;doubletmp4_attribute_(aligned(32);_m256dymm0,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();inti;#pragmaompparallelshared(ymm0,ymm1,ymm2)private(i,ymm3,tmp)#pragmaompforreduc

27、tion(+:pi)for(longinti=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;/doubletmp4_attribute_

28、(aligned(32);/_mm256_store_pd(tmp,ymm4);/pi+=tmp0+tmp1+tmp2+tmp3;returnpi*4.0;intmain()intdx;doublepai;dx=N;doublestart=omp_get_wtime();pai=get_pi(dx);doublefinish=omp_get_wtime();printf("%.8lfn",pai);printf("%lfn",finish-start);return0;通過對(duì)向量化代碼的分析，各個(gè)向量間的運(yùn)算是沒有任何依賴關(guān)系的，可以直接分線程并行運(yùn)算，

29、但需要注意最后要把各個(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ò)誤如下：adntnidrttnl-virtual-rtachtnei-/De£ktop$get-mavx03test2-handxtestZ-h«nd+c:Infun

30、ction*gct_pi1:testZ-hand.c:16:63:r<rexpected*)Jbeforetokentpragnaonpparallelshared(y«n6lynFilfynn2)private(ifym3,tnp6,t叩1,tnip2)tnp3)多次嘗試后，正確做法如下：以tmp4數(shù)組舉例：私有化描述如下:private(tmp);運(yùn)行結(jié)果如下圖:adHtnladni.ftl-Artual-Machtflc:-/Desktopsgecfopennp-navxO2test2KanoTEadninl4dninlL-virtual-nachtne:"/D

31、esktop5Ja*out3.14159265<3.adntnladnitnl-vtrtuaI-machine:"/Desktops=/a-out3.141592650.425336adriinlffladntnl-virtual-ndchine;"-/Desktops，/a.out3.141592650.4J1496手動(dòng)化結(jié)果明顯優(yōu)于自動(dòng)化結(jié)果。手動(dòng)化的修改更符合編寫的程序本身。二、矩陣-矩陣相乘的openmp優(yōu)化2.1編寫一個(gè)“矩陣-向量”或“矩陣矩陣”相乘的OpenMP弁行程序，或其他矩陣運(yùn)算相關(guān)程序。矩陣的驗(yàn)證均在1024*1024規(guī)模下完成矩陣-矩陣相乘的openmp代碼和串行代碼如下：#include<stdio.h>#include<omp.h>#include<time.h>#defineN1024#definen4intaNN;intbNN;intcNN;intdNN;intmain()inti,j,k;for(i=0;i<N;i+)for(j=0;j<N;j+)a皿=1；b皿=1；cij=0；d皿=0;doublestart1=clock();for(i=0;i<N;i+)for(j=0;j<N;j+)for