多线程与GPU编程

年轻人的第一次多线程。

这个工程源自我朋友在群中的提问：

抛硬币连续抛出 5 次正面朝上的期望是多少？

我没有啥数学思路，于是就写了个C语言程序去算，结果趋向于6.79726，和一个朋友的C语言程序结果一致，但是和另一个朋友的Python程序结果不一致，他的结果是62。

通过网上所给出的方法与结果(对于 NN 次连续一面的问题，次数的期望是 2N−22^N - 2 )，不难知道 62 是正确结果，我们遂开始研究。

通过测试，C语言的rand()函数所得到的随机数，奇偶在大数量级是均匀分布的，但连续 N 个奇数或者偶数的概率略低于实际概率。换用 C++ 的 mt19937 随机数算法（梅森旋转），将实验的循环次数提高到一亿次，所得平均值基本可以收敛到小数点后两位，但时间较慢，使用 i9 12900H 计算需要花费2分钟/1亿次。

C++ CPU多线程

于是我去学了C++的多线程。

通过C++的<thread>库，我们可以很容易的进行多线程并行。使用 i9 12900H 进行20线程计算一亿次迭代将只花费 2.5 秒左右。

需要注意的是，修改全局变量时不要忘记线程锁，主线程使用该全局变量时，应先使用.join()方法同步线程。

编译命令: g++ calc.cpp -o calc.exe -Wall -std=c++11 -O2

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
#include<iostream>
#include<random>
#include<thread>
#include<mutex>
using namespace std;
const int N = 5000000, T = 20;
typedef unsigned long long ull;
random_device rd;
ull result = 0;
mutex mtx;
void func() {
    mt19937 r_eng(rd());
    ull res = 0;
    int cnt = 0, tot = 0;
    for(int i = 0; i < N; i++) {
        cnt = 0;
        tot = 0;
        while(cnt < 5) {
            ++tot;
            cnt = r_eng() % 2 ? cnt + 1 : 0;
        }
        res += 1ull * tot;
    }
    mtx.lock();
    result += res;
    mtx.unlock();
}
int main() {
    thread* t[T];
    for(int i = 0; i < T; i++){
        t[i] = new thread(func);
    }
    for(int i = 0; i < T; i++){
        t[i]->join();
        delete t[i];
    }
    cout << (double)1 * result / N / T << endl;
    return 0;
}

CUDA C++ GPU多线程

于是我去学了C++通过CUDA的GPU编程。

通过GPU的计算单元计算这种简单的运算，时间将会更快，使用最高150瓦的 RTX 3070 Ti 进行计算，一亿次只需要花费 0.32 秒左右，相比单线程C++程序提高了将近四百倍。

此处需要注意的是，这里GPU调用的函数中使用随机数需要使用CUDA提供的方法。并且，一味提高核心数降低每核心线程数好像并不能显著提高计算速度，现在代码中应该算较优配置数值。

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
#include<iostream>
#include<thread>
#include<cstdio>
#include"curand_kernel.h"
using namespace std;
typedef unsigned long long ull;
const int M = 100000000;
const int NPT = 10000;
const int N = M / NPT;
__global__ void setup_kernel(curandState* state, unsigned int seed, int NUM) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx > NUM) return;
    curand_init(seed, idx, 0, &state[idx]);
}
__global__ void func(ull* results, curandState* globalState, int NUM) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx >= NUM) return;
    curandState localState = globalState[idx];
    int cnt, tot;
    ull result = 0;
    for(int i = 0; i < NPT; i++) {
        cnt = 0, tot = 0;
        while(cnt < 5) {
            ++tot;
            cnt = curand(&localState) % 2 ? cnt + 1 : 0;
        }
        result += tot;
    }
    results[idx] += result;
}
int main() {
    size_t threads_per_block = 256;
    size_t number_of_blocks = (N + threads_per_block - 1) / threads_per_block;
    ull* results;
    curandState* devStates;
    cudaMalloc(&devStates, number_of_blocks * threads_per_block * sizeof(curandState));
    srand(time(NULL));
    cudaMallocManaged(&results, sizeof(ull) * N);
    setup_kernel<<<number_of_blocks, threads_per_block>>>(devStates, rand(), N);
    func<<<number_of_blocks, threads_per_block>>>(results, devStates, N);
    cudaDeviceSynchronize();
    ull result = 0;
    for(int i = 0; i < N; i++) result += results[i];
    cudaFree(results);
    cout << (double)1 * result / M << endl;
    return 0;
}

总结

现在能使用C++编写多线程程序与进行GPU编程了，还思考了伪随机数的问题，收获斐然。