矩阵相乘的OpenMP实现
OpenMP简介
OpenMP是基于共享内存的编程,不同于MPI,因为是共享内存,所以它不需要将计算结果丢来丢去共享。事实上,我们可以用很少的编译指导语句实现并行,而不需要关心底层的操作。
举个简单的例子,如果程序有for循环,我们只要在for循环前加了一句#pragma omp parallel for
,在循环语句前后不相关的情况下,对for循环就实现了并行。如下两种写法:
#pragma omp parallel for
for()
或者
#pragma omp parallel
{
#pragma omp for
for()
}
在OpenMP中,我们常用omp_get_num_procs()
函数来获取可用处理器个数,用omp_get_thread_num()
函数来获得每个线程的ID,为了使用这两个函数,我们需要include <omp.h>
。另外,omp_get_num_threads()
返回当前并行区域中的活动线程个数,如果在并行区域外部调用,返回1。omp_set_num_threads(5)
设置进入并行区域时,将要创建的线程个数为5个。在指导语句中直接设置核心数也是可以的:#pragma omp parallel num_threads(6)
。
OpenMP的规约操作以形如#pragma omp parallel for reduction(+:sum)
,做编译指导,它的意思是告诉编译器:下面的for循环你要分成多个线程跑,但每个线程都要保存变量sum的拷贝,循环结束后,所有线程把自己的sum累加起来作为最后的输出,避免了写入脏数据。
#pragma omp critical
语句告诉我们各个线程并行执行语句,但当你们执行到critical里面时,要注意有没有其他线程正在里面执行,如果有的话,要等其他线程执行完再进去执行。这样就避免了数据相关的问题,但显而易见,它的执行速度会变低,因为可能存在线程等待的情况。
sections语句的用法如下:
#pragma omp parallel sections
#pragma omp sections
它告诉计算机,sections里面的内容要并行执行,具体分工上,每个线程执行其中的一个section,如果section数大于线程数,那么就等某线程执行完它的section后,再继续执行剩下的section。
改造程序:矩阵乘法OpenMP并行
先给出程序如下:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#define m 2000
float a[m][m];
float b[m][m];
float c[m][m];
int main() {
clock_t start, end;
int i, j, k;
// srand((unsigned)time(NULL));//根据时间来产生随机数种子
srand(0);
for (i = 0; i<m; i++) {
for (j = 0; j<m; j++) {
a[i][j] = (float)rand() / (RAND_MAX);//产生[0,1]随机数
b[i][j] = (float)rand() / (RAND_MAX);
c[i][j] = 0;//初始化c矩阵用来存结果
#if m <= 1 //打印一下呗
printf("%f ", a[i][j]);
if (j == m - 1) {
printf("\n");
}
#endif
}
}
// int gettimeofday(struct timeval*tv, struct timezone *tz);
omp_set_num_threads(5);
int pnum = omp_get_num_procs();
fprintf(stderr, "Thread_pnum = %d\n", pnum);
// gettimeofday(&st, NULL);
start = clock();
#pragma omp parallel shared(a,b,c) private(j,k)
{
#pragma omp for //schedule(dynamic)
for (i = 0; i<m; i++) {
// printf("i = %d,thread_num:%d\n", i, omp_get_thread_num());
// fprintf(stderr, "Thread# %d: pnum = %d\n", omp_get_thread_num(), pnum);
for (j = 0; j<m; j++) {
// printf("j = %d,thread_num:%d\n\n",j,omp_get_thread_num());
for (k = 0; k<m; k++) {
c[i][j] = c[i][j] + a[i][k] * b[k][j];
//printf("%d%d%d ",i,j,k);
// printf("thread_num:%d\n",omp_get_thread_num);
}
}
}
}
end = clock();
// printf("time=%d\n", difftime(end, start));
// printf("matrix multiply time: %0.6lf sec\n", et.tv_sec + et.tv_usec*1e-6 - st.tv_sec - st.tv_usec*1e-6);//tv_sec表示秒, tv_usec表示微妙
printf("matrix multiply time: %0.6lf sec\n", ((double)end - start)/CLK_TCK);
getchar();
return 0;
}
注意到,并行程序块中的j、k两个循环指标要设置为私有,只有这样,各个进程执行任务才能完整且不互相干扰。
结果如下:当使用一个线程时用时69.246s,两个线程44.594s,3个线程39.817s,4个线程38.333s,再往后,基本上就没什么提升了。
我的CPU核心是4个,所以4个进程之后速度提升就不明显了,所以我们一般设置的线程数不超过CPU核心数。因为我的核心数目不够,增加线程只会带来OpenMP共享内存的额外开销,得不偿失。这里我只对第一重循环做了并行,事实上,第二重循环表示矩阵和向量乘法,也是可以并行。甚至第三重循环,使用规约求和也是可以并行的。但是,由于我的计算机的核心数及其有限,做太多的并行化,并没有太多意义。
数据相关简单讨论
并不是所有的循环,所有的代码块都可以做并行,当做并行分割后,存在数据依赖,或者由于写写冲突导致脏数据的话,那么就不能直接简单地加一句代码来实现并行。因为我们谁也不知道哪个线程跑得快,不知道线程跑的速度不同是否会改变原来的串行逻辑关系。
数据相关包括,输出相关、流相关、反相关等,输出相关可表述为,当多个线程并行执行时,有可能多个线程同时对某变量进行了读写操作,从而导致不可预知的结果。流相关和反相关表示计算结果上的依赖。
例题:列举出下面循环中的流相关、输出相关、反相关
for(i=1;i¡10; i++){
s1: a[i]=b[i]+a[i];
s2: c[i+1]=a[i]+d[i];
s3: a[i-1]=2*b[i];
s4: b[i+1]=2*b[i];
}
容易看得出来,输出相关有:s1和s3,流相关有:s1和s2、s2和s3,反相关有:s1和s3、s1和s4、s3和s4。注意,要考虑循环依赖,需要将此过程多展开一层。