GPUNvidia CUDA 编程高级教程——利用蒙特卡罗法求解 的近似值

Posted 2022-12-07 从善若水

博主未授权任何人或组织机构转载博主任何原创文章，感谢各位对原创的支持！
博主链接

本人就职于国际知名终端厂商，负责modem芯片研发。
在5G早期负责终端数据业务层、核心网相关的开发工作，目前牵头6G算力网络技术标准研究。


博客内容主要围绕：
       5G/6G协议讲解
       算力网络讲解（云计算，边缘计算，端计算）
       高级C语言讲解
       Rust语言讲解

利用蒙特卡罗法求解 𝜋 的近似值

算法简介

估算 𝜋 有一个著名的技巧，那就是在单位面积内随机选择大量点，并计算落在单位圆内的点数。因为正方形的面积是 1，圆的面积是 𝜋/4 ，所以落在圆上的点的点数乘以 4，就是一个 𝜋 的良好近似值。

高度可并行

从并行编程的角度来看，该算法的良好特征之一是每个随机点都可以独立计算。我们只需要知道一个点的坐标，即可评估其是否落在圆内，对于点坐标 (𝑥,𝑦) 而言，如果 $$x^2+y^2<=1$$，那么点落在圆内，只要我们能处理好与计数器相关的任何竞态条件，表示圆内点数的计数器就可以递增。

单一 GPU 实现

我们来看看在单 GPU 上的 CUDA 中的实现情况。我们已提供实现情况的示例，如下所示：

#include <iostream>
#include <curand_kernel.h>

#define N 1024*1024

__global__ void calculate_pi(int* hits) 
    int idx = threadIdx.x + blockIdx.x * blockDim.x;

    // 初始化随机数状态（网格中的每个线程不得重复）
    int seed = 0;
    int offset = 0;
    curandState_t curand_state;
    curand_init(seed, idx, offset, &curand_state);

    // 在 (0.0, 1.0] 内生成随机坐标
    float x = curand_uniform(&curand_state);
    float y = curand_uniform(&curand_state);

    // 如果这一点在圈内，增加点击计数器
    if (x * x + y * y <= 1.0f) 
        atomicAdd(hits, 1);
    



int main(int argc, char** argv) 
    // 分配主机和设备值
    int* hits;
    hits = (int*) malloc(sizeof(int));

    int* d_hits;
    cudaMalloc((void**) &d_hits, sizeof(int));

    // 初始化点击次数并复制到设备
    *hits = 0;
    cudaMemcpy(d_hits, hits, sizeof(int), cudaMemcpyHostToDevice);

    // 启动核函数进行计算
    int threads_per_block = 256;
    int blocks = (N + threads_per_block - 1) / threads_per_block;

    calculate_pi<<<blocks, threads_per_block>>>(d_hits);
    cudaDeviceSynchronize();

    // 将最终结果复制回主机
    cudaMemcpy(hits, d_hits, sizeof(int), cudaMemcpyDeviceToHost);

    // 计算 pi 的最终值
    float pi_est = (float) *hits / (float) (N) * 4.0f;

    // 打印结果
    std::cout << "Estimated value of pi = " << pi_est << std::endl;
    std::cout << "Error = " << std::abs((M_PI - pi_est) / pi_est) << std::endl;

    // 清理
    free(hits);
    cudaFree(d_hits);

    return 0;

请注意，此代码仅用于指导目的，并不代表具有特别高的性能。具体原因如下：

• 我们将使用设备侧 API（属于cuRAND），直接在核函数中生成随机数。即使您不熟悉 cuRAND 也无妨，只需知道每个 CUDA 线程都有各自唯一的随机数即可。
• 我们让每个线程只计算一个值，所以计算强度很低。
• 在更改hits（“命中”）计数器时，我们将遇到许多原子操作的冲突。

即便如此，我们仍可用 100 万个样本点快速估算 𝜋 。与正确值相比，我们的计算误差应该仅约为 0.05%。

运行结果如下：

Estimated value of pi = 3.14319
Error = 0.000507708
CPU times: user 51.4 ms, sys: 16.5 ms, total: 67.9 ms
Wall time: 3.17 s

扩展到多个 GPU

有一个简单的方法可以将我们的示例扩展到多个 GPU，那就是使用管理多个 GPU 的单一主机进程。如果我们利用 M 个 GPU 对 N 个采样点进行计算，则可以将N/M采样点分配给每个 GPU，原则上可以M 倍地加快计算。

为了实施这一方法，我们要：

• 使用cudaGetDeviceCount确定可用 GPU 的数量。
• 以GPU数量为循环次数，在每次循环中使用cudaSetDevice指定执行代码的是哪个GPU。
• 在指定的 GPU 上执行分配给它的那部分工作。

  int device_count;
  cudaGetDeviceCount(&device_count);
  
  for (int i = 0; i < device_count; ++i) 
      cudaSetDevice(i);
      # Do single GPU worth of work.
  
  

代码实现

请注意，在此示例中，我们会给每个 GPU 一个不同的随机数生成器种子，以便每个 GPU 进行不同的工作。因此，我们的答案会有所不同。

#include <iostream>
#include <curand_kernel.h>

#define N 1024*1024

__global__ void calculate_pi(int* hits, int device) 
    int idx = threadIdx.x + blockIdx.x * blockDim.x;

    // 初始化随机数状态（网格中的每个线程不得重复）
    int seed = device;
    int offset = 0;
    curandState_t curand_state;
    curand_init(seed, idx, offset, &curand_state);

    // 在 (0.0, 1.0] 内生成随机坐标
    float x = curand_uniform(&curand_state);
    float y = curand_uniform(&curand_state);

    // 如果这一点在圈内，增加点击计数器
    if (x * x + y * y <= 1.0f) 
        atomicAdd(hits, 1);
    



int main(int argc, char** argv) 
    // 确定 GPU 数量
    int device_count;
    cudaGetDeviceCount(&device_count);

    std::cout << "Using " << device_count << " GPUs" << std::endl;

    // 分配主机和设备值（每个 GPU 一个）
    int** hits = (int**) malloc(device_count * sizeof(int*));
    for (int i = 0; i < device_count; ++i) 
        hits[i] = (int*) malloc(sizeof(int));
    

    int** d_hits = (int**) malloc(device_count * sizeof(int*));
    for (int i = 0; i < device_count; ++i) 
        cudaSetDevice(i);
        cudaMalloc((void**) &d_hits[i], sizeof(int));
    

    // 初始化点击次数并复制到设备
    for (int i = 0; i < device_count; ++i) 
        *hits[i] = 0;
        cudaSetDevice(i);
        cudaMemcpy(d_hits[i], hits[i], sizeof(int), cudaMemcpyHostToDevice);
    

    // 启动核函数进行计算
    int threads_per_block = 256;
    int blocks = (N / device_count + threads_per_block - 1) / threads_per_block;

    // 先启动所有核函数，以支持异步执行
    // 然后在所有设备上同步。
    for (int i = 0; i < device_count; ++i) 
        cudaSetDevice(i);
        calculate_pi<<<blocks, threads_per_block>>>(d_hits[i], i);
    

    for (int i = 0; i < device_count; ++i) 
        cudaSetDevice(i);
        cudaDeviceSynchronize();
    

    // 将最终结果复制回主机
    for (int i = 0; i < device_count; ++i) 
        cudaSetDevice(i);
        cudaMemcpy(hits[i], d_hits[i], sizeof(int), cudaMemcpyDeviceToHost);
    

    // 计算所有设备的点击总数
    int hits_total = 0;
    for (int i = 0; i < device_count; ++i) 
        hits_total += *hits[i];
    

    // 计算 pi 的最终值
    float pi_est = (float) hits_total / (float) (N) * 4.0f;

    // 打印结果
    std::cout << "Estimated value of pi = " << pi_est << std::endl;
    std::cout << "Error = " << std::abs((M_PI - pi_est) / pi_est) << std::endl;

    // 清理
    for (int i = 0; i < device_count; ++i) 
        free(hits[i]);
        cudaFree(d_hits[i]);
    
    free(hits);
    free(d_hits);

    return 0;

运行结果

我这边使用了4个GPU，根据测试环境不同，我们显示的结果也是不同的。

Using 4 GPUs
Estimated value of pi = 3.14072
Error = 0.000277734
CPU times: user 27.2 ms, sys: 16.1 ms, total: 43.3 ms
Wall time: 2.46 s

---