矩阵乘法的OpenMP实现及性能分析.doc

一. 实验目的
1) 用OpenMP实现最基本的数值算法“矩阵乘法”
2) 掌握for编译制导语句
3) 对并行程序进行简单的性能
二. 实验环境
1) 硬件环境：32核CPU、32G内存计算机；
2) 软件环境：Linux、Win2003、GCC、MPICH、VS2008；
4) Windows登录方式：通过远程桌面连接192.168.150.197，用户名和初始密码都是自己的学号。
三. 实验内容
1. 用OpenMP编写两个n阶的方阵a和b的相乘程序，结果存放在方阵c中，其中乘法用for编译制导语句实现并行化操作，并调节for编译制导中schedule的参数，使得执行时间最短，写出代码。
方阵a和b的初始值如下：

输入：
方阵的阶n、并行域的线程数
输出：
c中所有元素之和、程序的执行时间
提示：
a,b,c的元素定义为int型，c中所有元素之各定义为long long型。
Windows计时:
用<time.h>中的clock_t clock( void )函数得到当前程序执行的时间
Linux计时:
#include <sys/time.h>
timeval start,end;
gettimeofday(&start,NULL);
gettimeofday(&end,NULL);
cout<<"execution time:"<< (end.tv_sec-start.tv_sec)+(double)(end.tv_usec-start.tv_usec)/ 1000000<<"seconds" <<endl;

答：

在windows下使用Microsofe Visual Studio编程，源代码如下：

#include <omp.h>
#include <stdio.h>
#include <time.h>

#define NN 2000

int a[NN][NN], b[NN][NN];
long long c[NN][NN];

void solve(int n, int num_thread)
{
int i, j, t, k, time;
clock_t startTime, endTime;
long long sum;

omp_set_num_threads(num_thread);

for(i=0;i<n;i++)//对矩阵a和矩阵b进行初始化
{
t=i+1;
for(j=0;j<n;j++)
{
a[i][j]=t++;
b[i][j]=1;
}
}
startTime=clock();
sum=0;

#pragma omp parallel shared(a,b,c) private(i,j,k)
{
#pragma omp for schedule(dynamic)

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
c[i][j]=0;
for(k=0;k<n;k++)
{
c[i][j]+=a[i][k]*b[k][j];
}
}
}
}

for(i=0;i<n;i++)for(j=0;j<n;j++) sum+=c[i][j];

endTime=clock();
time=endTime-startTime;

printf("sum=%lld time=%dms\n",sum,time);
}

int main()
{
int n, num_thread;

while(scanf("%d%d",&n,&num_thread)!=EOF)
{
solve(n,num_thread);
}
return 0;
}

2. 分析矩阵相乘程序的执行时间、加速比和效率：方阵阶固定为1000，节点数分别取1、2、4、8、16和32时，为减少误差，每项实验进行5次，取平均值作为实验结果。
答： 串行执行时程序的执行时间为：T = 15.062s
加速比=顺序执行时间/并行执行时间
效率=加速比/节点数
表1 不同节点数下程序的执行时间（秒）
节点数
实验结果
1
2
4
8
16
32
第1次
16.640
8.172
4.078
2.125
1.093
0.594
第2次
16.422
8.156
4.172
2.141
1.078
0.578
第3次
16.406
8.266
4.078
2.125
1.094
0.563
第4次
16.781
8.172
4.079
2.109
1.094
0.563
第5次
16.422
8.171
4.078
2.125
1.093
0.578
平均值
16.5342
8.1874
4.0970
2.1250
1.0904
0.5752

图1 不同节点数下程序的执行时间

图2 不同节点数下程序的加速比

图3 不同节点数下程序的效率

执行时间的分析：
随着节点数的增加，程序的执行时间减少，大概可以从结果中得出，随着节点书的增加一倍，执行时间减少一半
加速比的分析：
随着节点数的增加，程序的加速比增加，大概可以从结果中得出，随着节点书的增加一倍，加速相应的增加接近一倍
效率的分析：
随着节点数的增加，程序的效率逐渐减少

3. 分析矩阵相乘程序的问题规模与效率的关系：固定节点数为4，让方阵阶从200到1600之间变化，每隔100取一个值。（为了减少时间，每项实验可只执行1次）
答：
表2 相同节点数下不同问题规模程序的执行时间与效率
方阵阶数
并行执
行时间
串行执
行时间
效率
200
0.015
0.047
0.783333
300
0.016
0.109
1.703125
400
0.063
0.297
1.178571
500
0.156
0.657
1.052885
600
0.406
1.64
1.009852
700
0.907
3.578
0.986218
800
1.609
6.36
0.988191
900
2.578
10.109
0.980314
1000
3.812
14.891
0.976587
1100
5.39
21.032
0.97551
1200
7.344
28.734
0.978145
1300
9.688
37.937
0.978969
1400
12.422
48.64
0.978908
1500
15.656
60.938
0.973077
1600
19.234
74.829
0.972614

图3.1 不同问题规模下程序的效率
问题规模与效率的关系分析：

随着问题规模的增加，程序的效率趋于稳定，但是略微有点下降。

嵌套循环中,如果外层循环迭代次数较少时,如果将来CPU核数增加到一定程度时,创建的线程数将可能小于CPU核数。另外如果内层循环存在负载平衡的情况下，很难调度外层循环使之达到负载平衡。
下面以矩阵乘法作为例子来讲述如何将嵌套循环并行化，以满足上述扩展性和负载平衡需求。
一个串行的矩阵乘法的函数代码如下：
/**矩阵串行乘法函数
@param int*a -指向要相乘的第个矩阵的指针
@param int row_a -矩阵a的行数
@param int col_a -矩阵a的列数
@param int *b –指向要想成的第个矩阵的指针
@param int row_b -矩阵b的行数
@param int col_b -矩阵b的列数
@param int *c -计算结果的矩阵的指针
@param int c_size -矩阵c的空间大小（总元素个数）
@return void – 无

*/
void Martrix_Multiply(int *a, int row_a,int col_a,
int*b,int row_b,int col_b,
int*c,int c_size)
{
If(col_a!=row_b||c_size<row_a*col_b)
{
return;
}
int i,j,k;
//#pragma omp for private(i,j,k)
for(i = 0;i<row_a;i++)
{
int row_i=i*col_a;
int row_c=i*col_b;
for(j=0;j<col_b;j++)
{
c[row_c+j]=0;
for(k=0;k<row_b;k++)
{
c[row_c+j]+=a[row_i+k]*b[k*col_b+j];
}
}
}
}
如果在外层循环前面加上OpenMP的for语句时，它就变成了一个并行的矩阵乘法函数，但是这样简单地将其并行化显然无法满足前面所述的扩展性需求。
其实可以采用一个简单地方法将最外层循环和第2层循环合并成一个循环，下面便是采用合并循环后的并行实现。
void Parallel_Matrix_Multiply(int *a,int row_a,int col_a,
int *b,int row_b,int col_b,
int *c,int c_size)
{
If(col_a!=row_b)
{
return;
}
int i,j,k;
int index;
int border=row_a*col_b;

i=0;
j=0;
//#pragma omp parallel private(i,j,k) num_threads(dtn(border,1))
for(index = 0;index<border;index++)
{
i=index/col_b;
j=index%col_b;

int row_i=i*col_a;
int row_c=i*col_b;

c[row_c+j]=0;
for(k=0;k<row_b;k++)
{
c[row_c+j]+=a[row_i+k]*b[k*col_b+j];
}
}
}

从上面代码可以看出，合并后的循环便捷border=row_a*col_b；即等于原来的两个循环边界之积，然后再循环中计算出原来的外层循环和第2层循环的迭代变量i和j，采用除法和取余来求出i和j的值。
需要值得注意的是，上面求i和j的值必须要保证循环迭代的独立性，即不能有循环迭代间的依赖关系。不能讲求i和j的值得过程优化成如下的形式
if(j==col_b)
{
j=0;
i++;
}
//.......此处代表实际的矩阵乘法代码
j++；
上面这种优化，省去了除法，效率高，但是只能在串行代码中使用，因为它存在循环迭代间的依赖关系，无法将其正确地并行