积分函数
//积分函数
double blog4::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;
}
数学函数
//数学函数
double blog4::f(double x)
{
return pow(x, 3);
}
OpenMP代码
//第二个案例
void blog4::Test2(int argc, char* argv[])
{
int thread_count = strtol(argv[1], NULL, 10); double global_result = 0.0; #pragma omp parallel num_threads(thread_count)
{
double my_result = 0.0;
my_result = Trap(0, 3, 1024);
#pragma omp critical
global_result += my_result;
} printf("%f\n", global_result);
}
归约
void blog4::Test2(int argc, char* argv[])
{
int thread_count = strtol(argv[1], NULL, 10); double global_result = 0.0; #pragma omp parallel num_threads(thread_count) reduction(+: global_result)
{
global_result += Trap(0, 3, 1024);
} printf("%f\n", global_result);
}
运行：
后面的2指定使用2个线程运行程序。
通过以上的例子，大概知道OpenMP只要添加一条很简单的parallel for指令，就能够并行化大量的for循环所组成的串行程序。但使用它有几天限制条件：
（1）、OpenMP只能并行化for循环，它不会并行while和do-while循环，而且只能并行循环次数在for循环外面就确定了的for循环。
（2）、循环变量只能是整型和指针类型（不能是浮点型）
（3）、循环语句只能是单入口单出口的。循环内部不能改变index，而且里面不能有goto、break、return。但是可以使用continue，因为它并不会减少循环次数。另外exit语句也是可以用的，因为它的能力太大，他一来，程序就结束了。
3.数据依赖问题
使用OpenMP在明面上有上面的这些限制规则，而需要关注的是一个更隐匿的问题：数据依赖问题。这其中涉及两种依赖情况：数据依赖和循环依赖。为了弄清楚这两种依赖的具体是什么，我们先来看一个例子。
计算前n个斐波那契（fibonacci）数：
然后加上OpenMP并行
通过测试，我们有时会得到：
1 1 2 3 5 8 13 21 34 55（正确）
然后也可能偶尔会得到：
1 1 2 3 5 8 0 0 0 0.（错误）
这究竟是发生了什么？似乎运行时，系统将fibo[2], fibo[3], fibo[4], 和 fibo[5]的计算分配给了一个线程，将fibo[6], fibo[7],fibo[8], 和 fibo[9]分配给了另外一个线程。所以在计算fibo[6]的时候，如果第一个线程在此线程之前已经完成，入口初始化为fibo[5]=8，得到的计算结果就是正确的；而如果在计算fibo[6]的时候，如果第一个线程在此线程之前还没有完成，入口初始化为fibo[5]=0，则得到的计算结果就是错误的。
从这里，我们看到两个要点：
（1）OpenMP编译器不检查被parallel for指令并行化的循环所包含的迭代间的依赖关系，而是由程序员来识别这些依赖关系。
（2）一个或更多个迭代结果依赖于其它迭代的循环，一般不能被OpenMP正确的并行化。
从求斐波那契（fibonacci）数的这个例子中，我们发现fibo[3], fibo[4]计算间的依赖关系称之为数据依赖。而fibo[5]的值在一个迭代中计算，它的结果是在之后的迭代中使用，这样的依赖关系称之为循环依赖。
4.怎么寻找循环依赖及解决办法
当我们试图使用一个parallel for指令时，需要小心循环依赖关系，而不用担心数据依赖。例如，在下列循环中：
int i = 0;
for (i = 0; i < n; i++)
{
x[i] = a + i*h;
y[i] = exp(x[i]);
}
在这个程序中y[i]数据依赖于x[i]的，这种依赖使用OpenMP并行化的话，是没有关系的。
在看一个例子，如下公式计算π：
使用能够在串行代码中实现的公式
//第5个案例
void blog4::Test5(int argc, char* argv[])
{
double factor = 1;
double sum = 0.0;
int n = 1000;
for (int i = 0; i < n; i++)
{
sum += factor / (2 * i + 1);
factor = -factor; }
double pi_approx = 4.0*sum;
printf("%f", pi_approx);
}
得到的结果3.140593，符合预期。我们在上面这个段代码中加上OpenMP的并行化指令
//第5个案例
void blog4::Test5(int argc, char* argv[])
{
double factor = 1;
double sum = 0.0;
int n = 1000;
int thread_count = strtol(argv[1], NULL, 10);//omp_get_num_threads();
cout << thread_count << endl;
#pragma omp parallel for num_threads(thread_count) reduction(+:sum)
for (int i = 0; i < n; i++)
{
sum += factor / (2 * i + 1);
factor = -factor;
//cout << i << "-" << factor << endl;
}
double pi_approx = 4.0*sum;
printf("%f", pi_approx);
}
运行程序
得到结果2.730014，不符合预期。这是存在一个循环依赖，不能保证第k次迭代中对factor的更新对第k+1次迭代有影响，所以导致结果不正确。通过分析，对代码进行相应的改动。
if (i % 2 == 0)
factor = 1.0;
else
factor = -1.0;
sum += factor / (2 * i + 1);
运行结果为：2.976587，依旧不正确，为什么呢？
缺省情况下，任何在循环前声明的变量，在线程间都是共享的，。为了消除这种循环依赖关系之外，还需要保证每个线程有它自己的factor副本，OpenMP的语句为我们考虑了，通过添加private子句到parallel指令中来实现这一目标。
//第5个案例
void blog4::Test5(int argc, char* argv[])
{
double factor = 1;
double sum = 0.0;
int n = 1000;
int thread_count = strtol(argv[1], NULL, 10);//omp_get_num_threads();
cout << thread_count << endl;
#pragma omp parallel for num_threads(thread_count) reduction(+:sum) private(factor)
for (int i = 0; i < n; i++)
{
if (i % 2 == 0)
factor = 1.0;
else
factor = -1.0;
sum += factor / (2 * i + 1);
//factor = -factor;
//cout << i << "-" << factor << endl;
//SumForNumber();
}
double pi_approx = 4.0*sum;
printf("%f", pi_approx);
}
这下终于正确了。