文章介绍了如何使用OpenMP进行共享内存编程，通过示例展示了并行计算梯形积分法，讨论了线程安全、数据依赖和归约子句在多线程环境中的应用，以及如何处理竞争条件。同时，提到了并行化冒泡排序和奇偶变换排序的挑战。
用OpenMP进行共享内存编程
Hello world
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
void Hello(void);
int main(int argc,char* argv[]) {
int thread_count = strtol(argv[1], NULL, 10);/*不使用第二个参数，因此只是传入一个NULL(空)指针*/
# pragma omp parallel num_threads(thread_count)
Hello();
return 0;
}
void Hello(void) {
int my_rank = omp_get_thread_num();
int thread_count = omp_get_num_threads();
printf("Hello from thread %d of %d\n", my_rank, thread_count);
}
- strtol 函数说明
long strtol(
const char* number p, /* in */
char** endp, /* out */
int base /* in */);
)
第一个参数是一个字符串，最后一个参数是字符串所表示的数的基数
Hello world升级版
- 输入：创建的线程数和需要打印的内容
- 输出：由各线程打印出内容
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
void Hello(char* x);
int main(int argc, char* argv[])
{
int thread_count = strtol(argv[1], NULL, 10);
char* x = argv[2];
# pragma omp parallel num_threads(thread_count)
Hello(x);
return 0;
}
void Hello(char* x)
{
int my_rank = omp_get_thread_num();
int thread_count = omp_get_num_threads();
printf("print %s %d of %d\n", x, my_rank, thread_count);
}
在控制台调用生成的exe文件可得到如下输出结果：
梯形积分法
估计
y
=
f
(
x
)
y=f(x)
y=f(x)与垂直线
x
=
a
、
x
=
b
x=a、x=b
x=a、x=b和
x
x
x轴所围成的区域的面积，方法是将区间
[
a
,
b
]
[a,b]
[a,b]分成
n
n
n个子区间并在每个子区间上使用一个梯形的面积来近似估计区域的面积。
每个子区间取同样的宽度 h h h，并且定义 h = ( b − a ) / n ， x i = a + i h ， i = 0 , 1 , . . . , n h=(b-a)/n，x_i=a+ih，i=0,1,...,n h=(b−a)/n，xi=a+ih，i=0,1,...,n，近似值为： h [ f ( x 0 ) / 2 + f ( x 1 ) + f ( x 2 ) + . . . + f ( x n − 1 ) + f ( x n ) / 2 ] h[f(x_0)/2+f(x_1)+f(x_2)+...+f(x_{n-1})+f(x_n)/2] h[f(x0)/2+f(x1)+f(x2)+...+f(xn−1)+f(xn)/2]
可以使用下列代码实现梯形积分的串行算法：
h = (b-a)/n;
approx = (f(a)+f(b))/2.0;
for (i = 1; i <= n-1; i++){
x_i = a+i*h;
approx += f(x_i);
}
approx = h*approx;
方法1
假设梯形数量远大于核数量，通过给每个线程分配连续的梯形块，能将区间
[
a
,
b
]
[a,b]
[a,b]划分成更大的子区间，每个线程对它的子区间简单地应用串行梯形积分法。
#include <iostream>
using namespace std;
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
void Trap(double a, double b, int n,double* global_result_p);
int main(int argc, char* argv[]) {
double global_result = 0.0;
double a, b;
int n;
int thread_count;
thread_count = strtol(argv[1], NULL, 10);
printf("Enter a,b and n\n");
cin >> a >> b >> n;
/* n需要被thread_count整除,否则会漏算*/
if (n % thread_count != 0) {
fprintf(stderr, "n must be evenly divisible by thread_count\n");
exit(0);
}
# pragma omp parallel num_threads(thread_count)
Trap(a, b, n, &global_result);
printf("With n = %d trapezoids, our estimate of the integral\n", n);
printf("from % f to % f = % .14e\n", a, b, global_result);
return 0;
}
double f(double x) {
return -x * x + 9;
}
void Trap(double a, double b, int n, double* global_result_p) {
double h, x, my_result;
double local_a, local_b;
int i, local_n;
int my_rank = omp_get_thread_num();
int thread_count = omp_get_num_threads();
h = (b - a) / n;
local_n = n / thread_count;
local_a = a + my_rank * local_n * h;
local_b = local_a + local_n * h;
my_result = (f(local_a) + f(local_b)) / 2.0;
for (i = 1; i < local_n - 1; i++) {
x = local_a + i * h;
my_result += f(x);
}
my_result = my_result * h;
# pragma omp critical
* global_result_p += my_result;
}
最后累加线程的结果这一步存在竞争条件，因此使用critical指令告诉编译器需要安排线程对相应代码块进行互斥访问，即一次只有一个线程能够执行下面的结构化代码。
我们以
y
=
−
x
2
+
9
y=-x^2+9
y=−x2+9在区间
[
0
,
3
]
[0,3]
[0,3]上的积分为例，线程数为5，分为100000个子区间。得到结果如下图所示：
-
变量的作用域：在OpenMP中，变量的作用域涉及在parallel块中能够访问该变量的线程集合。一个能够被线程组中所有线程访问的变量拥有共享作用域，而一个只能被单个线程访问的变量拥有私有作用域。在parallel指令前已经被声明的变量，拥有在线程组中所有线程间的共享作用域，而在块中声明的变量（例如函数的局部变量）中有私有作用域。另外，在parallel块开始处的共享变量的值，与该变量在parallel块之前的值一样。在parallel块完成后，该变量的值是块结束时的值。
-
归约子句：归约（reduction）就是将相同的规约操作符重复地应用到操作数序列来得到一个结果的计算。所有的操作的中间结果存储在同一个变量里：归约变量(reduction variable)。归约操作符（reduction operator）是一个二元操作（例如+、*）。
reduction 子句的语法是：
reduction(<operator>:<variable list>)
operator可能是操作符+、*、-、&、|、^、&&、|| 中的任意一个，但是减法不满足交换律和结合律，所以不能用在OpenMP归约子句中。同时，如果一个归约变量是float或double型数据，那么当使用不同数量的线程时，由于浮点数运算不满足结合律，即(a+b)+c可能不会准确地等于a+(b+c)，结果可能会有些许不同。
方法2
double Local_trap(double a, double b, int n);
替换
void Trap(double a, double b, int n, double* global_result_p);
除了没有临界区外Local_trap的函数体与Trap的函数体是一样的
#include <iostream>
using namespace std;
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
double Local_trap(double a, double b, int n);
int main(int argc, char* argv[]) {
double global_result = 0.0;
double a, b;
int n;
int thread_count;
thread_count = strtol(argv[1], NULL, 10);
printf("Enter a,b and n\n");
cin >> a >> b >> n;
/* n需要被thread_count整除*/
if (n % thread_count != 0) {
fprintf(stderr, "n must be evenly divisible by thread_count\n");
exit(0);
}
# pragma omp parallel num_threads(thread_count) \
reduction(+:global_result)
global_result += Local_trap(a, b, n);
printf("With n = %d trapezoids, our estimate of the integral\n", n);
printf("from % f to % f = % .14e\n", a, b, global_result);
return 0;
}
double f(double x) {
return -x * x + 9;
}
double Local_trap(double a, double b, int n) {
double h, x, my_result;
double local_a, local_b;
int i, local_n;
int my_rank = omp_get_thread_num();
int thread_count = omp_get_num_threads();
h = (b - a) / n;
local_n = n / thread_count;
local_a = a + my_rank * local_n * h;
local_b = local_a + local_n * h;
my_result = (f(local_a) + f(local_b)) / 2.0;
for (i = 1; i < local_n - 1; i++) {
x = local_a + i * h;
my_result += f(x);
}
my_result = my_result * h;
return my_result;
}
#include<omp.h>
#include <iostream>
using namespace std;
int ReductionTest() {
cout << "reduction输出：\n";
int i = 0, j = 10;
#pragma omp parallel for reduction(+:j) //结果与线程数有关
for (i = 0; i < 100; i++)
{
//j = 2; //如果不对j进行初始化，则最后j的值与线程数无关
cout << i << ":" << j << "\n";
j++;
cout << i << ":" << j << "\n";
}
cout << "最后j的值为:" << j << "\n";
return 0;
}
int main() {
ReductionTest();
}
parallel for 命令
并行化串行梯形积分法
h = (b-a)/n;
approx = (f(a)+f(b))/2.0;
# pragma omp parallel for num_threads(thread_count)\
reduction(+:approx)
for (i = 1; i <= n-1; i++){
approx += f(a+i*h);
}
approx = h*approx;
代码中，approx作为一个归约变量是必要的。如果不这样做，它将是一个普通的共享变量，循环体中的
approx += f(a+i*h);
将会是一个无保护的临界区。
parallel for 指令之后的结构化块必须是for循环。另外，运用parallel指令，系统通过在线程间划分循环迭代来并行化for 循环。在一个被parallel for指令并行化的循环中，循环变量的缺省作用域是私有的，每个线程组中的线程拥有它自己的i的副本。
警告
OpenMP只能并行化那些可以在如下情况下确定迭代次数的for循环：
- 由 for 语句本身（即for(…;…;…)）来确定。
- 在循环执行之前确定。
例如以下两种情况均不能被并行化
//无限循环
for ( ; ; ){
...
}
//break添加了另一个从循环退出的出口
for (i = 0; i < n; i++){
if(...)break;
...
}
运行时系统必须能够在执行前决定迭代的数量，但唯一的例外是：在循环体中可以有一个exit调用。
数据依赖性
考虑下述代码，计算前10个斐波那契（Fibonacci）数：
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
int main(int argc, char* argv[]) {
int n = 10;
int thread_count = 5;
int fibo[10] = { 0 };
fibo[0] = fibo[1] = 1;
# pragma omp parallel for num_threads(thread_count)
for (int i = 2; i < n; i++) {
fibo[i] = fibo[i - 1] + fibo[i - 2];
}
# pragma omp critical
for (int i = 0; i < n; i++) {
printf("%d ", fibo[i]);
}
return 0;
}
运行3次，发现只有第3次结果是正确的。在第一个运行结果中fibo[6]在fibo[4]和fibo[5]的值还未更新的情况下便进行了计算，因此fibo[6]的值为0。
以上例子说明一个或更多个迭代结果依赖于其他迭代的循环，一般不能被OpenMP正确地并行化。fibo[6]和fibo[5]计算间的依赖关系称为数据依赖。由于fibo[5]的值在一个迭代中计算，其结果在之后的迭代中使用，该依赖关系有时称为循环依赖(loop-carried dependence)。
π \pi π值估计
一个对
π
\pi
π的数值估计的方法是使用下列公式：
π
=
4
[
1
−
1
3
+
1
5
−
1
7
+
.
.
.
]
=
4
∑
k
=
0
∞
(
−
1
)
k
2
k
+
1
\pi=4[1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+...]=4\sum_{k=0}^{\infty} \frac{(-1)^k}{2k+1}
π=4[1−31+51−71+...]=4∑k=0∞2k+1(−1)k
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
int main() {
//factor 的类型是double，而不是int或long，为什么这点很重要？
double factor;
double sum = 0.0;
int n = 10000;
int thread_count = 10;
# pragma omp parallel for num_threads(thread_count)\
reduction(+:sum) private(factor)
for (int k = 0; k < n; k++) {
/*
if (k % 2 == 0)
factor = 1.0;
else
factor = -1.0;*/
factor = (k % 2 == 0) ? 1.0 : -1.0;
sum += factor / (2 * k + 1);
}
double pi_approx = 4.0 * sum;
printf("With n = %d and %d threads, Our estimate of pi = %f.", n, thread_count, pi_approx);
return 0;
}
如果第k次迭代被分配给一个线程，而第k+1次迭代被分配给另一个线程，则不能保证factor的值是正确的。因此我们将下述代码：
for (int k = 0; k < n; k++) {
sum += factor / (2 * k + 1);
factor = -factor
}
替换为：
for (int k = 0; k < n; k++) {
factor = (k % 2 == 0) ? 1.0 : -1.0;
sum += factor / (2 * k + 1);
}
从而消除了循环依赖。
缺省情况下，任何在循环前声明的变量（唯一的例外是循环变量）在线程间都是共享的。因此factor是共享的。代码中，private(factor)使得thread_count个线程中的每一个都有它自己的factor变量的副本，因此一个线程对factor的更新不会影响另一个线程的factor值。
最终得到代码的运行结果如下：
- 要记住的重要的一点是，一个有私有作用域的变量的值在parallel块或者parallel for 块的开始处是未指定的。它的值在parallel或paralle for块完成之后也是未指定的。
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
int main() {
int thread_count = 10;
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", my_rank, x);
x = 2 * my_rank + 2;
printf("Thread %d > before initialization, x = % d\n", my_rank, x);
}
printf("After parallel block, x = %d", x);
return 0;
}
当我们运行上述代码时，会报出如下错误
当我们将第一行的 printf(“Thread %d > before initialization, x = % d\n”, my_rank, x);注释后，则会得到如下运行结果。
default(none)子句
OpenMP提供了一个子句default，该子句显式地要求我们明确在块中使用地每个变量和已经在块之外声明的变量和已经在块之外声明的变量的作用域。(在一个块中声明的变量都是私有的，因为它们会被分给线程的栈。)
使用一个default(none)子句对 π \pi π的计算如下所示
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
int main() {
double factor;
double sum = 0.0;
int n = 10000;
int thread_count = 10;
int k;
# pragma omp parallel for num_threads(thread_count)\
default(none) reduction(+:sum) private(k, factor)\
shared(n)
for (k = 0; k < n; k++) {
factor = (k % 2 == 0) ? 1.0 : -1.0;
sum += factor / (2 * k + 1);
}
double pi_approx = 4.0 * sum;
printf("With n = %d and %d threads, Our estimate of pi = %f.", n, thread_count, pi_approx);
return 0;
}
排序
冒泡排序
- 串行冒泡排序核心代码
for (int list_length = n; list_length >= 2; list_length--)
for (int i = 0; i < list_length-1; i++)
if (arr[i] > arr[i + 1]) {
tmp = arr[i];
arr[i] = arr[i + 1];
arr[i + 1] = tmp;
}
在外部循环的任何一次迭代中，当前列表内容依赖于外部循环的前一次迭代。内部循环在第i次迭代中，被比较的元素依赖于第i-1次迭代。如果第i-1次迭代中，a[i-1]和a[i]没有互换，那么第i次迭代将比较a[i]和a[i+1]。如果第i-1次迭代交换了a[i-1]和a[i]，那么第i次迭代将比较原始的a[i-1]和a[i+1]。因此冒泡排序不能用parallel for进行并行化。
奇偶变换排序
- 串行奇偶排序核心代码
for (int phase = 0; phase < n; phase++){
if (phase % 2 == 0) {
for (int i = 1; i < n; i += 2) {
if (arr[i - 1] > arr[i]) swap(&arr[i - 1], &arr[i]);
}
}
else{
for (int i = 1; i < n - 1; i += 2) {
if (arr[i] > arr[i+1]) swap(&arr[i], &arr[i+1]);
}
}
}
列表arr存储n个整数，算法对它们进行升序排序。在一个“偶阶段”（phase % 2 == 0）里，每个偶下标元素arr[i]与它左边的元素arr[i-1]相比较。如果没有排好序就交换它们。在一个”奇阶段“里，每个奇下标元素与它右边的元素相比较。
- 奇偶排序的第一个OpenMP实现
#include <iostream>
using namespace std;
#include <omp.h>
void swap(int* a, int* b)
{
int temp;
temp = *a;
*a = *b;
*b = temp;
}
int main() {
int N = 10, m = 0, arr[10] = {10,9,8,7,6,5,4,3,2,1}, i, tmp;
int thread_count = 2;
/*
srand(1234); //设置随机数种子
for (int j = N;j>0;j--)
{
arr[m++] = rand();
}
*/
for (int phase = 0; phase < N; phase++){
if (phase % 2 == 0) {
# pragma omp parallel for num_threads(thread_count)\
default(none) shared(arr,N) private(i,tmp)
for (i = 1; i < N; i += 2) {
if (arr[i - 1] > arr[i]) swap(&arr[i - 1], &arr[i]);
}
}
else{
# pragma omp parallel for num_threads(thread_count)\
default(none)shared(arr, N) private(i, tmp)
for (i = 1; i < N - 1; i += 2) {
if (arr[i] > arr[i+1]) swap(&arr[i], &arr[i+1]);
}
}
}
printf("Ordered array:");
for (int i = 0; i < N; i++)
printf("%d ", arr[i]);
return 0;
}
考虑创建与合并线程的开销，上面的代码会在每一遍外部循环都创建和合并thread_count个线程。修改为只创建一次线程，并在每次内部循环的执行中重复使用它们，代码如下
- 奇偶排序的第二个OpenMP实现
#include <iostream>
using namespace std;
#include <omp.h>
void swap(int* a, int* b)
{
int temp;
temp = *a;
*a = *b;
*b = temp;
}
int main() {
int N = 10, m = 0, arr[10] = {10,9,8,7,6,5,4,3,2,1}, i, tmp;
int thread_count = 2;
int phase;
/*
srand(1234); //设置随机数种子
for (int j = N;j>0;j--)
{
arr[m++] = rand();
}
*/
# pragma omp parallel num_threads(thread_count)\
default(none) shared(arr,N) private(i,tmp,phase)
for (phase = 0; phase < N; phase++){
if (phase % 2 == 0) {
# pragma omp for
for (i = 1; i < N; i += 2) {
if (arr[i - 1] > arr[i]) swap(&arr[i - 1], &arr[i]);
}
}
else{
# pragma omp for
for (i = 1; i < N - 1; i += 2) {
if (arr[i] > arr[i+1]) swap(&arr[i], &arr[i+1]);
}
}
}
printf("Ordered array:");
for (int i = 0; i < N; i++)
printf("%d ", arr[i]);
return 0;
}
循环调度
生产者和消费者问题
缓存、缓存一致性、伪共享
线程安全性
利用OpenMP实现以下功能和算法：
- 图像直方图计算
- 实现直方图均衡化和直方图规定化算法
- 计算图像梯度的模，并获得梯度模的总和
扫一扫
951
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
