OpenMP 归约和reduction子句

简述归约
　　归约操作在MPI里也学过，不过那时候还不太熟悉这种操作。当时只知道MPI_Reduce可以把全局求和和集合通信封装起来，非常方便。实际上将相同的二元归约操作符重复地应用到一个序列上得到结果的计算过程都可以称为归约。

　　python里那个难理解的reduce()函数也就是归约：

1 >>> from functools import reduce
2 >>> def myfun(x,y):
3 ...     return x+y-1
4 ... 
5 >>> reduce(myfun,[1,2,3])
6 4
7 >>> 

reduction子句
　　OpenMP中的归约是parallel并行指令的reduction子句，在子句中指定归约操作符和归约变量。

　　归约操作符是序列中的两两元素做的运算，一定是一个二元运算符。归约变量则保存归约操作的中间结果。OpenMP用归约变量为每个线程创建一个私有的变量，用来存储自己归约的结果，所以归约的代码不需要critical保护也不会发生冲突：

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<omp.h>
 4 
 5 int main(int argc,char *argv[])
 6 {
 7     int sum=20;
 8     int thrdCnt=strtol(argv[1],NULL,10);
 9     //归约子句(归约操作符:归约变量)
10 #   pragma omp parallel num_threads(thrdCnt) reduction(+:sum)
11     {
12         int myRank=omp_get_thread_num();
13         sum+=myRank;
14         printf("%d->%d\n",myRank,sum);
15     }
16     printf("sum=%d\n",sum);//归约结果
17     return 0;
18 }

输出

 1 [lzh@hostlzh OpenMP]$ !gcc
 2 gcc -fopenmp -o test1.o test1.c
 3 [lzh@hostlzh OpenMP]$ ./test1.o 10
 4 9->9
 5 1->1
 6 2->2
 7 3->3
 8 4->4
 9 7->7
10 8->8
11 5->5
12 0->0
13 6->6
14 sum=65
15 [lzh@hostlzh OpenMP]$

可以看到在加法归约时，为每个线程创建的私有变量初始值是0，即使是0号线程也不例外。最后再把各自归约后的私有变量归约到归约变量上。

比较特殊的是，在减法时私有变量的初始值也是0，最后再把这些值按加法运算归约到归约变量上。这是因为减法不满足结合率，A-B-C-D-E-F和(A-B)-(C-D)-(E-F)是截然不同的。

除法运算也不满足结合率，但OpenMP的归约运算不包括除法运算。

虽然浮点数可以作归约，但如果对精度有要求，可能需要注意浮点数不满足结合率。

当进行乘法归约时，私有变量的初始值就是1了：

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<omp.h>
 4 
 5 int main(int argc,char *argv[])
 6 {
 7     int sum=-2;
 8     int thrdCnt=strtol(argv[1],NULL,10);
 9     //归约子句(归约操作符:归约变量)
10 #   pragma omp parallel num_threads(thrdCnt) reduction(*:sum)
11     {
12         int myRank=omp_get_thread_num();
13         if(myRank!=0)
14             sum*=myRank;
15         printf("%d->%d\n",myRank,sum);
16     }
17     printf("sum=%d\n",sum);//归约结果
18     return 0;
19 }

输出

1 [lzh@hostlzh OpenMP]$ !gcc
2 gcc -fopenmp -o test1.o test1.c
3 [lzh@hostlzh OpenMP]$ ./test1.o 4
4 3->3
5 1->1
6 2->2
7 0->1
8 sum=-12
9 [lzh@hostlzh OpenMP]$

归约使用限制

允许使用的归约操作符

　　只有+,*,-,&,|,^,&&,||八种。

对归约变量的操作

　　在reduction子句修饰的并行块中，在书写逻辑上只允许对归约变量做如下的写入操作(之所以说是书写逻辑上，因为实际操作的始终是私有变量)：

1 x = x op expr
2 x = expr op x (除了减法)   
3 x binop = expr
4 x++
5 ++x
6 x--
7 --x

　　在并行块中对归约变量(实际上是对私有变量)做这些操作往往使用的是归约本身的操作符，如果使用其它操作符，只要满足上面的条件就是允许的，不过最后各个线程得到的结果仍然会按reduction子句指定的归约操作符所应当的归约方式(如+和-即做加法，*即做乘法)进行归约，而不会真的按照在并行块中表面上对归约变量的操作去执行：

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<omp.h>
 4 
 5 int main(int argc,char *argv[])
 6 {
 7     int sum=-2;
 8     int thrdCnt=strtol(argv[1],NULL,10);
 9     //归约子句(归约操作符:归约变量)
10 #   pragma omp parallel num_threads(thrdCnt) reduction(*:sum)
11     {
12         int myRank=omp_get_thread_num();
13         if(myRank!=0)
14             sum*=myRank;
15         sum++;//实际是加在私有变量上!
16         printf("%d->%d\n",myRank,sum);
17     }
18     printf("sum=%d\n",sum);//归约结果
19     return 0;
20 }

输出

1 [lzh@hostlzh OpenMP]$ !gcc
2 gcc -fopenmp -o test1.o test1.c
3 [lzh@hostlzh OpenMP]$ ./test1.o 4
4 3->4
5 1->2
6 2->3
7 0->2
8 sum=-96
9 [lzh@hostlzh OpenMP]$

　　这个-96来自(1*3+1)*(1*1+1)*(1*2+1)*(1+1)*(-2)，总之理解好OpenMP归约的本质是在并行块里按归约操作符来生成指定初始值的私有变量，在并行块中看似对归约变量的操作是对私有变量的操作，最终再将各个线程计算好的私有变量按归约操作符的形式以特定的方式归约。

　　特别注意最后的归约时’-‘操作符做加操作！下面这个例子能再次证明这一点：

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<omp.h>
 4 
 5 int main(int argc,char *argv[])
 6 {
 7     int sum=30;
 8     int thrdCnt=strtol(argv[1],NULL,10);
 9     //归约子句(归约操作符:归约变量)
10 #   pragma omp parallel num_threads(thrdCnt) reduction(-:sum)
11     {
12         int myRank=omp_get_thread_num();
13         sum+=2;
14         printf("%d->%d\n",myRank,sum);
15     }
16     printf("sum=%d\n",sum);//归约结果
17     return 0;
18 }

输出

1 [lzh@hostlzh OpenMP]$ !gcc
2 gcc -fopenmp -o test1.o test1.c
3 [lzh@hostlzh OpenMP]$ ./test1.o 4
4 3->2
5 1->2
6 2->2
7 0->2
8 sum=38
9 [lzh@hostlzh OpenMP]$