并行计算(一)——OpenMP

一、简介

OpenMP是一种用于共享内存并行系统的多线程库，其支持C/C++、Fortran，并且目前大多数常用编译器，如VS内置编译器、gcc、icc等都提供了openmp的相关支持，以gcc为例编译时只需要添加-fopenmp选项即可完成OpenMP代码的编译。OpenMP中包含了一套编译器伪指令、运行时函数和一些环境变量。其通过对串行代码的很少的修改就可以实现串行代码的并行化(不过要想得到更好性能依旧需要仔细的设计)，并且可以*控制编译器是否忽略代码中的并行段，从而实现一套代码既可以串行亦可以并行，十分方便。

OpenMP提供的是多线程并行适用于多核CPU，不同线程间的数据可以直接共享，无需进行消息传递。利用OpenMP将任务分为多个子块，分发给不同的CPU核心从而提高多核CPU的利用率，加速程序的运行。

二、简单使用

OpenMP测试就必须先说一下环境了，测试机器是双核四线程Intel的CPU，OpenMP默认使用四个线程进行并行。当然我们也可以通过在并行代码块前使用omp_set_num_threads(int)函数来手动指定使用的线程数量。

矢量相加

矢量相加它又来了

#include <iostream>
using namespace std;
#define N 20000000
int main(){
    double *x,*y,*z;
    x=(double*)malloc(sizeof(double)*N);
    y=(double*)malloc(sizeof(double)*N);
    z=(double*)malloc(sizeof(double)*N);
    for(int i=0;i<N;i++) z[i]=x[i]+y[i];
    free(x);
    free(y);
    free(z);
    return 1;
}

这是一个简单的矢量相加程序，这个程序的用时我之前的文章里已经测过了，这里不再测试，后面结果展示中这个串行程序的时间都是从前面文章里拷贝的了。接下来我们用OpenMP并行它。使用OpenMP加速这种简单循环很容易，只需要在for循环前面添加一个#pragma omp parallel for即可。

#include <iostream>
#include "omp.h"
using namespace std;
#define N 20000000
int main()
{
	double *x,*y,*z;
	x=(double*)malloc(sizeof(double)*N);
	y=(double*)malloc(sizeof(double)*N);
	z=(double*)malloc(sizeof(double)*N);
	#pragma omp parallel for//自动并行后面的for循环
		for(int i=0;i<N;i++) z[i]=x[i]+y[i];
	cout<<z[0]<<'\n';
	free(x);
	free(y);
	free(z);
	return 1;
}

测试一下速度，这里就不测关闭优化时的性能了，gcc O3优化
未使用OpenMP

real	0m0.112s
user	0m0.040s
sys		0m0.072s

使用OpenMP

real	0m0.148s
user	0m0.095s
sys		0m0.251s

非常遗憾，OpenMP非但没有提高性能反而使性能有所下降。接下来我们就不靠OpenMP自带的循环优化来并行，而使用手动并行看一下效果。

想要手动并行需要使用的主要有#pragma omp parallel(用于告知编译器接下来的代码块并行执行)，omp_get_thread_num函数(获得当前线程的线程号)，omp_get_num_threads函数(获得并行的总线程数)。

#include <iostream>
#include "omp.h"
using namespace std;
#define N 20000000
int main()
{
	double *x,*y,*z;
	x=(double*)malloc(sizeof(double)*N);
	y=(double*)malloc(sizeof(double)*N);
	z=(double*)malloc(sizeof(double)*N);
	#pragma omp parallel//接下来的代码块为并行代码
    {
        int id=omp_get_thread_num(),//获取当前所在线程的线程编号
        	n=omp_get_num_threads();//获取并行的总线程数量
		for(int i=id*N/n;i<(id+1)*N/n;i++) z[i]=x[i]+y[i];
    }	
    cout<<z[0]<<'\n';
	free(x);
	free(y);
	free(z);
	return 1;
}

同样使用O3优化用时为

real	0m0.065s
user	0m0.053s
sys		0m0.137s

可以看到这时的用时就少于串行的用时了，大概快了40%，为什么4个线程没有达到4倍的性能。这是因为机器本身是双核的，四线程是通过一些其他技术多模拟出来了两个线程，理论上四线程比双线程的性能提升并不能达到真四核四线程的效果。另外这个任务本身耗时就不多，而多线程是有额外开销的，所以在一定程度上会导致四线程并行的性能达不到4倍。至于自动并行变慢原因，目前我也不是很明白，初步估计可能是自动并行的方式不对。

随机数初始化

这是我以前踩过的一个坑，这里讲一下。

串行代码如下

#include <iostream>
#include <cstdlib>
using namespace std;
#define N 80000000
int main()
{
	double *x;
	x=(double*)malloc(sizeof(double)*N);
    for(int i=0;i<N;i++) x[i]=double(rand())/RAND_MAX;
    cout<<x[0]<<'\n';
	free(x);
	return 1;
}

这个代码主要是完成数组的随机初始化。接下来手动并行它。

#include <iostream>
#include <cstdlib>
#include "omp.h"
using namespace std;
#define N 80000000
int main()
{
	double *x;
	x=(double*)malloc(sizeof(double)*N);
    //for(int i=0;i<N;i++) x[i]=double(rand())/RAND_MAX;
    #pragma omp parallel
    {
        int id=omp_get_thread_num(),n=omp_get_num_threads();
		for(int i=id*N/n;i<(id+1)*N/n;i++) x[i]=double(rand())/RAND_MAX;
    }
    cout<<x[0]<<'\n';
	free(x);
	return 1;
}

gcc O3优化，测试结果如下：
串行版本

real	0m0.935s
user	0m0.802s
sys		0m0.132s

并行版本

real	0m11.101s
user	0m14.880s
sys		0m23.758s

我们看到用时直接多了一个数量级，我自己都觉得是不是测反了，但是没错确实并行版的效率更低了，这是为什么呢？回答这个问题之前我们先自己做一个伪随机数生成器。

unsigned rand(int &seed)
{
	seed = seed * 1103515245 + 12345;
    return((unsigned)(seed/65536) % 32768);
}

然后我们用自己的伪随机数生成器再做一次并行的随机数初始化

#include <iostream>
#include <cstdlib>
#include "omp.h"
using namespace std;
#define N 80000000
unsigned myrand(unsigned &seed)
{
	seed = seed * 1103515245 + 12345;
    return((unsigned) seed% 32768);
}
int main()
{
	double *x;
    unsigned *seed;
	x=(double*)malloc(sizeof(double)*N);
    //for(int i=0;i<N;i++) x[i]=double(rand())/RAND_MAX;
    #pragma omp parallel
    {
        int id=omp_get_thread_num(),n=omp_get_num_threads();
        if(id==0) {
            seed=(unsigned*)malloc(sizeof(unsigned)*n);
            for(int i=0;i<n;i++) seed[i]=rand();
        }
        #pragma omp barrier//所有线程在此处同步后再向下执行
        for(int i=id*N/n;i<(id+1)*N/n;i++) x[i]=double(myrand(seed[id]))/32768;
    }
    cout<<x[0]<<'\n';
	free(x);
	free(seed);
	return 1;
}

这个程序用时是

real	0m0.194s
user	0m0.215s
sys		0m0.266s

这一下性能明显提高了，而且达到之前性能的4倍。这又是为什么呢？这就是前面说过的踩的坑，当初这个问题困扰了我近两个月，使我一度怀疑OpenMP就是个废物，根本没有用。直到后来我了解了伪随机数的生成算法我才发现问题的关键。标准库中的rand函数是使用一个静态变量来保存随机数种子的，每调用一次rand函数，函数就会修改一次这个种子，从而使得每一次生成的数字都不相同。那这和之前那个效率极低的并行程序又有什么关系呢？正是因为每一次调用rand都要修改种子，所以rand执行时就既要读一次种子又要写一次种子，而不同线程又是共有一个rand函数，也就是说共用一个种子，那么必然会出现不同线程同时读写种子的情况。这个时候为了避免读写冲突OpenMP会自动给种子加锁，保证不同线程不会同时读写种子，这就导致了rand函数实际执行中并没有并行，而是总在等待其他线程完成rand之后再去执行rand，程序依然还是串行的，而由于锁的额外开销就导致了并行版本的程序性能严重下降。当我们使用自己的随机数生成器计算时，由于给每个线程各自分配了一个种子，因此避免了不同线程共用一个种子时产生的冲突，所以程序可以真正的并行执行，性能也提高了。