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

2、ble 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第二次：tim

3、e=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>#define N 10000000double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i; for(i=0; i<dt; i

4、+) double 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= _

5、mm_set1_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; doub

6、le pai; 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.00772000S三次平均為：0.00783S以下是SSE單精度的代碼：#include<

7、stdio.h>#include<x86intrin.h>#include<time.h>#define N 10000000float 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

8、 i=0; i<=dt-4; i+=4) xmm3= _mm_set1_ps(float)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.

9、0;int main() int dx; float pai; double start,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.00423S1.

10、1.3 AVX向量?jī)?yōu)化版本設(shè)計(jì)任務(wù)：此部分需要給出單精度和雙精度兩個(gè)優(yōu)化版本注意：（1）測(cè)試均在劃分度為10的7次方下完成。（2）在編譯時(shí)需要加-mavx 編譯選項(xiàng)，才能啟用AVX指令集，否則默認(rèn)SSE指令集（3）理論上，向量版本對(duì)比SSE版本和串行版本有明顯加速，單精度版本速度明顯優(yōu)于雙精度，速度接近雙精度的兩倍。以下是AVX雙精度的代碼：#include<stdio.h>#include<x86intrin.h>#include<time.h>#define N 10000000/*double get_pi(int dt) double pi=0.0;

11、 double delta =1.0/dt; int i; for(i=0; i<dt; i+) double x=(double)i/dt; pi+=delta/(1.0+x*x); return pi*4;*/double get_pi_avx(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,de

12、lta,0.0); ymm4=_mm256_setzero_pd(); 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); double tmp4 _attribute_(aligne

13、d(32); _mm256_store_pd(tmp,ymm4); pi+=tmp0+tmp1+tmp2+tmp3; return pi*4.0;int main() int dx; double pai; double start,finish; dx=N; start=clock(); pai=get_pi_avx(dx); finish=clock(); printf("%.8lfn",pai); printf("%.8lfSn",(double)(finish-start)/CLOCKS_PER_SEC); return 0;時(shí)間運(yùn)行如下：第一次

14、：time=0.00720200S第二次：time=0.00659800S第三次：time=0.00670600S三次平均為：0.00683S以下是AVX單精度的代碼：時(shí)間運(yùn)行如下：第一次：time=0.00234200S第二次：time=0.00234200S第三次：time=0.00230000S三次平均為：0.002328S由以上實(shí)驗(yàn)統(tǒng)計(jì)得出結(jié)論：AVX-float=0.002328 SAVX-double=0.00683 SSSE-float=0.00423 SSSE-double= 0.00783 S基本符合規(guī)律：（以下為速度比較）AVX-float > AVX-double

15、SSE-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>#define N 1000000000double get_pi(int dt) double pi=0.0; double delta =1.0/dt; int i;#pragma omp parallel for reduction(+:pi) for(i=0;

16、i<dt; 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; 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é)果如下：

17、提速十分明顯。1.2.2 OpenMP并行化+SIMD向量化任務(wù)：實(shí)現(xiàn)OpenMP線程級(jí)和SIMD兩級(jí)并行自動(dòng)向量化代碼如下：#include<stdio.h>#include<omp.h>#define N 1000000000double 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; i<dt; i+) double x=(double)i/dt; pi+=delta/(1.0+x*

18、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(); printf("%.8lfn",pai); printf("%lfn",finish-start); return 0;注意：自動(dòng)向量化語(yǔ)句為#pragma omp parallel for simd for .使用編譯語(yǔ)句為：gcc -fopenmp -mavx -O3 .運(yùn)行結(jié)果如下圖：從結(jié)果看

19、出：有很明顯的提速。手動(dòng)向量化代碼如下：#include<stdio.h>#include<x86intrin.h>#include<omp.h>#define N 1000000000double get_pi(int dt) double pi=0.0; double delta =1.0/dt; double tmp4 _attribute_(aligned(32); _m256d ymm0,ymm1,ymm2,ymm3,ymm4; ymm0=_mm256_set1_pd(1.0); ymm1=_mm256_set1_pd(delta); ymm2=_

20、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) 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,

21、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 _attribute_(aligned(32); /_mm256_store_pd(tmp,ymm4); /pi+=tmp0+tmp1+tmp2+tmp3; return pi*4.0;int main() int dx; double pai; dx=N; d

22、ouble start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtime(); printf("%.8lfn",pai); printf("%lfn",finish-start); return 0;通過(guò)對(duì)向量化代碼的分析，各個(gè)向量間的運(yùn)算是沒(méi)有任何依賴關(guān)系的，可以直接分線程并行運(yùn)算，但需要注意最后要把各個(gè)線程的運(yùn)算結(jié)果累加。而線程的定義openmp的函數(shù)reduction是沒(méi)有辦法直接使用（+：）進(jìn)行累加，需要手動(dòng)完成。引入數(shù)組tmp用于將ymm3向量分割存放，并累加到pi變量，使

23、用openmp函數(shù)reduction(+:pi)對(duì)pi變量進(jìn)行累加（詳見(jiàn)代碼）解決的問(wèn)題：并行塊中如何私有化一個(gè)數(shù)組：直接將數(shù)組名稱寫(xiě)入private（）函數(shù)中。曾經(jīng)嘗試將數(shù)組各項(xiàng)都放入private()函數(shù)中，錯(cuò)誤如下：多次嘗試后，正確做法如下：以tmp4數(shù)組舉例：私有化描述如下：private(tmp);運(yùn)行結(jié)果如下圖：手動(dòng)化結(jié)果明顯優(yōu)于自動(dòng)化結(jié)果。手動(dòng)化的修改更符合編寫(xiě)的程序本身。二、 矩陣-矩陣相乘的openmp優(yōu)化2.1 編寫(xiě)一個(gè)“矩陣-向量”或“矩陣-矩陣”相乘的 OpenMP 并行程序，或其他矩陣運(yùn)算相關(guān)程序。矩陣的驗(yàn)證均在1024*1024規(guī)模下完成矩陣-矩陣相乘的openmp代碼和串行代碼如下：#include<stdio.h>#include<omp.h>#include<time.h>#define N 1024#define n 4 int aNN; int bNN; int cNN; int dNN;int main() int i,j,k; for(i=0;i<N;i+) for(j=0;j<N;j+) aij=1; bij=1; cij=0; dij=0; double start1=clock(); for(i=