CUDA 的 Mersenne Twister，用于任意数量的线程

Question

CUDA 对 Mersenne Twister （MT） 随机数生成器的实现仅限于 256 和 200 块/网格的最大线程/块数，即最大线程数为 51200 。

因此，无法启动使用 MT 的内核

kernel<<<blocksPerGrid, threadsPerBlock>>>(devMTGPStates, ...)

哪里

int blocksPerGrid = (n+threadsPerBlock-1)/threadsPerBlock;

n是线程总数。

使用MT进行threads > 51200的最佳方法是什么？

我的方法是否对blocksPerGrid和threadsPerBlock使用常量值，例如 在内核代码中<<<128,128>>>并使用以下代码：

__global__ void kernel(curandStateMtgp32 *state, int n, ...) { 
    int id = threadIdx.x+blockIdx.x*blockDim.x;
    while (id < n) {
        float x = curand_normal(&state[blockIdx.x]);
        /* some more calls to curand_normal() followed
           by the algorithm that works with the data */
        id += blockDim.x*gridDim.x; 
    }
}

我不确定这是否是正确的方式，或者它是否会以不希望的方式影响 MT 状态？

谢谢。

Answer 1

我建议您仔细阅读CULAND文档。

当每个块使用 256 个线程（最多 64 个块来生成数字）时，MT API 将最有效。

如果您需要更多，您有多种选择：

1. 只需从现有状态生成更多数字 - set （即 64块，256 个线程），并将这些数字分配给需要它们的线程。
2. 每个块使用多个状态（但这不允许你超过状态集中的总体限制，它只是解决了对单个块的需求。
3. 创建多个具有独立种子（因此具有独立状态集）的 MT 生成器。

一般来说，我认为您概述的内核没有问题，并且与上面的选择 1 大致一致。 但是，它不允许超过 51200 个线程。 （您的示例有 <<<128, 128>>> 所以 16384 个线程）

Answer 2

根据罗伯特的回答，下面我提供了一个完整的示例，介绍如何将 cuRAND 的 Mersenne Twister 用于任意数量的线程。我正在使用 Robert 的第一个选项从现有状态集生成更多数字，并将这些数字分配给需要它们的线程。

// --- Generate random numbers with cuRAND's Mersenne Twister
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cuda.h>
#include <curand_kernel.h>
/* include MTGP host helper functions */
#include <curand_mtgp32_host.h>
#define BLOCKSIZE   256
#define GRIDSIZE    64
/*******************/
/* GPU ERROR CHECK */
/*******************/
#define gpuErrchk(x) do { if((x) != cudaSuccess) { 
    printf("Error at %s:%dn",__FILE__,__LINE__); 
    return EXIT_FAILURE;}} while(0)
#define CURAND_CALL(x) do { if((x) != CURAND_STATUS_SUCCESS) { 
    printf("Error at %s:%dn",__FILE__,__LINE__); 
    return EXIT_FAILURE;}} while(0)
/*******************/
/* iDivUp FUNCTION */
/*******************/
__host__ __device__ int iDivUp(int a, int b) { return ((a % b) != 0) ? (a / b + 1) : (a / b); }
/*********************/
/* GENERATION KERNEL */
/*********************/
__global__ void generate_kernel(curandStateMtgp32 * __restrict__ state, float * __restrict__ result, const int N)
{
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    for (int k = tid; k < N; k += blockDim.x * gridDim.x)
        result[k] = curand_uniform(&state[blockIdx.x]);
}
/********/
/* MAIN */
/********/
int main()
{
    const int N = 217 * 123;
    // --- Allocate space for results on host
    float *hostResults = (float *)malloc(N * sizeof(float));
    // --- Allocate and initialize space for results on device 
    float *devResults; gpuErrchk(cudaMalloc(&devResults, N * sizeof(float)));
    gpuErrchk(cudaMemset(devResults, 0, N * sizeof(float)));
    // --- Setup the pseudorandom number generator
    curandStateMtgp32 *devMTGPStates; gpuErrchk(cudaMalloc(&devMTGPStates, GRIDSIZE * sizeof(curandStateMtgp32)));
    mtgp32_kernel_params *devKernelParams; gpuErrchk(cudaMalloc(&devKernelParams, sizeof(mtgp32_kernel_params)));
    CURAND_CALL(curandMakeMTGP32Constants(mtgp32dc_params_fast_11213, devKernelParams));
    //CURAND_CALL(curandMakeMTGP32KernelState(devMTGPStates, mtgp32dc_params_fast_11213, devKernelParams, GRIDSIZE, 1234));
    CURAND_CALL(curandMakeMTGP32KernelState(devMTGPStates, mtgp32dc_params_fast_11213, devKernelParams, GRIDSIZE, time(NULL)));
    // --- Generate pseudo-random sequence and copy to the host
    generate_kernel << <GRIDSIZE, BLOCKSIZE >> >(devMTGPStates, devResults, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());
    gpuErrchk(cudaMemcpy(hostResults, devResults, N * sizeof(float), cudaMemcpyDeviceToHost));
    // --- Print results
    //for (int i = 0; i < N; i++) {
    for (int i = 0; i < 10; i++) {
        printf("%fn", hostResults[i]);
    }
    // --- Cleanup
    gpuErrchk(cudaFree(devMTGPStates));
    gpuErrchk(cudaFree(devResults));
    free(hostResults);
    return 0;
}