OpenMP优化for循环的基础运用

OpenMP作为多线程并行优化API，其使用方式与C++自带的多线程使用方式有很大的不同。

在使用OpenMP时，我们是通过 #pragma omp+字句 所组成的命令对线程的行为进行控制，之后编译器会自动对这些命令进行分析与优化，将相关代码由串行变为并行。

整个过程中编译器已经替我们做了相当多的工作，大多数情况下只需要略微的改动就能将程序由串行转化为并行，从而达到成倍的性能提升。

预先准备

因为是让编译器进行优化，所以需要加上-fopenmp的编译选项来开启这一功能。

在未启用openmp的情况下编译器会直接忽视#pragma omp相关命令，将其当成普通程序编译，并不会发出提醒或报错。

如果要判断openmp是否成功启用的话，可以对宏 _OPENMP 定义进行检测

#ifdef _OPENMP
    puts("OpenMp available");
#else
    puts("OpenMp unavailable");
#endif

1.parallel简单并行

这里以Hello World为例：

#include <stdio.h>

int main()
{
    printf("Hello World!\n");
}

很明显，执行后的输出是

Hello World!

然后加上OpenMP命令

#include <stdio.h>

int main()
{
    #pragma omp parallel
    printf("Hello World!\n");
}

输出结果变成了

Hello World!
Hello World!
Hello World!
Hello World!
Hello World!
Hello World!
Hello World!
Hello World!

可以看到Hello World被输出了8次。

这是因为#pragma omp parallel之后的代码被转化为了多线程执行，而输出8次就代表他使用了8线程。

2.num_threads设置线程数

默认的线程数量和CPU有关，不过我们也可以自行指定。

一种方法是设置 OMP_NUM_THREADS 环境变量，

另一种则是使用 num_threads(n) 进行调整，括号里的n就代表线程数。

#include <stdio.h>

int main()
{
    int threads_num = 2; // 甚至可以用变量

    #pragma omp parallel num_threads(threads_num)
    printf("Hello World!\n");
}

现在就是调用2个线程运行了。

Hello World!
Hello World!

3.for循环

大多数情况下的并行需求都出现在for循环中，而OpenMP自然也有对for循环进行优化。

这里我们导入 <omp.h> 头文件，里面包含了OpenMP提供的函数方便调试和控制。

其中的 omp_get_thread_num() 可以返回执行这段代码的线程id。

用一个简单循环打印程序进行演示：

#include <stdio.h>
#include <omp.h>

int main()
{
    #pragma omp parallel for
    for (int i = 0; i < 10; i++)
        printf("i=%d thread:%d\n", i, omp_get_thread_num());
}

输出如下：

i=0 thread:0
i=1 thread:0
i=7 thread:5
i=5 thread:3
i=6 thread:4
i=8 thread:6
i=4 thread:2
i=9 thread:7
i=2 thread:1
i=3 thread:1

可以看出循环变为了乱序执行，这也反映出了循环并行的效果。

但同时也表明，进行多线程处理的for循环一定是独立并且前后不相干，否则这种乱序执行的特性就会导致BUG出现。

for (int i = 1; i < 10; i++)
{
    a[i] = a[i - 1] + 1;
}

像是上面这种进行数据修改并且前后依赖的循环就会因为乱序执行而出BUG。

4.三种private

既然有多个线程在同时运行，那么他们对变量是如何进行管理的呢？

默认条件下有两种情况：

创建线程之前就存在的变量——所有线程共享
线程运行时创建的局部变量——每个线程独有

既然这只是默认情况，那就代表我们其实可以用命令进行修改。

先写一个普通的二层嵌套循环，这里的i和j就是在多线程循环开始之后才进行创建的，属于情况2。

因此每个线程都有属于自己的i和j的拷贝，修改和访问时不会相互干扰，从而保证了循环的正常进行。

#include <stdio.h>

int main()
{
    #pragma omp parallel for
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            printf("i=%d j=%d\n", i, j);
}

运行一下

i=2 j=0
i=2 j=1
i=2 j=2
i=0 j=0
i=0 j=1
i=0 j=2
i=1 j=0
i=1 j=1
i=1 j=2

输出的结果也符合我们的预期。

现在我们把i、j的创建放到开启多线程之前。

就变成了情况1，此时i和j就变成了 所有线程共享 的变量。

#include <stdio.h>

int main()
{
    int i, j;

    #pragma omp parallel for
    for (i = 0; i < 3; i++)
        for (j = 0; j < 3; j++)
            printf("i=%d j=%d\n", i, j);
}

这时如果我们再次运行

i=1 j=0
i=1 j=1
i=1 j=2
i=0 j=0
i=2 j=0

这时~~期待已久的~~BUG就出现了。

这就是多个线程同时访问和修改相同变量导致的问题，

private

为此OpenMP也非常贴心地提供了 private() 命令来解决。

只需要加上 private() 并在括号里写上变量名，

就能让 每个线程都有该变量的独立拷贝，相当于一个同名的局部变量。

#include <stdio.h>

int main()
{
    int i, j;

    #pragma omp parallel for private(i, j)
    for (i = 0; i < 3; i++)
        for (j = 0; j < 3; j++)
            printf("i=%d j=%d\n", i, j);
}

虽然private声明变量后，每个线程都会生成一个相应的拷贝。

但这些线程并不会对他们进行初始化。

#include <stdio.h>

int main()
{
    int x = -1;

    #pragma omp parallel for private(x)
    for (int i = 0; i < 5; i++)
        printf("x=%d\n", x);
}

可以看到拷贝出的变量是什么值都有，但就是没有和原来的值相同的。

x=14908664
x=13903992
x=13909720
x=13905304
x=8

firstprivate

想要有初始化的话就需要用到 firstprivate()。

#include <stdio.h>

int main()
{
    int x = -1;

    #pragma omp parallel for firstprivate(x)
    for (int i = 0; i < 5; i++)
        printf("x=%d\n", x);
}

这样每个变量就都能初始化了

x=-1
x=-1
x=-1
x=-1
x=-1

虽然现在我们有了初始化，但循环结束后变量的值还是无法保留。

继续举个例子测试一下：

（private和firstprivate同理）

#include <stdio.h>

int main()
{
    int x = -1;
    printf("start x=%d\n", x);

    #pragma omp parallel for private(x)
    for (int i = 0; i < 5; i++)
        printf("x=%d\n", x = i);

    printf("final x=%d", x);
}

输出结果

start x=-1
x=1
x=2
x=3
x=0
x=4
final x=-1

可以看到使用private时，循环结束后x的值是不会保留的。

lastprivate

这时候 lastprivate() 就派上用场了，它的功能就是在private的基础上，能够在循环结束时保留变量的值。

改成lastprivate。

#include <stdio.h>

int main()
{
    int x = -1;
    printf("start x=%d\n", x);

    #pragma omp parallel for lastprivate(x)
    for (int i = 0; i < 5; i++)
        printf("x=%d\n", x = i);

    printf("final x=%d", x);
}

输出结果

start x=-1
x=3
x=2
x=1
x=4
x=0
final x=4

变量是成功保存了，但为什么明明最后一个输出的是3，保存的结果却是4呢？

这是因为lastprivate保存的变量是 逻辑上的最后的值。

从代码运行逻辑上来讲x最后的值是4，所以结果就是4。

说通俗点就是和单线程运行的结果相同，不受乱序执行的影响。

组合使用

从功能上可以看出，firstprivate是对功能private的扩展，二者是相互替代的关系，所以 不能同时使用（会编译失败）。

lastprivate和private同样也是相互替代的关系，依旧 不能同时使用（会编译失败）。

但是firstprivate和lastprivate在功能上有相互补充关系，所以 可以同时使用。

5.reduction

假如现在我们要做一个高斯最喜欢的数学题，进行1~100的求和。

使用循环累加并且加上OpenMP并行优化。

#include <stdio.h>

int main()
{
    int sum = 0;

    #pragma omp parallel for
    for (int i = 1; i <= 100; i++)
        sum += i;

    printf("sum=%d", sum);
}

稍微尝试几次之后就会发现，有时候输出的答案并不是5050。

到现在你应该可以很容易就能想到，这是多个线程同时访问sum时产生冲突导致的。

但在这种情况下用private就没办法对其进行求和，此时 reduction 就派上用场了。

reduction的命令格式是 reduction(operation : variable)，其中operation是操作类型，variable则是操作变量。

reduction的作用就是给每个线程创建一个独立的变量，在结束后根据操作类型进行归约。

默认操作包括

算数运算：+, *, -, max, min
逻辑运算：&&, ||
位运算：&, |, ^

我们需要对sum进行累加，所以应该使用 reduction(+ : sum)

#include <stdio.h>

int main()
{
    int sum = 0;

    #pragma omp parallel for reduction(+ : sum)
    for (int i = 1; i <= 100; i++)
        sum += i;

    printf("sum=%d", sum);
}

这样就能正常得到结果了。

6.对于首个for循环迭代器的限制

还有一个比较细节的地方就是OpenMP在优化for语句时，会自动把第一个循环的迭代器（也就是i）设为private。

比如在前面演示private时举例的双重for循环中，如果我们将

#pragma omp parallel for private(i, j)

修改为

#pragma omp parallel for private(j) // 不显式声明i为private

修改后并不会影响循环的正常运行。

而从下面这个例子则能够更好地证实这一点。

在代码中，我们将循环外ij的初值设为-1，并在循环内不断改变i和j的值，观察循环结束后二者的变化情况。

#include <stdio.h>

int main()
{
    int i = -1, j = -1;

    #pragma omp parallel for
    for (i = 0; i < 5; i++)
    {
        printf("i=j=%d\n", j = i);
    }

    printf("final i=%d j=%d\n", i, j);
}

运行结果为：

i=j=0
i=j=2
i=j=1
i=j=4
i=j=3
final i=-1 j=3

可以看出循环体内外的两个i的值是不同的，循环体里的i更像是for循环的局部变量。

这是OpenMP为了保证循环能够正常运作而进行的优化，

但这种默认且强制性的优化的代价是语法规则的拘束，一些正常可以通过编译的代码，放在parallel for的首个循环就会出现编译错误的情况。

#include <stdio.h>

int main()
{
    int i, j;

    #pragma omp parallel for // 去掉for优化后才能编译
    for (i = 0, j = 0; i < 5; i++) // 只是加了个j=0就报错了
    {
    }
}