【高性能计算】OpenMP/MPI优化技巧

0.前言

笔者最近参与了并行计算相关的比赛，赛题主要内容就是把一份C源码的程序利用2个节点、每节点64个核进行优化（当然也包括使用其他优化手段，但主要的加速在于多线程/多进程）。新手上路，和队友在OpenMP/MPI折腾了不少时间，现在把一些优化的技巧记录在这里。

优化都不是绝对的，具体哪种方式适用于代码，还是要就事论事的吧。

1.OpenMP的使用方式

OpenMP最容易被想到的使用方式莫过于对循环进行加速：

#pragma omp parallel for
for (int i = 0; i < n; ++i) {
  ...
}

这当然是最简单易行的想法，但可能会因为循环当中存在数据竞争或者程序不是标准的循环形式而难以使用，同时，直接这样写并不清楚每个线程执行的细节，加速效果碎线程数/核数的增加可能并不明显。如果使用更原始的并行制导块，可以弥补这些缺点：

#pragma omp parallel
{
  int rank = omp_get_thread_num();
  ...
}

这样需要手动根据线程号划分任务，代码量大一些，但是能够清楚地知道每一个线程执行的细节。实践上，在较复杂的代码上似乎这种结构也更加泛用，效率更高。

OpenMP在多线程库中的一个突出特点就是很容易使用#pragma omp parallel for等将单线程改造成多线程，但为了追求更高的效率，可能也需要在逐线程细化分配任务与循环并行化两种风格之间灵活切换。

2.OpenMP+MPI

MPI中的rank对应每个进程，单个进程可以包括多个核，不同的进程可以在不同的节点上。而OpenMP中的每个线程号对应每个线程，一个线程不能对应多个核。一般线程数等于核数；线程数小于核数意味着有几个核上没有执行任务，线程数大于核数也能够执行，但这意味着有几个核需要不停地轮换执行多个线程的内容，实践上效率几乎总是更低。并且，这里所调用的核都是指同一节点上的核（即同一个物理设备）。OpenMP不能跨节点分配任务，MPI才可以。

对于超线程技术，一个核可以执行两个线程，那又是另外一回事。

以2个节点、每个节点64核（CPU）为例。在运行的时候，同样是128核满核运作，我们可以指定2个节点、128个进程、每进程1个线程，又或者2个节点、8个进程、每进程16个线程，等等。这可以在集群提交作业的sbatch或srun等命令中找到相应参数。

但是，按照这种想法，只要把每个核都指定成一个进程进行运作（比如2个节点、128个进程）不就不需要使用OpenMP了吗？

理论上确实是这样。但是，MPI相比OpenMP虽然有可以跨节点的优势，但是也有初始化慢的劣势。MPI_Init()函数执行所用的时间是随进程数增多而增加的。比赛时发现，如果采取上面这种做法，在数据规模较大的情况下MPI_Init()的时间也有将近主要计算花费时间的1/5，在数据规模小的情况下更是有主要计算时间的5、6倍。这对于小规模数据尤其不友好。而OpenMP就基本没有初始化的时间（OpenMP根本就没有初始化的函数嘛）。

因此，理论上最占优势的划分模式是把每节点作为一个进程，把每个节点下的所有核都作为线程处理。也就是上例中的指定2个节点、2个进程、每进程64个线程。实践上，这样做的运行速度确实也足够快。

3.OpenMP使用的位置

这里仅仅是一个经验之谈：OpenMP用在越外层，加速越明显。比如f()函数中有一个循环，在这个循环中调用了g()函数，而g()函数中又有循环；最后我们在main()函数中的一个循环里面反复调用了f()。（听起来有点绕）这时候如果打算加#pragma omp parallel for，就有很多位置可以加。经验上最好的做法是加在main()的循环外面，而不是g()的循环外面。不过，笔者也并不知道这是为什么。

4.数据竞争与内存优化

经常在并行的区域，会出现对同一个变量的写操作。假如这个变量是一个数组：

int a[N];
#pragma omp parallel
{
  ...
  a[]=...;
  ...
}

（这代码看起来很虚假，但是权当伪代码看看罢）

于是出现了数据冲突。如果发生冲突的是单个变量，那可以用private子句非常容易地解决。但是数组等比较复杂的数据结构似乎不能这样做（至少我不知道OpenMP有提供简单的解决办法）。为了保证代码的正确性，我们可能会选择手动完成private所做的事情，即给每个线程都分配一份相同的空间：

int size = omp_get_num_threads();
int* a = (int*)malloc(sizeof(int) * size * N);
#pragma omp parallel
{	
	int rank = omp_get_thread_num();
	...
	a[rank * N + i] = ...;
	...
}

然后再对每个线程的结果规约处理，得到正确的结果。

但是上面这种写法可能很慢。因为a[]的空间是在并行区域外申请的，也就是说，后续所有的线程都必须到同一个地方进行读写操作，而这个地方不一定是这个线程所在的CPU的内存，可能是在其他CPU上。这样读写速度会非常之慢；如果并行区域存在大量的读写操作，运行时间就会拖长不少。

稍微改动一下，情况就会好很多：

#pragma omp parallel
{
  int a[N];
  ...
  a[]=...;
  ...
}

这样，每个线程仍然会得到私有空间，但是是在当前的CPU上的空间。访存速度会快很多。但是，这样就必须把规约操作提前到并行区域完成，因为离开并行区域后a[]数组就不存在了。需要稍微动动脑子改造一下代码。

其实也就是一时不幸写成了最上面的版本……如果自己写成这样要意识到访存有问题就是了。