利用 OpenMP CUDA MPI 并行计算 PI 的值，并与串行版本做简要的对比
# 准备
# 计算环境
- CUDA 11.0、clang 10/gcc 10
- manjaro 20
- CPU: i7-6700hp
- 显卡：gtx965m
- 内存：2133 16G
# 算法选择
- 蒙特卡洛方法 - 具体可参见 Wiki (opens new window)
大致的思路是，在 1/4 单位圆中，生成随机点，在根据位于圆内的点数与总数之比来估算。
# 串行
首先完成在 CPU 上的串行代码以供对比验证，代码如下（仅关键函数）
// 蒙特卡罗模拟，串行
double monteCarloSerial(size_t& length) {
double x = 0;
double y = 0;
uint64_t count = 0;
// 随机数生成
uniform_real_distribution<double> dis(0.0, 1.0);
random_device rd;
mt19937_64 gen(rd());
for (size_t i = 0; i < length; i++)
{
x = dis(gen);
y = dis(gen);
if (hypot(x, y) < 1.0) // hypot(x,y) 返回 x,y 的和的平方根
count++;
}
return 4.0 * count / (long double)length;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 并行
# openMP CPU
利用 openMP 实现线程，代码如下（仅关键函数）
double monteCarloOpenMP_C(size_t& length) {
double x = 0;
double y = 0;
uint64_t count = { 0 };
#pragma omp parallel num_threads(4) // 4 是我的本机 CPU 物理核心数
{
// 随机数生成
uniform_real_distribution<double> dis(0.0, 1.0);
random_device rd;
mt19937_64 gen(rd());
#pragma omp for reduction(+:count) private(x, y)
for (int j = 0; j < length; j++)
{
x = dis(gen);
y = dis(gen);
if (hypot(x, y) < 1.0)
count++;
}
}
return 4.0 * count / (long double)length;
}
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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# CUDA
利用 CUDA 11.0 模拟，不过对于线程和线程块的数量设置还不知道如何可以做出比较成熟的优化。
double monteCarloCUDA(size_t& length) {
dim3 blocks(4096, 1, 1);
dim3 threads(128, 1, 1);
unsigned int num_dots_pt = length;
unsigned int num_total = threads.x * blocks.x;
cout << endl << "随机点数：" << num_total * num_dots_pt << endl << endl;
cout << "\tThe true value of pi: "
<< setiosflags(std::ios::fixed) << setprecision(10) << PI << endl;
unsigned int* count_h = new unsigned int[num_total];
unsigned int* count_d;
cudaMalloc(&count_d, num_total * sizeof(unsigned int)); // 用于保存结果的显存
curandState_t* randomState_d;
cudaMalloc(&randomState_d, num_total * sizeof(curandState_t)); // 用于保存随机信息的显存
monteCarloCUDAKernel <<<blocks, threads >>> (count_d, num_dots_pt, randomState_d);
cudaMemcpy(count_h, count_d, num_total * sizeof(unsigned int), cudaMemcpyDeviceToHost);
unsigned int count = 0;
for (unsigned int i = 0; i < num_total; i++)
count += count_h[i];
cudaFree(count_d);
cudaFree(randomState_d);
return 4.0 * (double)count / num_total / num_dots_pt;
}
__global__ void
monteCarloCUDAKernel(unsigned int* g_data, unsigned int length, curandState_t* states) {
unsigned int idx = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int count = length;
float x = 0, y = 0;
curand_init(idx, 0, 0, &states[idx]);
for (unsigned int i = 0; i < length; i++)
{
x = curand_uniform(&states[idx]);
y = curand_uniform(&states[idx]);
int temp = floorf(hypotf(x, y));
count = count - temp;
}
g_data[idx] = count;
}
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
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
# 编译及运行
# 编译
nvcc -Xcompiler -fopenmp -O3 Pi_OpenMP_CUDA.cu -o Pi_OpenMP_CUDA.out
1
- -Xcompiler -fopenmp : 开启 OpenMP
- -O3 : 启用优化
# 运行
./Pi_OpenMP_CUDA.out <Number>
1
模拟点数因子：实际计算值为 `Number * 4096 * 128`，4096 和 128 分别为 CUDA 线程块和线程数目预设参数
# 结果展示
$ ./Pi_OpenMP_CUDA.out 1e2
随机点数：52428800
The true value of pi: 3.1415926536
CUDA : The simulated value of pi: 3.1413547516 Relative error: 0.007573% Takes 399.552352 ms
Serial : The simulated value of pi: 3.1413627625 Relative error: 0.007318% Takes 2115.599377 ms
OpenMP CPU : The simulated value of pi: 3.1415091705 Relative error: 0.002657% Takes 531.666598 ms
1
2
3
4
5
6
7
8
9
2
3
4
5
6
7
8
9
$ ./Pi_OpenMP_CUDA.out 1e3
随机点数：524288000
The true value of pi: 3.1415926536
CUDA : The simulated value of pi: 3.1416261444 Relative error: 0.001066% Takes 432.088711 ms
Serial : The simulated value of pi: 3.1416330643 Relative error: 0.001286% Takes 21494.102259 ms
OpenMP CPU : The simulated value of pi: 3.1416608887 Relative error: 0.002172% Takes 5326.408407 ms
1
2
3
4
5
6
7
8
9
2
3
4
5
6
7
8
9
# 完整代码
# TODO
- 完善代码
- 添加 openMP GPU 示例
- 添加 tpp 示例
- 添加 MPI 示例