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

Minimize the number of cuRAND states stored during a MC simulation

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

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

我知道一种可能的方法是以 grid stride loop 方式处理每个线程的多个数组元素。在这里我想研究一个不同的方法。

我知道对于我的目标计算能力 (3.5),每个 SM 的最大常驻线程数是 2048 中的 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;
}

根据您使用的 RNG 类型,有几种可能的方法。在您的情况下,如果您对 rng 类型有一点自由,并且如果您的设备具有良好的双精度,则您可能希望完全删除存储 rng 状态的概念并调用 init 和 skip ahead 技术。唯一需要的元素是种子和可以根据迭代和模拟 ID 计算的索引。

参见 Skip Ahead part of cuRand documentation, and see that most curand_init method 接受一个偏移量参数。在某些情况下,考虑到 RNG 状态结构的性质和 init 的小成本,最好调用 cuda_init 并在可能驻留在寄存器 space 中的状态数据结构上使用适当的偏移量比 load/store 每个随机值提取时来自全局内存的状态数据结构。

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

我从你的代码开始 shell 并将其变成经典的 "hello world" MC 问题:根据正方形面积与内切圆面积之比随机计算 pi生成点以估计面积。

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

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

  2. 创建一个更小的一维网格,以及每个线程的状态,并让每个线程计算多个随机数points/tests,使得相同数量的points/tests与情况 1 一样生成。

  3. 创建一个与方法 2 大小相同的网格,但也创建一个 "unique resident thread index",然后提供足够的状态来覆盖唯一的常驻线程。每个线程将有效地使用使用 "resident thread index" 提供的状态,而不是普通的全局唯一线程索引。计算与情况 1 中相同的 points/tests 个数。

"unique resident thread index"的计算可不是小事。在主机代码中:

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

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

在设备代码中:

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

从结果的角度来看,可以做出以下观察:

  1. 与其他答案中所述相反,初始化时间 并非 无关紧要。它不应该被忽视。对于这个简单的问题,init kernel 运行-time 超过了计算kernel 运行 的时间。因此,我们最重要的收获是我们应该寻找方法来最小化随机生成器状态的创建。我们当然不想无谓地重新进行运行的初始化。因此,对于这个特定的 code/test.

  2. ,我们可能应该根据这些结果放弃方法 1
  3. 从计算内核运行时间的角度来看,没有什么比另一个内核更好的了。

  4. 从代码复杂度的角度来看,方法 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,给出了概念的附加描述级别