最小化 MC 模拟期间存储的 cuRAND 状态数量

2024-05-19

我目前正在 CUDA 中编写蒙特卡罗模拟。因此,我需要生成lots使用随机数cuRAND图书馆。 每个线程处理一个巨大的元素floatarray(示例中省略)并在每次内核调用时生成 1 或 2 个随机数。

通常的方法(参见下面的示例)似乎是为每个线程分配一个状态。这些状态将在后续内核调用中重用。 然而,当线程数量增加(在我的应用程序中高达 10⁸)时,这不能很好地扩展,因为它成为程序中的主要内存(和带宽)使用。

我知道一种可能的方法是在一个线程中处理多个数组元素网格步幅循环 https://devblogs.nvidia.com/parallelforall/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/时尚。在这里我想研究一种不同的方法。

我知道对于我的目标计算能力 (3.5),每个 SM 的最大驻留线程数是2048 https://en.wikipedia.org/wiki/CUDA#Version_features_and_specifications, i.e.下例中的 2 个块。
无论线程总数如何,每个多处理器是否可以仅使用 2048 个状态?生成的所有随机数仍应在统计上独立。

我认为如果每个驻留线程都关联一个以 2048 为模的唯一索引,则可以完成此操作,然后可以使用该索引从数组中获取状态。这样的索引存在吗?

更普遍,还有其他方法可以减少 RNG 状态的内存占用吗?

示例代码(每个线程一个状态)

#include <cuda.h>
#include <curand.h>
#include <curand_kernel.h>
#include <assert.h>

#define gridSize 100
#define blockSize 1024
#define iter 100

__global__ void rng_init(unsigned long long seed, curandState * states) {
  const size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
  curand_init(seed, Idx, 0, &states[Idx]);
}

__global__ void kernel(curandState * states) {
  const size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
  const float x = curand_uniform(&states[Idx]);
  // Do stuff with x ...
}

int main() {
  curandState * states;
  assert(cudaMalloc(&states, gridSize*blockSize*sizeof(curandState)) == cudaSuccess);
  rng_init<<<gridSize,blockSize>>>(clock(), states);
  assert(cudaDeviceSynchronize() == cudaSuccess);

  for (size_t i = 0 ; i < iter ; ++i) {
    kernel<<<gridSize,blockSize>>>(states);
    assert(cudaDeviceSynchronize() == cudaSuccess);
  }
  return 0;
}

简而言之,在您的问题中,您提到使用网格跨步循环作为折叠所需状态的方法。我认为这种方法或类似的方法(如@talonmies建议的)是最明智的方法。选择一种方法,将线程数减少到保持机器繁忙/充分利用所需的最低限度。然后让每个线程计算多个操作,重新使用提供的随机生成器状态。

我从你的代码 shell 开始,把它变成了一个经典的“hello world”MC 问题:根据正方形面积与内切圆面积之比计算 pi,随机生成点来估计面积。

那么我们将考虑3种方法:

  1. 创建一个大型一维网格以及每个线程的状态,以便每个线程计算一个随机点并对其进行测试。

  2. 创建一个更小的一维网格以及每个线程的状态,并允许每个线程计算多个随机点/测试,以便生成与情况 1 相同数量的点/测试。

  3. 创建一个与方法 2 大小相同的网格,同时创建一个“唯一驻留线程索引”,然后提供足够的状态来覆盖唯一驻留线程。每个线程将有效地使用通过“常驻线程索引”提供的状态,而不是普通的全局唯一线程索引。计算与情况 1 相同数量的点/测试。

“唯一驻留线程索引”的计算并不是一件小事。在主机代码中:

  1. 我们必须确定理论上可以驻留的最大块数。我为此使用了一个简单的启发式方法,但可以说还有更好的方法。我只需将每个多处理器的最大驻留线程数除以每个块所选线程的数量。为此我们必须使用整数除法。

  2. 然后初始化足够的状态以覆盖最大块数乘以 GPU 上的 SM 数。

在设备代码中:

  1. 确定我在哪个 SM 上。
  2. 使用该 SM,执行原子测试以检索每个 SM 的唯一块号。
  3. 将 1 和 2 的结果相结合,再加上threadIdx.x变量来创建“唯一驻留线程索引”。这将成为我们典型的线程索引,而不是通常的计算。

从结果的角度来看,可以得出以下结论:

  1. 与其他答案中所述相反,初始化时间是not微不足道。它不应该被忽视。对于这个简单的问题,初始化内核运行时间超过了计算内核运行时间。因此,我们最重要的结论是,我们应该寻求方法来最小化随机生成器状态的创建。我们当然不想不必要地重新运行初始化。因此,对于这个特定的代码/测试,我们可能应该根据这些结果放弃方法 1。

  2. 从计算内核运行时的角度来看,一个内核没有什么值得称赞的。

  3. 从代码复杂度的角度来看,方法 2 显然没有方法 3 复杂,但性能却大致相同。

对于这个特定的测试用例,方法 2 似乎是获胜者,正如 @talonmies 所预测的那样。

下面是这 3 种方法的一个工作示例。我并不是说这对每个案例、代码或场景都没有缺陷或具有指导意义。这里有很多变化的部分,但我相信上面的 3 个结论对于这个案例是有效的。

$ cat t1201.cu
#include <cuda.h>
#include <curand.h>
#include <curand_kernel.h>
#include <assert.h>
#include <iostream>
#include <stdlib.h>

#define blockSize 1024

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

#define MAX_SM 64

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}


__device__ unsigned long long count = 0;
__device__ unsigned int blk_ids[MAX_SM] = {0};

__global__ void rng_init(unsigned long long seed, curandState * states) {
  const size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
  curand_init(seed, Idx, 0, &states[Idx]);
}

__global__ void kernel(curandState * states, int length) {
  const size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
  for (int i = 0; i < length; i++){
    const float x = curand_uniform(&states[Idx]);
    const float y = curand_uniform(&states[Idx]);
    if (sqrtf(x*x+y*y)<1.0)
      atomicAdd(&count, 1ULL);}
}

static __device__ __inline__ int __mysmid(){
  int smid;
  asm volatile("mov.u32 %0, %%smid;" : "=r"(smid));
  return smid;}

__device__ int get_my_resident_thread_id(int sm_blk_id){
  return __mysmid()*sm_blk_id + threadIdx.x;
}

__device__ int get_block_id(){
  int my_sm = __mysmid();
  int my_block_id = -1;
  bool done = false;
  int i = 0;
  while ((!done)&&(i<32)){
    unsigned int block_flag = 1<<i;
    if ((atomicOr(blk_ids+my_sm, block_flag)&block_flag) == 0){my_block_id = i; done = true;}
    i++;}
  return my_block_id;
}

__device__ void release_block_id(int block_id){
  unsigned int block_mask = ~(1<<block_id);
  int my_sm = __mysmid();
  atomicAnd(blk_ids+my_sm, block_mask);
}

__global__ void kernel2(curandState * states, int length) {

  __shared__ volatile int my_block_id;
  if (!threadIdx.x) my_block_id = get_block_id();
  __syncthreads();
  const size_t Idx = get_my_resident_thread_id(my_block_id);
  for (int i = 0; i < length; i++){
    const float x = curand_uniform(&states[Idx]);
    const float y = curand_uniform(&states[Idx]);
    if (sqrtf(x*x+y*y)<1.0)
      atomicAdd(&count, 1ULL);}
  __syncthreads();
  if (!threadIdx.x) release_block_id(my_block_id);
  __syncthreads();
}



int main(int argc, char *argv[]) {
  int gridSize = 10;
  if (argc > 1) gridSize = atoi(argv[1]);
  curandState * states;
  assert(cudaMalloc(&states, gridSize*gridSize*blockSize*sizeof(curandState)) == cudaSuccess);
  unsigned long long hcount;
  //warm - up
  rng_init<<<gridSize*gridSize,blockSize>>>(1234ULL, states);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  //method 1: 1 curand state per point
  std::cout << "Method 1 init blocks: " << gridSize*gridSize << std::endl;
  unsigned long long dtime = dtime_usec(0);
  rng_init<<<gridSize*gridSize,blockSize>>>(1234ULL, states);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  unsigned long long initt = dtime_usec(dtime);
  kernel<<<gridSize*gridSize,blockSize>>>(states, 1);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  dtime = dtime_usec(dtime) - initt;
  cudaMemcpyFromSymbol(&hcount, count, sizeof(unsigned long long));
  std::cout << "method 1 elapsed time: " << dtime/(float)USECPSEC << " init time: " << initt/(float)USECPSEC << " pi: " << 4.0f*hcount/(float)(gridSize*gridSize*blockSize) << std::endl;
  hcount = 0;
  cudaMemcpyToSymbol(count, &hcount, sizeof(unsigned long long));
  //method 2: 1 curand state per gridSize points
  std::cout << "Method 2 init blocks: " << gridSize << std::endl;
  dtime = dtime_usec(0);
  rng_init<<<gridSize,blockSize>>>(1234ULL, states);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  initt = dtime_usec(dtime);
  kernel<<<gridSize,blockSize>>>(states, gridSize);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  dtime = dtime_usec(dtime) - initt;
  cudaMemcpyFromSymbol(&hcount, count, sizeof(unsigned long long));
  std::cout << "method 2 elapsed time: " << dtime/(float)USECPSEC << " init time: " << initt/(float)USECPSEC << " pi: " << 4.0f*hcount/(float)(gridSize*gridSize*blockSize) << std::endl;
  hcount = 0;
  cudaMemcpyToSymbol(count, &hcount, sizeof(unsigned long long));
  //method 3: 1 curand state per resident thread
  // compute the maximum number of state entries needed
  int num_sms;
  cudaDeviceGetAttribute(&num_sms, cudaDevAttrMultiProcessorCount, 0);
  int max_sm_threads;
  cudaDeviceGetAttribute(&max_sm_threads, cudaDevAttrMaxThreadsPerMultiProcessor, 0);
  int max_blocks = max_sm_threads/blockSize;
  int total_state = max_blocks*num_sms*blockSize;
  int rgridSize = (total_state + blockSize-1)/blockSize;
  std::cout << "Method 3 sms: " << num_sms << " init blocks: " << rgridSize << std::endl;
  // run test
  dtime = dtime_usec(0);
  rng_init<<<rgridSize,blockSize>>>(1234ULL, states);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  initt = dtime_usec(dtime);
  kernel2<<<gridSize,blockSize>>>(states, gridSize);
  assert(cudaDeviceSynchronize() == cudaSuccess);
  dtime = dtime_usec(dtime) - initt;
  cudaMemcpyFromSymbol(&hcount, count, sizeof(unsigned long long));
  std::cout << "method 3 elapsed time: " << dtime/(float)USECPSEC << " init time: " << initt/(float)USECPSEC << " pi: " << 4.0f*hcount/(float)(gridSize*gridSize*blockSize) << std::endl;
  hcount = 0;
  cudaMemcpyToSymbol(count, &hcount, sizeof(unsigned long long));
  return 0;
}

$ nvcc -arch=sm_35 -O3 -o t1201 t1201.cu
$ ./t1201 28
Method 1 init blocks: 784
method 1 elapsed time: 0.001218 init time: 3.91075 pi: 3.14019
Method 2 init blocks: 28
method 2 elapsed time: 0.00117 init time: 0.003629 pi: 3.14013
Method 3 sms: 14 init blocks: 28
method 3 elapsed time: 0.001193 init time: 0.003622 pi: 3.1407
$

对于方法 3,给出了概念的附加级别描述here https://stackoverflow.com/a/59110744/1695960

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

最小化 MC 模拟期间存储的 cuRAND 状态数量 的相关文章

  • 在 k 个 bin 之间随机分配一个整数

    我正在寻找一个高效的Python函数 它可以随机分配一个整数k垃圾箱 也就是说 某个函数allocate n k 将产生一个k sized整数数组总和为n 例如 allocate 4 3 可以产生 4 0 0 0 2 2 1 2 1 etc
  • MySQL:你能指定一个随机限制吗?

    有没有办法在 SQL MySQL 中随机化限制数字 我希望能够做的是在查询中获取随机数量的结果以在插入子查询中使用 而无需任何服务器端脚本 我希望能够作为假设说明运行的查询是 SELECT id FROM users ORDER BY RA
  • 如何在C++11中让每个线程使用自己的RNG

    我正在 C 11 中使用新的随机数生成器 虽然有不同的看法 但由此看来thread https stackoverflow com questions 8813592 c11 thread safety of random number g
  • 模行为背后的数学

    Preamble 这个问题与 P RNG 的行为无关rand 它是关于使用均匀分布的两个值的幂对模 介绍 我知道不应该使用模数 将一个值从一个范围转换为另一个范围 例如从 0 到 5 之间的值rand 功能 会有偏差 这里解释了https
  • CUDA 模型 - 什么是扭曲尺寸?

    最大工作组大小和扭曲大小之间有什么关系 假设我的设备有 240 个 CUDA 流处理器 SP 并返回以下信息 CL DEVICE MAX COMPUTE UNITS 30 CL DEVICE MAX WORK ITEM SIZES 512
  • 固定长度的随机数

    我想生成一个 0 9 数字且长度 5 的随机整数 我尝试这样做 function genRand min max for var i 1 i lt 5 i var range max min 1 return Math floor Math
  • 使用 GPU 进行 Matlab 卷积

    我用gpuArray尝试了matlab的卷积函数conv2 convn 例如 convn gpuArray rand 100 100 10 single gpuArray rand 5 single 并将其与 cpu 版本 convn ra
  • 使用 cudamalloc()。为什么是双指针?

    我目前正在浏览有关的教程示例http code google com p stanford cs193g sp2010 http code google com p stanford cs193g sp2010 学习CUDA 演示的代码 g
  • 第一个随机数始终小于其余随机数

    我碰巧注意到 在 C 中 使用 std rand 方法调用的第一个随机数大多数时候都明显小于第二个随机数 关于 Qt 实现 第一个几乎总是小几个数量级 qsrand QTime currentTime msec qDebug lt lt q
  • 在 R 中,如何让 PRNG 在平台之间给出相同的浮点数?

    在 R 4 1 1 中运行以下代码会在平台之间产生不同的结果 set seed 1 x lt rnorm 3 3 print x 22 0 83562861241004716 intel windows 0 8356286124100471
  • java代码的等效vb代码

    谁能告诉我这段Java代码到底做了什么 SecureRandom random SecureRandom getInstance SHA1PRNG byte bytes new byte 20 synchronized random ran
  • 使用 Cuda 并行读取多个文本文件

    我想使用 CUDA 在多个文件中并行搜索给定字符串 我计划使用 pfac 库来搜索给定的字符串 问题是如何并行访问多个文件 示例 我们有一个包含 1000 个文件的文件夹 需要搜索 这里的问题是我应该如何访问给定文件夹中的多个文件 应该动态
  • 如何获取常量内存中的统计数据

    我有一个函数 它会创建一些随机的数值结果 我知道 结果将是 a 小 a b 约 50 范围内的整数a b 我想创建一个执行上述函数 1000000 次的函数 并计算每个结果出现的频率 该函数使用随机生成器来生成结果 问题是 我不知道如何在常
  • CUDA 代码会损坏 GPU 吗?

    在测试包含内存错误的 CUDA 时 我的屏幕被冻结了 重新启动后我无法再检测到显卡 我的代码是否有可能物理损坏该卡 这发生在 Ubuntu 14 04 下 我不知道该卡的型号 因为我无法检测到它 但我记得它是一张相当新的卡 感谢所有的评论我
  • cudaMalloc使用向量>进行管理 > C++ - NVIDIA CUDA

    我正在通过 NVIDIA GeForce GT 650M GPU 为我创建的模拟实现多线程 为了确保一切正常工作 我创建了一些辅助代码来测试一切是否正常 在某一时刻 我需要更新变量向量 它们都可以单独更新 这是它的要点 device int
  • 内存高效的随机数迭代器,无需替换

    我觉得这应该很容易 但经过多次搜索和尝试后我找不到答案 基本上 我有大量的物品 我想以随机顺序进行采样 而不需要更换 在本例中 它们是二维数组中的单元 我用于较小数组的解决方案不会转换 因为它需要对内存数组进行改组 如果我必须采样的数量很小
  • rand() 的实现

    我正在用 C 编写一些嵌入式代码 需要使用 rand 函数 不幸的是 控制器的库不支持 rand 我需要一个快速的简单实现 但更重要的是空间开销很小 可以产生相对高质量的随机数 有谁知道使用哪种算法或示例代码 编辑 它用于图像处理 因此 相
  • 一维纹理内存访问比一维全局内存访问更快吗?

    我正在测量标准纹理和 1Dtexture 内存访问之间的差异 为此 我创建了两个内核 global void texture1D float doarray int size int index calculate each thread
  • CUDA 中的 JPEG 库

    我正在尝试在 CUDA 中压缩和解压缩图像 到目前为止我已经找到了这个库 http sourceforge net projects cuj2k source navbar http sourceforge net projects cuj
  • 如何从尖点库矩阵格式获取原始指针

    我需要从尖点库矩阵格式获取原始指针 例如 cusp coo matrix

随机推荐