矩阵乘法的OpenMP并行程序

绝不收费

于 2024-07-17 22:49:32 发布

阅读量676

点赞数 10

分类专栏：并行计算实验文章标签： c++

本文链接：https://blog.csdn.net/weixin_60705585/article/details/140507576

版权

并行计算实验专栏收录该内容

4 篇文章 0 订阅

订阅专栏

深圳大学并行计算实验一

一、实验目的

1. 学会编写简单的OpenMP程序；

2. 掌握for编译制导语句；

3. 对并行程序进行简单的性能分析；

二、实验环境

1. 硬件环境：64核CPU、256GB内存的共享内存并行计算平台；

2. 软件环境：Ubuntu Linux、gcc、g++（g++ -O3 -fopenmp -o a.out a.cc）；

3. 远程登录：本地PowerShell中执行ssh bxjs@hpc.szu.edu.cn；

4. 传输文件：本地PowerShell中执行scp c:\a.cpp bxjs@hpc.szu.edu.cn:/home/bxjs/或ftp://hpc.szu.edu.cn。

三、实验内容

1. 用OpenMP编写两个n阶方阵a和b的并行相乘程序，结果存放在方阵c中，其中矩阵乘法部分用for编译制导语句实现并行化操作。为了验证结果正确性，将矩阵乘法的串行计算结果存放在方阵d中，并比较是否与c相等。

代码实现思路：

首先，我们考虑用double类型的二维数组构建a、b、c、d四个矩阵，用于存储四个矩阵的值。接着，我们考虑对a、b进行随机的初始化。每次进行矩阵乘法前，对a、b进行随机初始化，才能避免分配的内存中的值对实验产生影响，更好地验证并行化结果的正确性。对此，我采用了cstdlib库中rand()函数获取一个0到RAND_MAX之间的随机整数，再将随机数乘以1.0变成double类型，再除以RAND_MAX，最终获得一个0到1之间的随机的浮点数，再将该随机浮点数赋值给a、b矩阵。同时，对c、d矩阵进行初始化，为0。

然后，我考虑使用了OpenMP库中，omp_set_num_threads(int tp)函数用于设置并行区域的线程数为tp。再用omp_get_wtime()函数，获取当前线程距离线程线程开始的时间差，赋值给double类型的变量t0，作为并行域计算的起始时间。接着，使用“#pragma omp parallel for”编译指导语句，构建一个for循环并行域，然后在并行域中实现一个三层for循环，实现a矩阵和b矩阵的矩阵乘法，并将结果存储到c矩阵中，实现矩阵乘法的并行化。

接着，再次使用omp_get_wtime()函数，获取当前线程距离线程线程开始的时间差，赋值给double类型的变量t1，作为并行域计算的结束时间。继而，通过t1-t0就可以获得该并行域的运行时间。紧接着，通过使用一个三层for循环，串行实现a矩阵和b矩阵的矩阵乘法，并将结果存储到d矩阵中。

最后，采用两层for循环比较c矩阵和d矩阵的中数值，判断是否相等，来判断矩阵乘法并行化的正确性。但是，考虑到矩阵的数组是浮点类型，所以在判断c矩阵和d矩阵的中数值是否相等时，为了避免精度丢引起判断错误，只需保证误差的一定小的范围内，即可认为数值相等。

初步代码实现：

//初始化设置 
int i, j, k;
for (i = 0; i < n; i++)
	for (j = 0; j < n; j++){
		a[i][j] = rand() * 1.0 / RAND_MAX;
		b[i][j] = rand() * 1.0 / RAND_MAX;
		c[i][j] = d[i][j] = 0; 
	}
omp_set_num_threads(tp);
double t0 = omp_get_wtime();
// 并行代码
#pragma omp parallel for private(i, j, k) shared(a, b, c)
for (i = 0; i < n; i++)
	for (j = 0; j < n; j++)
		for (k = 0; k < n; k++){
			#pragma omp atomic
			c[j][k] += a[j][i] * b[i][k];
		}
double t1 = omp_get_wtime();
cout << "num_thread: "<<tp<<", Parallel time: " << t1 - t0 << " seconds" << endl;
// 验证代码
for (i = 0; i < n; i++)
	for (j = 0; j < n; j++)
		for (k = 0; k < n; k++)
			d[i][j] += a[i][k] * b[k][j];
for (i = 0; i < n; i++)
	for (j = 0; j < n; j++)
		if (abs(c[i][j] - d[i][j]) > PanDuan_WuChai){
			cout << "Results are not equal!" << endl;
			return 0;
		}
cout << "Results are equal!" << endl;

2. 测试并行程序在不同线程数下的执行时间和加速比（与线程数=1时的执行时间相比）。其中，n固定为1000，线程数分别取1、2、4、8、16、32、64时，为减少误差，每项实验进行5次，取平均值作为实验结果。

为了避免代码多次修改，这里可以将对上述1中的代码进行封装成函数“double juzheng(int tp)”，传入参数tp，为程序并行域的线程数，并将并行域的运行时间作为返回值进行返回。

同时，为了避免程序的多次间断的运行，导致程序运行时CPU的状态差，导致实验获取的结果误差较大，这里考虑，创建数组提前存储“1、2、4、8、16、32、64”，然后通过for循环将数值传入函数“double juzheng(int tp)”中，获取运行时间，并计算出对应5次实验的平均值。从而，避免了程序的多次间断的运行带来的误差。具体的代码在下列“四、代码描述”中可见。

四、代码描述

(这部分当时没有写分析被扣分了)

#include <omp.h>
#include <iostream>
#include <cmath>
using namespace std;
const int n = 1000;//矩阵大小
double a[n][n], b[n][n], c[n][n], d[n][n];//矩阵 
const int ThreadNum_CanShu=7;//线程数组的大小
const int set_thread_num[]={1,2,4,8,16,32,64};//线程数
const int XunHuan_CanShu=5;//循环次数  
const double PanDuan_WuChai=1e-10;//判断 c、d 相等的允许误差 
//tp个线程的矩阵乘法
double juzheng(int tp){
	//初始化设置 
	int i, j, k;
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++){
			a[i][j] = rand() * 1.0 / RAND_MAX;
			b[i][j] = rand() * 1.0 / RAND_MAX;
			c[i][j] = d[i][j] = 0; 
		}
	omp_set_num_threads(tp);
	double t0 = omp_get_wtime();
	// 并行代码
          #pragma omp parallel for private(i, j, k) shared(a, b, c)
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			for (k = 0; k < n; k++){
				c[i][j] += a[i][k] * b[k][j]; 
			}
	double t1 = omp_get_wtime();
	cout << "num_thread: "<<tp<<", Parallel time: " << t1 - t0 << " seconds" << endl;
	// 验证代码
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			for (k = 0; k < n; k++)
				d[i][j] += a[i][k] * b[k][j];
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			if (abs(c[i][j] - d[i][j]) > PanDuan_WuChai){
				cout << "Results are not equal!" << endl;
				cout<<c[i][j]<<" "<<d[i][j]<<endl;
				return 0;
			}
	cout << "Results are equal!" << endl;
	return t1-t0;
} 
int main()
{
	cout<<"CPU核心数："<<omp_get_num_procs()<<endl;
	cout<<"最大线程数："<<omp_get_thread_limit()<<endl; 
	cout<<"----------------------------------"<<endl; 
	for(int i=0;i<ThreadNum_CanShu;i++){
		double pingjunzhi=0;
		for(int j=0;j<XunHuan_CanShu;j++){
			double temp=juzheng(set_thread_num[i]);
			pingjunzhi+=temp;
		}
		cout<<"omp_set_num_threads="<<i<<" 平均运行时间："<< pingjunzhi/XunHuan_CanShu<<endl;
		cout<<"----------------------------------"<<endl; 
	} 
	return 0;
}

五、实验结果和分析

通过在文件目录下输入“ftp://hpc.szu.edu.cn”，将本地的a.cc的代码文件赋值到“/2021150233”的目录下，通过“g++ -O3 -fopenmp -o a.out a.cc”将a.cc文件编译为a.out文件。然后通过“./a.out”运行a.out文件，最终执行的结果如下图1所示，并将执行的结果记录进表1中，同时计算出加速比。

图1 a.out文件的执行结果

表1 并行程序在不同线程数下的执行时间（秒）和加速比

线程数执行时间/s	1	2	4	8	16	32	64
第1次	35.5256	24.6552	11.5390	6.4151	3.1810	1.8856	0.9603
第2次	44.6757	24.9403	12.0915	6.4780	3.5551	1.7985	0.8535
第3次	44.6589	23.8688	12.3917	6.2319	2.7552	1.8854	0.9076
第4次	44.8472	24.9445	11.8045	5.7337	3.2462	1.8707	1.0183
第5次	45.1704	22.4681	11.3752	6.5860	28865	1.6242	0.9514
平均值	42.9756	24.1754	11.8404	6.2889	3.1248	1.8129	0.9382
加速比	1.0000	1.7777	3.6296	6.8336	13.7530	23.7054	45.8064

根据上述表1的平均值和加速比数据，绘制出下面的折线图：

图2 平均值和加速比的折线图

根据对上述，表1和图2的分析可知，随着线程数的增加，矩阵乘法的运行时间逐渐下降，加速比逐渐上升。即，在CPU核心数量的范围内，通过增加线程数，分割任务进行并行运算可以实现对复杂计算的优化。

同时，观测到加速比与线程数存在线性关系，但是并不是预期的“加速比==线程数”的关系，而是“加速比==a*线程数(0<a<1)”的关系。这主要是由于随着线程数的增加，CPU对任务调度的开销也随之增加，导致出现“加速比==a*线程数(0<a<1)”的关系。

六、实验结论

在实验过程中，我尝试对“#pragma omp parallel后跟 #pragma omp for”的编译指导语句实现的功能进行研究时，发现一条与此极为相似的编译指导语句“#pragma omp parallel for”，然后我对这两条编译指导语句进行了研究：

#pragma omp parallel for和#pragma omp parallel后跟#pragma omp for在功能上是相同的，都是用来并行化for循环的。如果只需要并行化一个for循环，那么#pragma omp parallel for会更简洁。如果你需要在并行区域中执行更复杂的操作，那么#pragma omp parallel后跟#pragma omp for可能会更合适。但是，你需要注意的是，使用#pragma omp parallel后跟#pragma omp for时，必须确保所有的线程都能执行到#pragma omp for，否则可能会导致程序的行为不正确。

在实验过程中，我还尝试了对并行运算中的“c[i][k] += a[i][j] * b[j][k]”操作进行研究，我在思考：

这句代码应该不是原子操作的啊？从汇编的角度看，这行代码在x86指令集下，最少也需要两条指令才能完成，amd指令集就更不止两条指令了，那行代码基本就应该不算原子操作的？原子操作的判断是通过代码编译成的汇编指令数为一条来判断的？

在探索这些问题的过程中，我也收获颇丰：

原子操作确实可以通过判断具体的高级语言的代码是否被编译成单独的一条汇编指令来判断是否是原子操作，这是因为每条汇编指令都对应一条机器指令，对能CPU进行操作的集合。但是优于我是使用高级语言进行程序的编写，中间需要有编译器将我代码编译为汇编代码，但是这个过程却是我不能控制的。编译是可能出于优化，采用多条指令实现本来可以只用一条指令实现的操作。所以在高级语言的层面，除非通过编译器提供的方法或者本来的规定的语言的规则，否则，我们是无法轻易判断一个操作是原子操作的。

根据上述的探索，可以断定并行运算中的“c[i][k] += a[i][j] * b[j][k]”操作肯定不是原子操作，但是我通过各种手段，我尝试将并行域的线程数开到很大很大（比n*n大），目的很明确，就是使得可以有不同线程可以同时操作i和k相等的c[i][k]。因为通过查找资料了解到这个“c[i][k] += a[i][j] * b[j][k]”操作正常情况下的汇编代码涉及到读取c[i][k]，然后加上乘积，然后写回。

而这其中，如果读取成功后，在进行乘积和写回的过程中，如果被中断，那么就很可能出现‘脏读’‘脏写’等问题，会导致最后的矩阵乘法的结果错误。但是当我把线程开到比n*n大很多后，还是没有发现有矩阵乘法的结果错误的输出。然后，我通过调试不断调试和查询资料后，惊奇的发现：

#paragam omp for 默认情况下，只会对最近的一层循环进行并行优化，不会对嵌套的多层循环进行进行并行优化。也就是这个例子中的线程开出超过n后，也只会根据最外层的i循环分配出最多n个线程的并行工作，这种情况下，就根本无法实现多个线程同时操作i和k相等的c[i][k]的情况，得到上面推理的结果。

然后，我就探索如何能使得不同线程可以同时操作i和k相等的c[i][k]，出现上面的预期结果，在这个过程中我找到了两种方法：

调换“四、代码描述”中所示的i、j、k的循环逻辑为c[j][k] += a[j][i] * b[i][k]
在编译指导语句的末尾加上collapse(3)（前提：n层循环之间不能有其他语句；循环域要规则，例如，不能第二层从j=i开始）会根据参数将循环进行尽可能的展开，能展开几层就展开几层，最多展成一层大循环，对这个大循环进行并行优化。

最后，我对omp中如何设置代码为原子操作进行了探索，发现：

可以使用#pragma omp atomic指令来确保一段代码作为原子操作执行。但是#pragma omp atomic只能保证单个内存位置的读-改-写操作是原子的。如果需要保护更复杂的代码段，需要使用#pragma omp critical或者#pragma omp lock。但是因为可能涉及到对内存的锁定等问题，所以可能会牺牲一定的性能。

通过本次实验，我掌握了OpenMP编程的基本知识和应用，了解了并行计算的优势和影响因素。实验结果表明，通过合理设置线程数，可以有效提高计算性能，但也需要考虑线程数增加导致的调度开销。这为今后的并行计算和优化提供了

绝不收费

关注

10
点赞
踩
6

收藏

觉得还不错? 一键收藏
0
评论
矩阵乘法的OpenMP并行程序

同时，为了避免程序的多次间断的运行，导致程序运行时CPU的状态差，导致实验获取的结果误差较大，这里考虑，创建数组提前存储“1、2、4、8、16、32、64”，然后通过for循环将数值传入函数“double juzheng(int tp)”中，获取运行时间，并计算出对应5次实验的平均值。所以在高级语言的层面，除非通过编译器提供的方法或者本来的规定的语言的规则，否则，我们是无法轻易判断一个操作是原子操作的。根据对上述，表1和图2的分析可知，随着线程数的增加，矩阵乘法的运行时间逐渐下降，加速比逐渐上升。
复制链接

扫一扫