OpenMP并行实验报告(同名2333)(共16页)

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

2、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第二次：time=0.02446500S第三次：time=0.02402800S三次平均为：0.02508S1.1.2 SSE向

3、量优化版本设计任务：此部分需要给出单精度和双精度两个优化版本。注意：（1）测试均在划分度为10的7次方下完成。以下是SSE双精度的代码：#include#include#include#define N 10000000double 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; pi+=delta/(1.0+x*x); return pi*4;double get_pi_sse(size_t dt) double pi=0.0; double

4、 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_pd(double)i*delta); xmm3= _mm_add_pd(xmm3,xmm2); xmm3= _mm_mul_pd(xmm3,xmm3); xmm3= _mm_add_pd(xmm0,xmm3)

5、; 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; double start,finish; dx=N; start=clock(); pai=get_pi_sse(dx); finish=clock(); printf(%.8lfn,pai); printf(

6、%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); return 0;时间运行如下：第一次：time=0.00837500S第二次：time=0.00741100S第三次：time=0.00772000S三次平均为：0.00783S以下是SSE单精度的代码：#include#include#include#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_set

7、1_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(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_p

8、s(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 start,finish; dx=N; start=clock(); pai=get_pi_sse(dx); finish=clock(); printf(%.8fn,pai); printf(%.8lfSn,(double)(finish-start)/CLOCKS_PER_SEC); retur

9、n 0;时间运行如下：第一次：time=0.00406100S第二次：time=0.00426400S第三次：time=0.00437600S三次平均为：0.00423S1.1.3 AVX向量优化版本设计任务：此部分需要给出单精度和双精度两个优化版本注意：（1）测试均在划分度为10的7次方下完成。（2）在编译时需要加-mavx 编译选项，才能启用AVX指令集，否则默认SSE指令集（3）理论上，向量版本对比SSE版本和串行版本有明显加速，单精度版本速度明显优于双精度，速度接近双精度的两倍。以下是AVX双精度的代码：#include#include#include#define N 1000000

10、0/*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; 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); ym

11、m2=_mm256_set_pd(delta*3,delta*2,delta,0.0); ymm4=_mm256_setzero_pd(); for(long int i=0; i AVX-double SSE-float SSE-double serial1.2 积分计算圆周率的OpenMP优化1.2.1 OpenMP并行化任务：在串行代码的基础上进行OpenMP并行优化注意：测试在划分度为10的9次方下完成。参考代码：#include#include#define N 1000000000double get_pi(int dt) double pi=0.0; double delta =

12、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; /double start,finish; dx=N; double start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtime(); printf(%.8lfn,pai); printf(%lfn,fin

13、ish-start); return 0;运行结果如下图：串行结果如下：提速十分明显。1.2.2 OpenMP并行化+SIMD向量化任务：实现OpenMP线程级和SIMD两级并行自动向量化代码如下：#include#include#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; idt; i+) double x=(double)i/dt; pi+=delta/(1.

14、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(); printf(%.8lfn,pai); printf(%lfn,finish-start); return 0;注意：自动向量化语句为#pragma omp parallel for simd for .使用编译语句为：gcc -fopenmp -mavx -O3 .运行结果如下图：从结果看出：有很明显的提速。手动向量化代码如下：

15、#include#include#include#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=_mm256_set_pd(delta*3,delta*2,delta,0.0); ymm4=_mm256_setzero_pd();

16、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,ymm3); ymm3= _mm256_add_pd(ymm0,ymm3); ymm3= _mm256_div_pd(ymm1,ymm3);

17、/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; double start=omp_get_wtime(); pai=get_pi(dx); double finish=omp_get_wtim

18、e(); printf(%.8lfn,pai); printf(%lfn,finish-start); return 0;通过对向量化代码的分析，各个向量间的运算是没有任何依赖关系的，可以直接分线程并行运算，但需要注意最后要把各个线程的运算结果累加。而线程的定义openmp的函数reduction是没有办法直接使用（+：）进行累加，需要手动完成。引入数组tmp用于将ymm3向量分割存放，并累加到pi变量，使用openmp函数reduction(+:pi)对pi变量进行累加（详见代码）解决的问题：并行块中如何私有化一个数组：直接将数组名称写入private（）函数中。曾经尝试将数组各项都放入pr

19、ivate()函数中，错误如下：多次尝试后，正确做法如下：以tmp4数组举例：私有化描述如下：private(tmp);运行结果如下图：手动化结果明显优于自动化结果。手动化的修改更符合编写的程序本身。二、 矩阵-矩阵相乘的openmp优化2.1 编写一个“矩阵-向量”或“矩阵-矩阵”相乘的 OpenMP 并行程序，或其他矩阵运算相关程序。矩阵的验证均在1024*1024规模下完成矩阵-矩阵相乘的openmp代码和串行代码如下：#include#include#include#define N 1024#define n 4 int aNN; int bNN; int cNN; int dNN;int main() int i,j,k; for(i=0;iN;i+) for(j=0;jN;j+) aij=1; bij=1; cij=0; dij=0; double start1=clock(); for(i=0;iN;i+) for