openmp理解与实践–常见问题解决
缘起
并行计算已经有些年没有碰过了。以前做CFD时，利用网格的分块考虑任务并行，所以主要利用mpi实现。
其实在cfd中mpi并行的确是比较适合，各个进程间维护自己的内存空间，利用有限的通信来实现交互，这是一种能够很自然理解的模型。
最近实现一个扑克赢率计算的小工具，在6人情况下，5个对手，当给出公共牌时，如果采用枚举遍历计算，
需要循环：
Binomial(45,2)*Binomial(43,2)*Binomial(41,2)*Binomial(39,2)*Binomial(37,2)
这样一个实现串型程序可能需要计算好几天，于是在想能否利用并行来进行加速。然而由于这个任务中的循环是强耦合在一起的，很难采用某种方式进行任务划分，所以利用mpi来实现多个进程的并行计算可能并不是非常合适，于是想能不能利用其它方式来实现，于是又想起openmp来，尽管我以前没有用过，但还有点映象，所以学习着来用。
显然openmp对于一定编程基础的人来说，似乎实现起来应该是很简单的事情。
我也是这么自认为算是能编程的，所以在csdn上收了几个帖子来看，的确不是非常复杂的。
做个简单例子，加上头文件，加上线程并行编译标识语句，编译运行，的确ok了。
但搬到赢率计算程序上，麻烦就来了，计算结果跟串行的不一致。怎么办？
我还是根据当前的理解快点解决掉问题，再看多几个帖子，看看别人的介绍和示例，然后再去尝试，
也用各种简单示例来实践，搞了将近一天还是存在问题。我想还是不能太偷懒，有时真的匀速则不达，还是老老实实去openmp官网，进行系统学习一下吧，不然靠示例是很难解决问题的，而且有的示例的实现与其解释是不对应的，有些概念的解释也是含混不清的。显然有些帖子是没有深入的理解openmp，也没有说到点上，这些就不谈了。
我的理解
我相信问题的关键在于自身对于openmp的理解，所以我在官网上下来几个tutorial来看，的确在循序渐进的讲解下，我应该是悟了。
不讲openmp的深层次原理，我总结openmp的理解如下：
-
利用openmp并行是线程的并行（线程包含在进程内，一个进程可以有多个进程。多个进程的并行可以把cpu的多个核用起来，多个线程的并行也可以
将cpu的多个核用起来） -
openmp是共享内存的，也就是说没有特别指定为线程私有的，所有的程序的东西都是共享的（而基于mpi的多进程并行，则是进程私有的，各个进程间不共享数据空间）。（这一点是解决问题的关键，特别是使用数组遇到一些稀奇的事情时，比如增加一句输出导致结果正确，而注释掉后结果由出错了。这个方面的错误大多与数组内存使用相关，在c、c++中常见，fortran也有。）
-
从“#pragma omp parallel” 开始程序进入并行区域，并行区域如果用{}包围，那么并行的代码就是包含在其中的内容，如果使用“#pragma omp parallel for”只针对for循环，那么并行的区域仅包含仅跟该标记的for循环的内部代码。在并行代码中如果使用“#pragma omp single”标记其后的一行代码或者用{}包围的代码，则表示这些代码仅在一个线程中运行，这个线程不一定是主线程。
-
为了程序运行结果正确，尽量减少各个线程运行的代码中的共享变量更新。可以使用private指令让各进程维护一个副本变量，但最好的方式还是使用临时变量而不要使用全局的共享变量，因为使用private维护变量副本时，有时你以为维护了，但其实没有，比如数组就容易出现这种问题，一般的变量可能在进入线程会设置成空，但数组不是自动将所有索引下的内存空间都处理的，所以如果不在线程内做赋值，很有可能会出错。
-
如果不可避免的要在各个线程内对全局的共享变量进行更新，那么要使用openmp提供的机制来进行处理，比如利用Reduction（规约）机制，利用critical机制，利用 atomic 机制。其中：
Reduction机制会在各个线程中创建规约列表中的变量副本，并使用规约列表中的运算符进行局部更新，最后将各局部副本合并到全局的共享变量中。
critical机制是利用一种互斥机制，同一时间只能在一个线程访问，称为临界区。 （Mutual exclusion: Only one thread at a time can enter a critical region.）。
atomic 机制也是利用一种互斥机制，同一时间只能一个线程访问，但只针对一个内存地址的更新，比如“X += tmp”。（Atomic provides mutual exclusion but only applies to the update of a memory location） -
数组是可以用private标记成副本的，但如前所述，在private中只记录了数组名，数组的大小信息其实线程是不清楚的，所以在线程执行代码中要进行赋值。
-
嵌套的在内部的循环一般不会自动并行化，而会采用串行运行。出发使用nested命令，我没有深入尝试，如有需要可以去了解。
-
openmp还有很多更复杂的机制，但我暂时用不到了，所以也没有深入去学习，有需要的可以去了解。
实践与代码讲解
先看一个简单的串行代码，我们主要是要对其中的三重循环进行并行优化，注意到最内层的循环与其他层是没有关系的，而第二层的遍历遍历k是与第一层的遍历变量j相关的。其计算输出的sumc结果为55。
#include <iostream>
#include <vector>
#include <string>
#include <time.h>
int main(int argc, const char * argv[]) {
clock_t start, end;
double cpu_time_used;
start = clock();
long long suma=0;
long long sumc=0;
for(long long j=0;j<7;j++)
{
long long sumb=0;
for(long long k=j+1;k<6;k++)
{
for(long long i=0;i<10000000;i++)
{
suma+=i;
}
sumb+=k;
}
sumc+=sumb;
}
end = clock();
std::cout<<"sumc="<<sumc << std::endl;
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
std::cout<<"time="<<cpu_time_used << std::endl;
std::cout << std::endl;
return 0;
}
要进行并行化，最简单的方式是采用“#pragma omp parallel for”标记，并使用规约来进行全局变量的更新，比如：
注意其中的两处修改：
#include <omp.h>-
#pragma omp parallel for reduction(+:sumc)
具体的作用前面已经解释过。
由于采用g++编译，所以命令加上选项 -fopenmp，比如：g++ test-p.cpp -o testpa5.exe -fopenmp -O3
在使用g++编译情况下，这一并行的代码比前面的串行代码运行还要慢一些，如果使用编译优化选项-O3，可能差异更明显可见，所以有时候并行并非一定会带来效率的提升，要看具体的问题来处理。
#include <iostream>
#include <vector>
#include <string>
#include <time.h>
#include <omp.h>
int main(int argc, const char * argv[]) {
clock_t start, end;
double cpu_time_used;
start = clock();
long long suma=0;
long long sumc=0;
#pragma omp parallel for reduction(+:sumc) //这句使用了规约操作没有问题，而且加速了很多
for(long long j=0;j<7;j++)
{
//printf("j=%d id=%d\n",j,omp_get_thread_num());
long long sumb=0;
for(long long k=j+1;k<6;k++)
{
for(long long i=0;i<200000000;i++)
{
suma+=i;
}
sumb+=k;
}
sumc+=sumb;
//printf("sumc=%d thread id=%d %d\n",sumc,omp_get_thread_num(),j);
}
end = clock();
std::cout<<"sumc="<<sumc << std::endl;
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
std::cout<<"time="<<cpu_time_used << std::endl;
std::cout << std::endl;
return 0;
}
第二种方式是不使用自动的并行，而自己来设计并行，如下面的代码所示
注意其中的修改：
-
#define MAX_THREADS 4主动定义线程的总数 -
omp_set_num_threads(MAX_THREADS);主动设定线程组，使并行区域根据这个设定运行，就是在未修改线程组前，后面的并行区域都按照最大线程数MAX_THREADS进行运行。 -
#pragma omp parallel private(j,tid)设置并行区域，后面跟着代码进入并行，并使用private指定了各线程需要维护副本的变量j，tid。
其实更好的处理是不使用全局的j和tid，而使用临时的j和tid，即在并行区域代码中定义这两个变量。 -
tid = omp_get_thread_num()各线程将其id赋值给tid，用于后面的并行任务分配。 -
for(j=tid;j<7;j+=MAX_THREADS)设计分配各个线程运行的代码，显然因为这样的设定，各个线程运行的任务就给定了。
对于线程0，那么运行的j循环包含j=0，j=4。
对于线程1，那么运行的j循环包含j=1，j=5。
对于线程2，那么运行的j循环包含j=2，j=6。
对于线程3，那么运行的j循环包含j=3。 - 全局更新的变量sumc，在不使用
#pragma omp critical或#pragma omp atomic时也没有出错，但最好还是用一下。
#include <iostream>
#include <vector>
#include <string>
#include <time.h>
#include <omp.h>
int main(int argc, const char * argv[]) {
clock_t start, end;
double cpu_time_used;
#define MAX_THREADS 4
start = clock();
long long suma=0;
long long sumc=0;
omp_set_num_threads(MAX_THREADS);
long long j;
int tid;
int numthreads;
#pragma omp parallel private(j,tid)
{
tid = omp_get_thread_num();//线程的id标记
for(j=tid;j<7;j+=MAX_THREADS)
{
//printf("j=%d thread id=%d\n",j,omp_get_thread_num());
long long sumb=0;
for(long long k=j+1;k<6;k++)
{
//printf("j=%d k=%d,id=%d\n",j,k,omp_get_thread_num());
for(long long i=0;i<200000000;i++)
{
suma+=i;
}
sumb+=k;
}
//printf("sumb=%d thread id=%d\n",sumb,omp_get_thread_num());
//#pragma omp critical
//#pragma omp atomic
sumc+=sumb;
}
//printf("sumc=%d thread id=%d\n",sumc,omp_get_thread_num());
}
end = clock();
std::cout<<"sumc="<<sumc << std::endl;
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
std::cout<<"time="<<cpu_time_used << std::endl;
std::cout << std::endl;
return 0;
}
下面再来看一个例子，循环的创建不同数量的线程组进行并行计算：
其中创建了6个线程组进行比较，分别使用1/2/3/4/5/6个线程进行计算。
#include <iostream>
#include <vector>
#include <string>
#include <time.h>
#include <omp.h>
int main(int argc, const char * argv[]) {
clock_t start, end;
double cpu_time_used;
start = clock();
for(int j1=1;j1<=6;j1++){ //c尝试多个不同线程数的线程组来运行比较
omp_set_num_threads(j1);
long long suma;
long long sumc;
suma=0;
sumc=0;
long long j;
int tid;
int numthreads;
#pragma omp parallel private(j,tid)
{
tid = omp_get_thread_num();//线程的id标记
numthreads = omp_get_num_threads();
//printf("number threads=%d id=%d\n",numthreads,omp_get_thread_num());
for(j=tid;j<7;j+=numthreads)
{
//printf("j=%d thread id=%d %d\n",j,omp_get_thread_num(),j1);
long long sumb=0;
for(long long k=j+1;k<6;k++)
{
for(long long i=0;i<200000000;i++)
{
suma+=i;
}
sumb+=k;
}
//printf("sumb=%d thread id=%d %d\n",sumb,omp_get_thread_num(),j1);
#pragma omp atomic
sumc+=sumb;
}
#pragma omp single
{
printf("sumc=%d thread id=%d %d\n",sumc,omp_get_thread_num(),j1);
end = clock();
std::cout<<"sumc="<<sumc << std::endl;
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
std::cout<<"time="<<cpu_time_used <<"\n"<< std::endl;
}
}
}
std::cout << std::endl;
return 0;
}
最后来看以下全局定义的数组的副本维护情况
注意到：
对数组j进行了副本维护，但显然在超出数据范围后，各线程仍然能够访问和使用。
所以在使用数组时要特别注意。这也是c++、fortran这类编译语言的特点，不对数组的范围进行维护，一切都要代码编写者自己注意。
比如：
#include <iostream>
#include <vector>
#include <string>
#include <omp.h>
int main(int argc, const char * argv[]) {
int k = 20;
int j[10];
j[0] = 20;
#pragma omp parallel for private(k,j)
for ( k=0; k < 11; k++)
{
//printf("k->%d\t", k);
j[k] = k*k;
printf("j[%d]->%d\t",k,j[k]);
}
printf("\n");
for ( k=0; k < 11; k++)
{
printf("j[%d]->%d\t",k,j[k]);
}
printf("\n");
#pragma omp parallel for
for ( k=0; k < 11; k++)
{
//printf("k->%d\t", k);
j[k] = k*k;
printf("j[%d]->%d\t",k,j[k]);
}
printf("\n");
for ( k=0; k < 11; k++)
{
printf("j[%d]->%d\t",k,j[k]);
}
printf("\n");
return 0;
}
小结
针对直接的需求问题，在学习基础上给出了openmp的理解，并进行了代码实践，应该说解决了openmp的初步使用问题。
因为没有需求的原因，也没有继续深入下去应用更为复杂的功能。简单其实也是我追求的，尽可能使用一些简单的东西，对于后期的理解来说是更有好处的，因为很有可能用完这段之后很久都不会再用了，发出来下来当做记录吧。
参考：
- https://www.openmp.org/resources/tutorials-articles/
- Hands-On Introduction to OpenMP, Mattson and Meadows, from SC08 (Austin)
- https://blog.csdn.net/qq_40379678/article/details/107788716?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-5&spm=1001.2101.3001.4242
- 其它一些csdn不少文章，也浏览过，没有一一记录下来，有兴趣自行搜索。