CS61c_18 — Thread-Level Parallelism II

引

在上一篇中，我们介绍了thread level parallelism的基本概念和多核CPU架构以及超线程，现在我们将落实到具体的编程实践中，看看如何使用并行化来加速程序的运行。
当前，众多编程语言支持并行化编程，比如 Go、C、Java、CUDA 等，它们针对不同领域进行了特别优化。本文将以 C 语言为例，介绍一些常用的并行化编程工具和库，比如 OpenMP 和 Pthreads，帮助我们更方便地实现并行计算。

值得注意的是，为什么在 C 中我们需要使用 intrinsics 来嵌入 SIMD 指令呢？ 这是因为 GCC 等编译器目前仍难以在通用场景下高效地将普通 C 代码自动并行化，所以我们需要手动告诉编译器哪些部分可以并行、哪些指令可以用 SIMD 加速。
不过，Inter等公司正在努力升级编译器，将会有更多的并行化优化。20+ years to translate C into good (fast!) assembly, How long to translate C into good (fast!) parallel assembly ? 呃呃，想要做到通用并行化确实很困难，所以各大语言基本上都是针对特定场景进行优化(比如，Scienctific computing/Machine learning, Webserver, Input/Output). 如果你有志为并行化编译器做贡献，大有前途啊小伙子！

OpenMP

Parallel Loops

下面是一个使用 OpenMP 来并行化 for 循环的示例代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
#include <stdio.h>
#include <omp.h>
int main() {
	omp_set_num_threads(4);
	int a[] = {0,1,2,3,4,5,6,7,8,9};
	int N = sizeof(a)/sizeof(int);

	#pragma omp parallel for
	for (int i=0; i<N; i++){
		printf("thread %d,i = %2d\n", omp_get_thread_num(), i);
		a[i] = a[i] + 10 * omp_get_thread_num();
	}

	for (int i=0; i<N; i++){
		printf("%2d ", a[i]);
	}
	printf("\n");
}

// use gcc -fopenmp for.c -o for.out to compile

你会发现线程0处理了i=0,1,2 , 线程1处理了i=3,4,5 , 线程2处理了i=6,7 , 线程3处理了i=8,9. 这就是OpenMP的默认调度策略，将一块大的共享内存分成若干个不同的区间，每个线程负责处理其中一个区间。并且OpenMP支持动态分配，比如这里10并不是4的倍数，它会自动分配成3，3，2，2.
需要注意，如果你多次运行./for.out, 你会发现每次各个线程输出顺序是不一样的，甚至有错位的情况，这就是并行化的非确定性。如果你的程序依赖于某些特定的输出顺序，请务必使用适当的同步机制来保证正确性！
请慎用time()或clock()来测量并行程序的时间！！ OpenMP提供了专门的函数omp_get_wtime()来测量并行程序的运行时间！

Computing π

下面我们将尝试使用OpenMP来并行计算π的近似值。问题来了，如何计算π呢？可以想一想，π在哪里出现过？比如圆的面积公式，周长公式，弧度制等。于是，有以下两种经典的方法：

对于一个单位圆，面积是π, 对于第一象限的四分之一部分，面积是$\frac{π}{4}$, $y = \sqrt{1-x^2}$, 因此我们可以用积分的方法来逼近这个面积，从而得到π的近似值。即 $$π = 4 \int_0^1 \sqrt{1-x^2} dx$$
另一个方法是弧度制，已知$\tan(\frac{π}{4}) = 1$, $\tan(0) = 0$, $\arctan(x)$的导数为$\frac{1}{1+x^2}$, 于是，利用函数积分： $$ \int_0^1 \frac{1}{1+x^2} dx = \arctan(1) - \arctan(0) = \frac{π}{4} $$ 故： $$ π = 4 \int_0^1 \frac{1}{1+x^2} dx $$

以上两种方法都可以通过数值积分来近似计算π的值，不过考虑到$\sqrt{1-x^2}$的计算比较慢，所以下面我们的代码采用方法二。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
// 先写个单线程版本

#include <stdio.h>

void main() {
	const long num_steps = 10;
	double step = 1.0/((double)(num_steps));
	double sum = 0.0;
	for (int i=0; i<num_steps; i++){
		double x = (i+0.5) * step;
		sum += 4.0 * step/(1.0 + x*x);
	}
	printf("pi = %6.12f\n",sum);
}

// 输出结果：pi = 3.142425985001

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// 考虑并行多线程
#include <stdio.h>
#include <omp.h>

void main() {
	const long num_steps = 10;
	omp_set_num_threads(4);
	double step = 1.0/((double)num_steps);
	double sum = 0.0;

	#pragma omp parallel for 
	for (int i=0; i<num_steps; i++){
		double x = (i+0.5) * step;
		sum += 4.0*step/(1.0+x*x);
	}

	printf("pi = %6.12f\n", sum);
}

/*
这是最初想法，运用了之前提到的#pragma parallel for来并行化for循环, 
然而，这段代码存在一个严重问题：多个线程同时访问和修改共享变量sum，可能导致结果错误！
实测结果：pi = 2.146336486161 , 确实出现严重错误！
*/

解决办法：由于问题产生原因是多个线程对同一个共享变量sum进行修改，那么我们就为每一个线程分配一个独立的sum变量，最后再将它们加起来。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/* 
由于我们是对不同thread分配一个sum变量，所以每步计算时需要知道当前thread id, 
然后修改对应的sum[id]!
*/

/*
这是教学视频中提供的写法，使用了#pragma omp parallel来创建一个并行区域，先确定哪个thread，
然后以THREAD_NUM为步长来分配循环迭代。最后再将每个线程的sum加起来得到最终结果。
*/

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

void main() {
	const int NUM_THREADS = 4;
	const long num_steps = 10;
	double step = 1 / ((double)num_steps);
	double sum[NUM_THREADS];
	for (int i=0; i<NUM_THREADS; i++) sum[i] = 0.0;

	omp_set_num_threads(NUM_THREADS);

#pragma omp parallel
	{
		int id = omp_get_thread_num();
		for (int i=id ; i<num_steps; i+= NUM_THREADS){
			double x = (i+0.5) * step;
			sum[id] += 4.0 * step/(1 + x*x);
			printf("i =%3d, id =%3d\n", i, id);
		}
	}
	double pi = 0;
	for(int i=0; i<NUM_THREADS; i++) pi += sum[i];
	printf("pi = %6.12f\n", pi);
}

/*
当然，用之前提到过的#pragma omp parallel for来写也完全可以，在for循环内部计算前确定当前线程id即可！
*/

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

void main() {
	const int NUM_THREADS = 4;
	const long num_steps = 10;
	double step = 1 / ((double)num_steps);
	double sum[NUM_THREADS];
	for (int i=0; i<NUM_THREADS; i++) sum[i] = 0.0;

	omp_set_num_threads(NUM_THREADS);

	#pragma omp parallel for
	for (int i=0 ; i<num_steps; i+=1 ){
		int id = omp_get_thread_num();
		double x = (i+0.5) * step;
		sum[id] += 4.0 * step/(1 + x*x);
		printf("i =%3d, id =%3d\n", i, id);
	}

	double pi = 0;
	for(int i=0; i<NUM_THREADS; i++) pi += sum[i];
	printf("pi = %6.12f\n", pi);
}


/*两种结果均正确，成功实现并行计算π. 输出结果如下：

方法一：
i =  0, id =  0
i =  4, id =  0
i =  8, id =  0
i =  3, id =  3
i =  7, id =  3
i =  2, id =  2
i =  6, id =  2
i =  1, id =  1
i =  5, id =  1
i =  9, id =  1
pi = 3.142425985001

方法二：
i =  0, id =  0
i =  1, id =  0
i =  2, id =  0
i =  8, id =  3
i =  9, id =  3
i =  6, id =  2
i =  7, id =  2
i =  3, id =  1
i =  4, id =  1
i =  5, id =  1
pi = 3.142425985001

如果从内存隔离角度考虑，方法二会好不少。
*/

对于可以开的最大Thread数量，你可以使用lscpu去查看你的CPU有多少个核心(上一篇中已经详细介绍)，或者使用OpenMP提供的函数omp_get_max_threads()来获取当前系统支持的最大线程数。比如，我偷偷用社团的开发机同时开128个线程核心，在num_steps=1000000000的情况下，能从1.5s优化到0.5s！
当然，num_steps越大，最终π的估计越精准。
可能有人疑惑了，如果能开无限多个线程，是不是运行时间趋近于0 ？对于这个问题，答案是否定的，因为最后我们还得将sum[i]进行求和，而这个过程是竞争态的，必须串行！如果没有锁的保护，那么两个线程同时读当前sum值，并且一前一后写回，那么先写回的sum[i]就被忽略了！！！因此，在并行编程中，竞争态是一个非常常见的问题，必须通过适当的同步机制来解决！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// 演示一下竞争态错误

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

void main() {
	const int NUM_THREADS = 100;
	const long num_steps = 10000;
	double step = 1 / ((double)num_steps);
	double sum[NUM_THREADS];
	for (int i=0; i<NUM_THREADS; i++) sum[i] = 0.0;
	double pi = 0.0;

	omp_set_num_threads(NUM_THREADS);

	#pragma omp parallel
	{
		int id = omp_get_thread_num();
		for (int i=id ; i<num_steps; i+= NUM_THREADS){
			double x = (i+0.5) * step;
			sum[id] += 4.0 * step/(1 + x*x);
			// printf("i =%3d, id =%3d\n", i, id);
		}
		pi += sum[id];
	}
	printf("pi = %6.12f\n", pi);
}

// pi = 0.376985602277 ？！！

前面提到最后的求和过程必须串行，如果我们硬要在并行块内进行求和，那么会得到奇怪的π值。一个程序最后的优化瓶颈其实是它的串行部分，无论你开多少线程，串行部分的时间都不会减少，所以并行化的效果是有限的，这就是所谓的Amdahl’s Law。
需要注意的是，在设置THREAD_NUM=100时，num_steps设为10000或更少一点才能看到明显的错误，因为如果num_steps太大，那么不同线程迭代次数很多，完成迭代的时间差可能会比较大，导致竞争态的概率降低，反而得到正确结果🤣

Synchronization

在前面的代码中，我们在尝试对sum进行并行计算时出现了竞争态问题。在生活中，一些公共资源同样会出现类似的问题，比如一份共享文档，如果没有适当的权限控制，那么多个用户同时编辑同一处内容时，会出现冲突，导致文档内容混乱。为了解决这个问题，一个很自然的想法是take turns! 当我在编辑文档时，其他人？等着！当我编辑完了，其他人才开始编辑。这种思想在计算机中被称为"Lock"。
计算机使用Locks来控制对共享资源的访问，通常使用一个 int lock 变量来实现，即 0 for unlocked, 1 for locked.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// wait for lock released
while (lock != 0) ;
// lock == 0 now (unlocked)

// set lock
lock = 1;

// access shared resource ...
// e.g. π
// sequential execution! (Amdahl ...)

// release lock
lock = 0;

上面是一个lock的示例，很简单也看似正确。然而，如果有两个线程时会发生什么呢？假设线程A和线程B同时执行到while (lock != 0)这行代码，此时lock=0，两个线程都会跳出循环，接着令lock=1，两个线程都认为自己获得了锁，接着访问共享资源，导致竞争态问题！因此，这个简单的lock实现是错误的。

而且，这个问题在C语言层面是无法解决的。我们需要保证一个线程在发现lock=0并且将lock设置为1的过程中，其他线程无法访问lock变量，这这这，用C的语法肯定做不到呀！我们需要一些更底层的“原子操作”来解决这个问题！
下一篇我们将深入探讨如何解决这个问题！