curand_init的offset参数

Question

我在了解的offset参数时遇到问题curand_init。

该cuRAND指南说：

offset参数用于跳过序列。如果offset = 100，则生成的第一个随机数将是序列中的第100个。这允许同一程序的多次运行继续从相同序列生成结果，而不会出现重叠。

这似乎是图中所示的“先行者”概念

T. Bradley, J. du Toit, R. Tong, M. Giles, P. Woodhams, "Parallelization Techniques for Random Number Generators", GPU Computing Gems, Emerald Edition.

考虑以下代码：

#include <stdio.h>
#include <curand.h>
#include <curand_kernel.h>

#define BLOCKSIZE 256

/**********/
/* iDivUp */
/**********/
int iDivUp(int a, int b){ return ((a % b) != 0) ? (a / b + 1) : (a / b); }

/***********************/
/* CUDA ERROR CHECKING */
/***********************/
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
    if (code != cudaSuccess) 
    {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); 
        if (abort) exit(code);
    }
}

/********************************************************/
/* KERNEL FUNCTION FOR TESTING RANDOM NUMBER GENERATION */
/********************************************************/
__global__ void testrand1(unsigned long seed, float *a, int N){
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    curandState state;
    if (idx < N) {
        curand_init(seed, idx, 0, &state);
        a[idx] = curand_uniform(&state);
    }
}

/********/
/* MAIN */
/********/
int main() {

    const int N = 10;

    float *h_a  = (float*)malloc(N*sizeof(float));
    float *d_a; gpuErrchk(cudaMalloc((void**)&d_a, N*sizeof(float)));

    testrand1<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(1234, d_a, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaMemcpy(h_a, d_a, N*sizeof(float), cudaMemcpyDeviceToHost));

    for (int i=0; i<N; i++) printf("%i %f\n", i, h_a[i]);

    getchar();
}

运行此代码将生成：

0 0.145468
1 0.820181
2 0.550399
3 0.294830
4 0.914733
5 0.868979
6 0.321921
7 0.782857
8 0.011302
9 0.285450

现在，如果我使用

curand_init(seed, idx+2, 0, &state);

我收到：

0 0.550399
1 0.294830
2 0.914733
3 0.868979
4 0.321921
5 0.782857
6 0.011302
7 0.285450
8 0.781606
9 0.233840

这意味着生成的序列从与关联的随机数生成器（RNG）序列的第三个元素开始seed = 1234。根据cuRAND指南中的定义，我会说

curand_init(seed, idx+2, 0, &state);

应该等于

curand_init(seed, idx, 2, &state);

不幸的是，如果我使用上述行，则会收到：

0 0.870710
1 0.511765
2 0.782640
3 0.620706
4 0.554513
5 0.214082
6 0.118647
7 0.993959
8 0.104572
9 0.231619

与上述不同。

在11同一指南的页面上，讨论时curand_init()，其写为：

设置的状态将是从种子状态2^67 # sequence + offset调用后curand()的状态。

我将其解释为与特定种子相关offset的2^67可用序列之一中的偏移量。

任何人都可以帮助我了解offsetin 中参数的用法curand_init吗？

Answer 1

curand_init(seed, idx+2, 0, &state);

应该不等同于：

curand_init(seed, idx, 2, &state);

因为使用相同种子和不同序列号产生的序列将没有统计相关值。请注意，此语句没有基于offset参数的其他限定。

生成的整个序列的周期大于2 ^ 190。在该总体序列中，存在由curand_init调用的第二个参数标识的子序列。这些子序列彼此独立（并且在统计上不相关）。这些子序列在整个序列中大约相距2 ^ 67。偏移量参数（第3个参数）选择该子序列中的起始位置。由于offset参数不能大于2 ^ 67，因此无法单独使用offset参数使生成的数字重叠。

但是，如果选择的话，提供的skipahead_sequence 功能可以使我们做到这一点。考虑对代码和示例运行进行以下修改：

$ cat t557.cu
#include <stdio.h>
#include <curand.h>
#include <curand_kernel.h>

#define BLOCKSIZE 256

/**********/
/* iDivUp */
/**********/
int iDivUp(int a, int b){ return ((a % b) != 0) ? (a / b + 1) : (a / b); }

/***********************/
/* CUDA ERROR CHECKING */
/***********************/
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
    if (code != cudaSuccess)
    {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
        if (abort) exit(code);
    }
}

/********************************************************/
/* KERNEL FUNCTION FOR TESTING RANDOM NUMBER GENERATION */
/********************************************************/
__global__ void testrand1(unsigned long seed, float *a, int N){
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    curandState state;
    if (idx < N) {
        curand_init(seed, idx, 0, &state);
        a[(idx*2)] = curand_uniform(&state);
        if(idx%2)
          skipahead_sequence(1, &state);
        a[(idx*2)+1] = curand_uniform(&state);

    }
}

/********/
/* MAIN */
/********/
int main() {

    const int N = 10;

    float *h_a  = (float*)malloc(2*N*sizeof(float));
    float *d_a; gpuErrchk(cudaMalloc((void**)&d_a, 2*N*sizeof(float)));

    testrand1<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(1234, d_a, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaMemcpy(h_a, d_a, 2*N*sizeof(float), cudaMemcpyDeviceToHost));

    for (int i=0; i<2*N; i++) printf("%i %f\n", i, h_a[i]);

}
$ nvcc -arch=sm_20 -o t557 t557.cu
$ ./t557
0 0.145468
1 0.434899
2 0.820181
3 0.811845
4 0.550399
5 0.811845
6 0.294830
7 0.557235
8 0.914733
9 0.557235
10 0.868979
11 0.206681
12 0.321921
13 0.206681
14 0.782857
15 0.539587
16 0.011302
17 0.539587
18 0.285450
19 0.739071
$

现在，内核使每个线程按其顺序生成2个数字，但是在生成2个数字skipahead_sequence之间将使用备用线程。由于我选择1作为的第一个参数skipahead_sequence，因此函数调用的效果就像我要求输入1 * 2 ^ 67数字，即好像我已经调用了curand_uniform2 ^ 67次一样。这意味着线程对将生成的第二个数字相同，这正是我们在输出中的索引3和5、7和9、11和13等处看到的。

这是一些ASCII艺术：

主序列（长度2 ^ 190）（由seed参数确定）：

|0     ...   2^67-1 2^67 ...     2^68-1 2^68 ... ... ... 2^190|0 ... 
                                                              (main seq repeats)

序列（每个序列的长度为2 ^ 67）（由sequence参数确定）：

|seq0  ...   2^67-1|seq1   ...   2^68-1|seq2   ... 
                           ^
                   |offset |                     (determined by offset parameter)
                           |
                           RNG begins here for given seed, sequence(=seq1), offset