OpenMP Reductions

遇到竞争写情况怎么办

critical section

最简单的解决方案是通过声明一个critical部分来消除竞争。

double result = 0;
#pragma omp parallel num_threads(ndata)
{
  double local_result;
  int num = omp_get_thread_num();
  if (num==0)      local_result = f(x);
  else if (num==1) local_result = g(x);
  else if (num==2) local_result = h(x);
#pragma omp critical
  result += local_result;
}

double result = 0;
#pragma omp parallel
{
   double local_result;
#pragma omp for
   for (i=0; i<N; i++) {
    local_result = f(x,i);
#pragma omp critical
   result += local_result;
} // end of for loop
}

原子操作/加锁

性能是不好的，变串行了

#pragma omp atomic
 pi += sum;

static omp_lock_t lock;
void omp_init_lock(&lock)：初始化互斥器
void omp_destroy_lock(omp_lock*)：销毁互斥器
void omp_set_lock(omp_lock*)：获得互斥器
void omp_unset_lock(omp_lock*)：释放互斥器
void omp_test_lock(omp_lock*): 试图获得互斥器，如果获得成功则返回true，否则返回false

reduction clause 子句

将其添加到一个omp并行区域有如下效果。

• OpenMP将为每个线程制作一个reduction变量的副本，初始化为reduction操作的身份，例如$1$用于乘法。
• 然后，每个线程将其reduce到其本地变量中。
• 在并行区域结束时，本地结果被合并，再次使用reduction操作，合并到全局变量。

多个变量的情况

reduction(+:x,y,z)
reduction(+:array[:])

对于复杂结构体

如果代码过于复杂，还是建议复制全局变量来手工实现，最后再合并。

//错误示例
double result,local_results[3];
#pragma omp parallel
{
  int num = omp_get_thread_num();
  if (num==0)      local_results[num] = f(x)
  else if (num==1) local_results[num] = g(x)
  else if (num==2) local_results[num] = h(x)
}
result = local_results[0]+local_results[1]+local_results[2]

虽然上面这段代码是正确的，但它可能是低效的，因为有一个叫做虚假共享的现象。即使线程写到不同的变量，这些变量也可能在同一个缓存线上。这意味着核心将浪费大量的时间和带宽来更新对方的缓存线副本。

可以通过给每个线程提供自己的缓存线来防止错误的共享。

// 不是最好
double result,local_results[3][8];
#pragma omp parallel
{
  int num = omp_get_thread_num();
  if (num==0)      local_results[num][1] = f(x)
// et cetera
}

最好的方法给每个线程一个真正的局部变量，并在最后用一个critial部分对这些变量进行求和。

double result = 0;
#pragma omp parallel
{
  double local_result;
  local_result = .....
#pragam omp critical
  result += local_result;
}

默认的归约操作

Arithmetic reductions: \(+,*,-,\max,\min\)

Logical operator reductions in C: & && | || ^

归约变量的初始值

初始化值大多是不言而喻的，比如加法的0和乘法的1。对于min和max，它们分别是该类型的最大和最小可表示值。

用户自定义reduction的声明与使用

语法结构如下

#pragma omp declare reduction
    ( identifier : typelist : combiner )
    [initializer(initializer-expression)]

例子1: 取int最大

int mymax(int r,int n) {
// r is the already reduced value
// n is the new value
  int m;
  if (n>r) {
    m = n;
  } else {
    m = r;
  }
  return m;
}
#pragma omp declare reduction \
  (rwz:int:omp_out=mymax(omp_out,omp_in)) \
  initializer(omp_priv=INT_MIN)
  m = INT_MIN;
#pragma omp parallel for reduction(rwz:m)
  for (int idata=0; idata<ndata; idata++)
    m = mymax(m,data[idata]);

openmp减法归约浮点运算有精度损失

如何对vector归约

累加

#include <algorithm>
#include <vector>

#pragma omp declare reduction(vec_float_plus : std::vector<float> : \
                              std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<float>())) \
                    initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))

std::vector<float> res(n,0);
#pragma omp parallel for reduction(vec_float_plus : res)
for(size_t i=0; i<m; i++){
    res[...] += ...;
}

编辑：原始initializer很简单：initializer（omp_priv = omp_orig）。但是，如果原始副本没有全零，结果将是错误的。因此，我建议使用更复杂的initializer，它总是创建零元素向量。

求最大值

#pragma omp declare reduction(vec_double_max : std::vector<double> : \
                          std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), [](double a, double b) {return std::max(a,b);}))     \
                    initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))

#pragma omp parallel for reduction(vec_double_max:maxlab)
for( int i = 0; i < sz; i++ )
{
   maxlab[klabels[i]] = max(maxlab[klabels[i]],distlab[i]);
}

std::transform

在指定的范围内应用于给定的操作，并将结果存储在指定的另一个范围内。

需要进一步的研究学习

1. 对vector的归约

2. 泥菩萨:
   你这么改，开-g，在vtune里面看汇编

泥菩萨:
看有没有vmm指令

遇到的问题

暂无

开题缘由、总结、反思、吐槽~~

写IPCC发现：openmp没想象中简单，

参考文献

https://stackoverflow.com/questions/43168661/openmp-and-reduction-on-stdvector

https://pages.tacc.utexas.edu/~eijkhout/pcse/html/omp-reduction.html

http://www.cplusplus.com/forum/general/201500/