OpenMP并行实验报告

1、并行实验报告 一、 积分计算圆周率1.1 积分计算圆周率的向量优化1.1.1 串行版本的设计任务：理解积分求圆周率的方法，将其用C代码实现。注意：理论上，dx越小，求得的圆周率越准确；在计算机中由于表示的数据是有精度范围的，如果dx太小，积分次数过多，误差积累导致结果不准确。以下为串行代码：#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;时间运行如下：第一次：time=0.02674000S第二次：tim

3、e=0.02446500S第三次：time=0.02402800S三次平均为：0.02508S1.1.2 SSE向量优化版本设计任务：此部分需要给出单精度和双精度两个优化版本。注意：（1）测试均在划分度为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;时间运行如下：第一次：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;时间运行如下：第一次：time=0.00406100S第二次：time=0.00426400S第三次：time=0.00437600S三次平均为：0.00423S1.

10、1.3 AVX向量优化版本设计任务：此部分需要给出单精度和双精度两个优化版本注意：（1）测试均在划分度为10的7次方下完成。（2）在编译时需要加-mavx 编译选项，才能启用AVX指令集，否则默认SSE指令集（3）理论上，向量版本对比SSE版本和串行版本有明显加速，单精度版本速度明显优于双精度，速度接近双精度的两倍。以下是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;时间运行如下：第一次

14、：time=0.00720200S第二次：time=0.00659800S第三次：time=0.00670600S三次平均为：0.00683S以下是AVX单精度的代码：时间运行如下：第一次：time=0.00234200S第二次：time=0.00234200S第三次：time=0.00230000S三次平均为：0.002328S由以上实验统计得出结论：AVX-float=0.002328 SAVX-double=0.00683 SSSE-float=0.00423 SSSE-double= 0.00783 S基本符合规律：（以下为速度比较）AVX-float > AVX-double

15、SSE-float > SSE-double > serial1.2 积分计算圆周率的OpenMP优化1.2.1 OpenMP并行化任务：在串行代码的基础上进行OpenMP并行优化注意：测试在划分度为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;运行结果如下图：串行结果如下：

17、提速十分明显。1.2.2 OpenMP并行化+SIMD向量化任务：实现OpenMP线程级和SIMD两级并行自动向量化代码如下：#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;注意：自动向量化语句为#pragma omp parallel for simd for .使用编译语句为：gcc -fopenmp -mavx -O3 .运行结果如下图：从结果看

19、出：有很明显的提速。手动向量化代码如下：#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;通过对向量化代码的分析，各个向量间的运算是没有任何依赖关系的，可以直接分线程并行运算，但需要注意最后要把各个线程的运算结果累加。而线程的定义openmp的函数reduction是没有办法直接使用（+：）进行累加，需要手动完成。引入数组tmp用于将ymm3向量分割存放，并累加到pi变量，使

23、用openmp函数reduction(+:pi)对pi变量进行累加（详见代码）解决的问题：并行块中如何私有化一个数组：直接将数组名称写入private（）函数中。曾经尝试将数组各项都放入private()函数中，错误如下：多次尝试后，正确做法如下：以tmp4数组举例：私有化描述如下：private(tmp);运行结果如下图：手动化结果明显优于自动化结果。手动化的修改更符合编写的程序本身。二、 矩阵-矩阵相乘的openmp优化2.1 编写一个“矩阵-向量”或“矩阵-矩阵”相乘的 OpenMP 并行程序，或其他矩阵运算相关程序。矩阵的验证均在1024*1024规模下完成矩阵-矩阵相乘的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=