OpenMP π值估计

π值估计

1.数学背景

image-20230109150836877

我们能够在串行代码下实行这个公式：

1
2
3
4
5
6
7
double factor = 1.0;
double sum = 0.0;
for(k=0 ; k < n; k++){
    sum += factor /(2*k+1);
    factor = - factor;
}
pi_approx = 4.0*sum;

2.OpenMP并行化

为了使用OpenMP来并行化，可以首先倾向于这样做：

1
2
3
4
5
6
7
8
1 double factor = 1.0;
2 double sum = 0.0;
3 #pragma omp parallel for num_threads(thread_count) reduction(+:sum)  //对sum进行求和归约
4 for(k=0 ; k < n; k++){
5     sum += factor /(2*k+1);
6     factor = - factor;
7 }
8 pi_approx = 4.0*sum;

然而，第k次迭代中对第6行的factor的更新和接下来的第k + 1次迭代中对第5行的sum的累加是一个循环依赖（数据依赖）。如果第k次迭代被分配·到一个线程，而第k + 1次迭代被分配给另一个线程，则我们不能保证第6行中factor的值是正确的。

在这种情况下我们能通过检查系数来解决这个问题：

image-20230109161031989

可以看到：在第k次迭代，factor的值应该是image-20230109161206995。如果k是偶数，那么值是+1；如果k是奇数，值是-1。

3.消除循环依赖（数据依赖）

因此，如果将下述代码：

1
2
sum += factor /(2*k+1);
factor = - factor;

替换为：

1
2
3
4
5
if(k % 2 ==0)			//通过奇偶性，来独立factor消除循环依赖（数据依赖）
    factor = 1.0;
else
    factor = -1.0;
sum += factor/(2*k+1);

这样就消除了循环依赖（数据依赖）。

4.作用域

​ 然而，事情仍然不是完全正确的。如果在我们的系统上使用两个线程运行程序，并设n=1000，那么结果仍然是错误的。例如，

image-20230109164120205

另一方面，如果只有一个线程运行程序，我们总是得到：

image-20230109164204808

为什么会有这种错误。在一个已经被parallel for指令并行化的块中，缺省情况下任何在循环前声明的变量（唯一的例外是循环变量）在线程间都是共享的。因此factor被共享（被所有线程所共享）。例如，线程0可能会给他赋值1，但在它能用这个值更新sum前，线程1可能又给他赋值为-1了。因此，除了消除计算factor时的循环依赖（数据依赖）外，我们还需要保证每个线程有它自己的factor副本，就是说，为了使代码正确，我们需要保证factor有私有作用域（简单来说就是保证当前线程的factor的值不能被其他线程修改，也只有当前线程能更新和使用factor）。通过添加一个private子句到parallel指令中来实现这一目标。

1
2
3
4
5
6
7
8
9
10
11
 double factor = 1.0;
 double sum = 0.0;
#pragma omp parallel for num_threads(thread_count) reduction(+:sum) private(factor)  //对sum进行求和归约
 for(k=0 ; k < n; k++){
     if(k % 2 ==0)			//通过奇偶性，来独立factor消除循环依赖（数据依赖）
    	factor = 1.0;
	else
    	factor = -1.0;
     sum += factor /(2*k+1);
 }
 pi_approx = 4.0*sum;

在private子句内列举的变量，在每个线程上都有一个私有副本被创建。因此，在我们的例子中，thread_count个线程中的每一个都有它自己的factor变量的副本，因此一个线程对factor的更新不会影响另一个线程的factor值。

​ 要记住的重要的一点是，一个有私有作用域的变量的值在parallel块或者parallel for块的开始处是未指定的。它的值在parallel或parallel for块完成之后也是未指定的。例如，下列代码中的第一个printf语句的输出是非确定的，因为在它被现实初始化之前就打印了私有变量x。类似地，最终的printf输出也是非确定的，因为他在parallel块完成之后打印x。

1
2
3
4
5
6
7
8
9
int x = 5;
#pragma omp parallel num_threads(thread_count) private(x)
{
    int my_rank = omp_get_thread_num();
    printf("Thread %d > before initialization,x = %d\n",myrank,x);
    x = 2*my_rank + 2;
    printf("Thread %d > after initialization,x = %d\n",my_rank,x);
}
printf("After parallel block, x = %d\n",x);

5.关于作用域的更多问题

​ 关于变量factor的问题是常见问题中的一个。我们通常需要考虑在parallel块或parallel for块中的每个变量的作用域。因此，与其让OpenMP决定每个变量的作用域，还不如让程序员明确块中每个变量的作用域。事实上，OpenMP提供了一个子句default，该子句显示地要求我们这样做。如果我们添加子句：

1
default(none)

到parallel或parallel for指令中，那么编译器将要求我们明确在这个块中使用的每个变量和已经在块之外声明的变量的作用域。（在一个块中声明的变量都是私有的，因为它们会被分配给线程的栈。）

​ 例如，使用一个default（none）子句，对π的计算将如下所示。

1
2
3
4
5
6
7
8
#pragma omp parallel for num_threads(thread_count) default(none) reduction(+:sum) private(k,factor)  //对sum进行求和归约
 for(k=0 ; k < n; k++){
     if(k % 2 ==0)			//通过奇偶性，来独立factor消除循环依赖（数据依赖）
    	factor = 1.0;
	else
    	factor = -1.0;
     sum += factor /(2*k+1);
 }

在这个例子中，我们在for循环中使用4个变量。由于default子句，我们需要明确每个变量的作用域。正如我们已经注意到的，sum是一个归约变量（同时拥有私有和共享作用域的属性）。我们也已经注意到factor和循环变量中k应该有私有作用域。从未在parallel或parallel for块中更新的变量，如这个例子中的n，能够被安全的共享。与私有变量不同，共享变量在块内具有在parallel或parallel for块之前的值，在块之后的值与块内的最后一个值相同。因此，如果n在块之前被初始化为1000，则在parallel for语句中他将保持这个值。因为在for循环中值没有改变，所有在循环结束后它将保持这个值。

6.总结

  1. 分析数学背景，解决循环依赖（数据依赖）
  2. 判断变量的作用域

7.参考资料

并行程序导论 （美）Peter S.Pacheco

HPC
OpenMP

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！