3. 设备API概述
要使用设备API，在定义使用mcRAND设备函数的内核文件时，需要包含文件mcrand_kernel.h。设备API包括用于 pseudorandom generation 和 quasirandom generation 函数。
3.1. 伪随机序列
伪随机序列的功能支持位（bit）生成和分布式生成。
3.1.1. 使用XORWOW和MRG32k3a的位生成
__device__ unsigned int
mcrand (mcrandState_t *state)
在调用mcrand_init()之后，mcrand()会返回一个周期大于 \(2^{190}\) 的伪随机数序列。如果每次调用mcrand()时使用相同的初始状态，并且在mcrand()多次调用之间不修改状态，那么生成的序列总是相同的。
__device__ void
mcrand_init (
unsigned long long seed, unsigned long long sequence,
unsigned long long offset, mcrandState_t *state)
mcrand_init()函数使用给定的种子、序列号和序列内偏移量来设置由调用者分配的初始状态。 不同的种子保证产生不同的起始状态和不同的序列。 相同的种子总是产生相同的状态和相同的序列。 设置的状态将是从种子状态开始经过 \(2^{67}\) * sequence + offset 次调用mcrand()后的状态。
虽然使用不同种子生成的序列通常不具有统计相关性，但某些特定的种子选择可能会产生例外。 使用相同种子和不同序列号生成的序列不会有统计相关性。
为了生成质量最高的并行伪随机数，每个实验都应该分配一个唯一的种子。在实验中，每个计算线程都应该分配一个唯一的序号。如果一次实验跨越了多次内核启动，那么建议内核启动之间的线程使用相同的种子，并以单调递增的方式分配序号。如果启动相同的线程配置，启动之间可以在全局内存中保存随机状态，以省去状态设置时间。
3.1.2. 使用MTGP32生成器的位生成
MTGP32 生成器改编自广岛大学开发的代码。
在该算法中，样本是为多个序列而生成，每个序列基于一组计算得到的参数。mcRAND使用了为32位生成器预先生成的200个参数集，周期为211214 。也可以生成其他参数集，并使用它们。每个参数集(序列)都有一个状态结构，该算法允许对200个序列中的每个序列进行多达256个并发线程(在单个块内)的线程安全生成和状态更新。
需要注意的是，两个不同的块不能对同一状态进行安全操作。还要注意的是，在一个块内，给定状态最多可以操作256个线程。
对于MTGP32生成器，提供了两个主机函数，用于帮助设置设备内存中不同序列的参数，以及设置初始状态。
__host__ mcrandStatust mcrandMakeMTGP32Constants(mtgp32paramsfastt params[],
mtgp32kernelparamst *p)
该函数将预生成格式(mtgp32_params_fast_t)的参数集数据重新组织为内核函数使用的格式(mtgp32_kernel_params_t)，并将其复制到设备内存中。
__host__ mcrandStatus_t
mcrandMakeMTGP32KernelState(mcrandStateMtgp32_t *s,
mtgp32_params_fast_t params[],
mtgp32_kernel_params_t *k,
int n,
unsigned long long seed)
该函数根据指定的参数集和种子初始化的n个状态，并将它们复制到由s指示的设备内存中。请注意，如果您使用预生成的状态，n的最大值为200。
mcRAND MTGP32生成器提供了两个内核函数来生成随机位。
__device__ unsigned int
mcrand (mcrandStateMtgp32_t *state)
该函数计算线程索引，并为该索引生成结果并更新状态。线程索引t的计算如下:
t= (blockDim.x * blockDim.y * threadIdx.z) + (blockDim.x * threadIdx.y) + threadIdx.x
该函数可以在单个内核启动中重复调用，但需要满足以下约束条件：
它只能在具有 256 个或更少线程的块中安全调用。
一个给定的状态不能被多个块使用。
一个给定的块可以使用多个状态生成随机数。
在代码的某个给定点，块中的所有线程要么全部调用此函数，要么全部不调用。
__device__ unsigned int
mcrandmtgp32specific(mcrandStateMtgp32_t *state, unsigned char index,
unsigned char n)
此函数根据线程特定的索引生成一个结果并更新状态，然后将状态中的偏移量向前推进n个位置。mcrand_mtgp32_specific可以在一个内核启动中多次调用，但需要满足以下约束条件：
对于给定的状态，最多256个线程可以调用这个函数。
在一个块内，对于给定的状态，如果有n个线程调用该函数，索引必须从0运行到n-1。索引不必与线程号匹配，并且可以根据调用程序的要求在线程之间分布。在代码的某个特定点，必须使用所有的索引从0到n-1，或者不使用任何索引。
一个给定的状态不能被多个块使用。
一个给定的块可以使用多个状态生成随机数。
图 3.1 MTGP32块和线程操作
图 3.1 是MTGP32中块和线程对生成器状态进行操作的示意图。每一行代表一个由 32 位整数 s(n) 组成的循环状态数组。操作数组的线程被标识为T(m)。所示的特定情况与主机API的内部实现相匹配，它启动了64个块，每个块有256个线程。每个块根据一组唯一的参数P(n)操作不同的序列。一个完整的MTGP32序列状态由351个32位整数定义。每个线程T(m)操作其中一个整数s(n+m)，将其与s(n+m+1)和一个取样元素s(n+m+p)（其中p <= 95）组合，它将新状态存储在状态数组的位置s(n+m+351)。线程同步后，基础索引n会根据已更新状态的线程数量进行推进。为了避免被覆盖，数组本身的长度必须至少为256 + 351个整数。实际上，为了索引的效率，它的大小为1024个整数。
对于可以操作给定状态数组的块中的线程数量的限制，是为了确保在需要作为取样状态之前，状态s(n+351)已经被更新。如果有一个线程T(256)，它可以使用s(n+256+95)，即在零号线程更新s(n+351)之前使用s(n+351)。如果一个应用程序需要在一个块中调用超过256个线程的MTGP32生成器函数，它必须使用多个MTGP32状态，可以通过使用多个参数集，或者使用具有不同种子的多个生成器来实现。还要注意，生成器函数在每次调用结束时同步线程，因此在一个块中调用256个线程来使用生成器是最有效的。
3.1.3. 使用Philox_4x32_10生成器的位生成
__device__ unsigned int
mcrand (mcrandState_t *state)
在调用mcrand_init()之后，mcrand()返回一个带句点的伪随机数序列 \(2^{128}\) 。如果每次以相同的初始状态调用mcrand()，并且在mcrand()多次调用之间没有修改状态，那么生成的序列总是相同的。
__device__ void
mcrand_init (
unsigned long long seed, unsigned long long subsequence,
unsigned long long offset, mcrandState_t *state)
mcrand_init()函数使用给定的种子、子序列和偏移量设置由调用者分配的初始状态。不同的种子保证产生不同的起始状态和不同的序列。子序列和偏移量共同定义了周期为 \(2^{128}\) 的序列中的偏移量。偏移量定义了长度为 \(2^{64}\) 的子序列中的偏移量。当子序列的最后一个元素被生成时，下一个随机数就是连续子序列的第一个元素。相同的种子总是产生相同的状态和序列。
虽然使用不同种子生成的序列通常不具有统计相关性，但某些特定的种子选择可能会产生例外。
为了获得最高质量的并行伪随机数生成，每个实验应分配一个唯一的种子值。在一个实验中，每个计算线程应分配一个唯一的ID号码。 如果一个实验跨越多个内核启动，建议在内核启动之间的线程给予相同的种子，并以单调递增的方式分配ID号码。如果启动相同配置的线程，可以在启动之间在全局内存中保留随机状态，以省去状态设置时间。
3.1.4. 分布
__device__ float
mcrand_uniform (mcrandState_t *state)
这个函数返回一个均匀分布在0.0到1.0之间的伪随机浮点数序列。它的返回范围是0.0到1.0，包含1.0，不包含0.0。分布函数可以使用基本生成器提供的任意数量的无符号整数值。使用的值的数量不保证是固定的。
__device__ float
mcrand_normal (mcrandState_t *state)
该函数返回一个具有均值为0.0和标准差为1.0的正态分布浮点数。此结果可以进行缩放和移动，以产生具有任意均值和标准差的正态分布值。
__device__ float
mcrand_log_normal (mcrandState_t *state, float mean, float stddev)
该函数基于具有给定均值和标准差的正态分布，返回一个单精度对数正态分布浮点数。
__device__ unsigned int
mcrand_poisson (mcrandState_t *state, double lambda)
该函数根据给定的lambda，返回一个泊松分布无符号整数。 用于从均匀分布的结果中得出泊松结果的算法取决于lambda的值和生成器的类型。有些算法在一个输出中会产生多个样本。
__device__ double
mcrand_uniform_double (mcrandState_t *state)
__device__ double
mcrand_normal_double (mcrandState_t *state)
__device__ double
mcrand_log_normal_double (mcrandState_t *state, double mean, double stddev)
上面的三个函数是双精度版本的mcrand_uniform()、mcrand_normal()和mcrand_log_normal()。
对于伪随机生成器，双精度函数使用多次调用mcrand()来生成53个随机位。
__device__ float2
mcrand_normal2 (mcrandState_t *state)
__device__ float2
mcrand_log_normal2 (mcrandState_t *state)
__device__ double2
mcrand_normal2_double (mcrandState_t *state)
__device__ double2
mcrand_log_normal2_double (mcrandState_t *state)
上述函数每次调用生成两个正态分布或对数正态分布的伪随机结果。由于底层实现使用了Box-Muller变换，这通常比每次调用生成一个结果更高效。
__device__ uint4
mcrand4 (mcrandStatePhilox4_32_10_t *state)
__device__ float4
mcrand_uniform4 (mcrandStatePhilox4_32_10_t *state)
__device__ float4
mcrand_normal4 (mcrandStatePhilox4_32_10_t *state)
__device__ float4
mcrand_log_normal4 (mcrandStatePhilox4_32_10_t *state, float mean, float stddev)
__device__ uint4
mcrand_poisson4 (mcrandStatePhilox4_32_10_t *state, double lambda)
__device__ uint4
mcrand_discrete4 (mcrandStatePhilox4_32_10_t *state, mcrandDiscreteDistribution_t discrete_distribution)
__device__ double2
mcrand_uniform2_double (mcrandStatePhilox4_32_10_t *state)
__device__ double2
mcrand_normal2_double (mcrandStatePhilox4_32_10_t *state)
__device__ double2
mcrand_log_normal2_double (mcrandStatePhilox4_32_10_t *state, double mean, double stddev)
上述函数每次调用生成四个单精度或两个双精度的结果。由于底层实现使用了Philox生成器，这通常比每次产生调用只生成一个结果更高效。
3.2. 准随机序列
尽管默认的生成器类型是来自XORWOW的伪随机数，但可以使用以下函数生成基于 Sobol’ 32位整数的 Sobol’序列：
__device__ void
mcrand_init (
unsigned int *direction_vectors,
unsigned int offset,
mcrandStateSobol32_t *state)
__device__ void
mcrand_init (
unsigned int *direction_vectors,
unsigned int scramble_c,
unsigned int offset,
mcrandStateScrambledSobol32_t *state)
__device__ unsigned int
mcrand (mcrandStateSobol32_t *state)
__device__ float
mcrand_uniform (mcrandStateSobol32_t *state)
__device__ float
mcrand_normal (mcrandStateSobol32_t *state)
__device__ float
mcrand_log_normal (
mcrandStateSobol32_t *state,
float mean,
float stddev)
__device__ unsigned int
mcrand_poisson (mcrandStateSobol32_t *state, double lambda)
__device__ double
mcrand_uniform_double (mcrandStateSobol32_t *state)
__device__ double
mcrand_normal_double (mcrandStateSobol32_t *state)
__device__ double
mcrand_log_normal_double (
mcrandStateSobol32_t *state,
double mean,
double stddev)
mcrand_init()函数用于初始化准随机数生成器的状态。它没有种子参数，只有方向向量和偏移量。对于Scrambled Sobol生成器，还有一个额外的参数scramble_c，它是Scrambled序列的初始值。对于mcrandStateSobol32_t类型和mcrandStateScrambledSobol32_t类型，方向向量是一个包含32个无符号整数值的数组。对于mcrandStateSobol64_t类型和mcrandStateScrambledSobol64_t类型，方向向量是一个包含64个无符号双长整数值的数组。对于32位Sobol生成器，偏移量和Scrambled序列的初始常数的类型是无符号整数。对于64位Sobol生成器，这些参数的类型是无符号双长整数。对于mcrandStateSobol32_t类型和mcrandStateScrambledSobol32_t类型，序列的长度恰好为 \(2^{32}\) 个元素，每个元素为32位。对于mcrandStateSobol64_t类型和mcrandStateScrambledSobol64_t类型，序列的长度恰好为 \(2^{64}\) 个元素，每个元素为64位。每次调用mcrand()都会返回下一个准随机元素。调用mcrand_uniform()返回从0.0到1.0的准随机浮点数或双精度数，其中1.0包含在内，0.0不包含在内。类似地，调用mcrand_normal()返回均值为0.0，标准差为1.0的正态分布浮点数或双精度数。调用mcrand_log_normal()返回从指定均值和标准差的正态分布派生的对数正态分布浮点数或双精度数。所有生成函数都可以使用任何类型的Sobol生成器进行调用。
例如，生成填充单位立方体的准随机坐标需要跟踪三个准随机生成器。 这三个生成器都从偏移量为0开始，维度分别为0、1和2。对于每个生成器状态，调用mcrand_uniform()一次将生成 \(x\), \(y\), 和 \(z\) 坐标。方向向量的表可以通过mcrandGetDirectionVectors32()和mcrandGetDirectionVectors64()函数在主机上访问。在使用之前，需要将所需的方向向量复制到设备内存中。
用于生成准随机序列的正态分布函数使用反累积密度函数来保持准随机序列的维数。因此，不存在像伪随机生成器那样一次生成多个结果的函数。
双精度Sobol32函数以双精度返回结果，使用底层生成器的32位内部精度。
双精度Sobol64函数以双精度返回结果，使用底层生成器的64位样本的高阶53位作为内部精度。
3.3. 跳过（Skip-Ahead）
有几个函数可以从生成器状态跳过。
__device__ void
skipahead(unsigned long long n, mcrandState_t *state)
__device__ void
skipahead(unsigned int n, mcrandStateSobol32_t *state)
使用这个函数等价于调用mcrand() \(n\) 次而不使用返回值，但速度更快。
__device__ void
skipahead_sequence(unsigned long long n, mcrandState_t *state)
这个函数相当于调用mcrand() \(n \cdot 2^{67}\) 次而不使用返回值，速度更快。
3.4. 用于离散分布的设备API
离散分布（如泊松分布）需要额外的API在主机端进行预处理，以生成特定分布的直方图。在泊松分布的情况下，不同的 lambda 值会产生不同的直方图。这些分布的最佳性能将在具有至少48KB L1缓存的GPU上体现。
mcrandStatus_t
mcrandCreatePoissonDistribution(
double lambda,
mcrandDiscreteDistribution_t *discrete_distribution)
mcrandCreatePoissonDistribution()函数用于创建给定lambda表达式下的泊松分布直方图。
__device__ unsigned int
mcrand_discrete (
mcrandState_t *state,
mcrandDiscreteDistribution_t discrete_distribution)
该函数根据给定离散分布直方图的分布返回单个离散无符号整型分布。
mcrandStatus_t
mcrandDestroyDistribution(
mcrandDiscreteDistribution_t discrete_distribution)
mcrandDestroyDistribution()函数用于清理与直方图相关的结构。
3.5. 性能说明
调用mcrand_init()比调用mcrand()或mcrand_uniform()慢。较大的偏移量对mcrand_init()的调用比较小的偏移量需要更多时间。保存和恢复随机生成器状态比重复计算起始状态要快得多。
如下所示，生成器状态可以在内核启动之间存储在全局内存中，使用本地内存进行快速生成，然后再存回全局内存中。
__global__ void example(mcrandState *global_state)
{
mcrandState local_state;
local_state = global_state[threadIdx.x];
for(int i = 0; i < 10000; i++) {
unsigned int x = mcrand(&local_state);
...
}
global_state[threadIdx.x] = local_state;
}
随机生成器状态的初始化通常需要比随机数生成更多的寄存器和本地内存。将mcrand_init()和mcrand()的调用分开成独立的内核可以获得最佳性能。
设置状态可能是一个昂贵的操作。加快设置的一种方法是为每个线程使用不同的种子和一个常数序列号为0。如果需要创建许多生成器，这种方法尤其有帮助。虽然设置速度更快，但该方法对生成序列的数学属性提供的保证较少。 如果从种子初始化生成器状态的哈希函数与生成器的周期性之间存在不良交互，那么对于某些种子值，可能会出现高度相关的输出线程。 我们还不知道有任何有问题的值；但即使存在，也很可能是罕见的。
3.6. 设备API示例(Device API Examples)
本示例使用mcRAND设备API，借助XORWOW或MRG32k3a生成器生成伪随机数。对于整数，它计算具有低位集的部分。对于均匀分布的实数，它计算出大于0.5的部分。对于正态分布的实数，它计算出在平均值的一个标准差范围内的部分
/*
* 这个程序使用设备MCRAND API来计算伪随机整数中具有低位设置的的比例。
* 然后生成统一结果，计算有多少大于0.5。
* 然后，它会生成正态结果，计算有多少在平均值的一个标准差内。
*/
#include <stdio.h>
#include <stdlib.h>
#include <mc_runtime.h>
#include <mcrand_kernel.h>
#define MC_CALL(x) do { if((x) != mcSuccess) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
__global__ void setup_kernel(mcrandState *state)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
/* 每个线程得到相同的种子，一个不同的序列号，
没有偏移量 */
mcrand_init(1234, id, 0, &state[id]);
}
__global__ void setup_kernel(mcrandStatePhilox4_32_10_t *state)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
/* 每个线程得到相同的种子，一个不同的序列号，
没有偏移量 */
mcrand_init(1234, id, 0, &state[id]);
}
__global__ void setup_kernel(mcrandStateMRG32k3a *state)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
/* 每个线程得到相同的种子，一个不同的序列号，
没有偏移量 */
mcrand_init(0, id, 0, &state[id]);
}
__global__ void generate_kernel(mcrandState *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
int count = 0;
unsigned int x;
/* 将状态复制到本地内存中以提高效率 */
mcrandState localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand(&localState);
/* Check if low bit set */
if(x & 1) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_kernel(mcrandStatePhilox4_32_10_t *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
int count = 0;
unsigned int x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStatePhilox4_32_10_t localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand(&localState);
/* Check if low bit set */
if(x & 1) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_uniform_kernel(mcrandState *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
float x;
/* 将状态复制到本地内存中以提高效率 */
mcrandState localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand_uniform(&localState);
/* Check if > .5 */
if(x > .5) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_uniform_kernel(mcrandStatePhilox4_32_10_t *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
float x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStatePhilox4_32_10_t localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand_uniform(&localState);
/* Check if > .5 */
if(x > .5) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_normal_kernel(mcrandState *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
float2 x;
/* 将状态复制到本地内存中以提高效率 */
mcrandState localState = state[id];
/* 生成伪随机的正态分布 */
for(int i = 0; i < n/2; i++) {
x = mcrand_normal2(&localState);
/* Check if within one standard deviaton */
if((x.x > -1.0) && (x.x < 1.0)) {
count++;
}
if((x.y > -1.0) && (x.y < 1.0)) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_normal_kernel(mcrandStatePhilox4_32_10_t *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
float2 x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStatePhilox4_32_10_t localState = state[id];
/* 生成伪随机正态分布 */
for(int i = 0; i < n/2; i++) {
x = mcrand_normal2(&localState);
/* Check if within one standard deviaton */
if((x.x > -1.0) && (x.x < 1.0)) {
count++;
}
if((x.y > -1.0) && (x.y < 1.0)) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_kernel(mcrandStateMRG32k3a *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
unsigned int x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStateMRG32k3a localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand(&localState);
/* Check if low bit set */
if(x & 1) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_uniform_kernel(mcrandStateMRG32k3a *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
double x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStateMRG32k3a localState = state[id];
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand_uniform_double(&localState);
/* Check if > .5 */
if(x > .5) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
__global__ void generate_normal_kernel(mcrandStateMRG32k3a *state,
int n,
unsigned int *result)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int count = 0;
double2 x;
/* 将状态复制到本地内存中以提高效率 */
mcrandStateMRG32k3a localState = state[id];
/* 生成伪随机的正态分布 */
for(int i = 0; i < n/2; i++) {
x = mcrand_normal2_double(&localState);
/* Check if within one standard deviaton */
if((x.x > -1.0) && (x.x < 1.0)) {
count++;
}
if((x.y > -1.0) && (x.y < 1.0)) {
count++;
}
}
/* 将状态复制回全局内存 */
state[id] = localState;
/* 存储结果 */
result[id] += count;
}
int main(int argc, char *argv[])
{
const unsigned int threadsPerBlock = 64;
const unsigned int blockCount = 64;
const unsigned int totalThreads = threadsPerBlock * blockCount;
unsigned int i;
unsigned int total;
mcrandState *devStates;
mcrandStateMRG32k3a *devMRGStates;
mcrandStatePhilox4_32_10_t *devPHILOXStates;
unsigned int *devResults, *hostResults;
bool useMRG = 0;
bool usePHILOX = 0;
int sampleCount = 10000;
bool doubleSupported = 0;
int device;
struct mcDeviceProp properties;
/* 检查是否支持双精度 */
MC_CALL(mcGetDevice(&device));
MC_CALL(mcGetDeviceProperties(&properties,device));
if ( properties.major >= 2 || (properties.major == 1 && properties.minor >= 3) ) {
doubleSupported = 1;
}
/* 检查MRG32k3a选项 (默认为XORWOW) */
if (argc >= 2) {
if (strcmp(argv[1],"-m") == 0) {
useMRG = 1;
if (!doubleSupported){
printf("MRG32k3a requires double precision\n");
printf("^^^^ test WAIVED due to lack of double precision\n");
return EXIT_SUCCESS;
}
}else if (strcmp(argv[1],"-p") == 0) {
usePHILOX = 1;
}
/* 允许覆盖样本计数 */
sscanf(argv[argc-1],"%d",&sampleCount);
}
/* 为主机上的结果分配空间 */
hostResults = (unsigned int *)calloc(totalThreads, sizeof(int));
/* 为设备上的结果分配空间 */
MC_CALL(mcMalloc((void **)&devResults, totalThreads *
sizeof(unsigned int)));
/* 将结果设置为0 */
MC_CALL(mcMemset(devResults, 0, totalThreads *
sizeof(unsigned int)));
/* 为设备上的prng状态分配空间 */
if (useMRG) {
MC_CALL(mcMalloc((void **)&devMRGStates, totalThreads *
sizeof(mcrandStateMRG32k3a)));
}else if(usePHILOX) {
MC_CALL(mcMalloc((void **)&devPHILOXStates, totalThreads *
sizeof(mcrandStatePhilox4_32_10_t)));
}else {
MC_CALL(mcMalloc((void **)&devStates, totalThreads *
sizeof(mcrandState)));
}
/* 设置prng状态 */
if (useMRG) {
setup_kernel<<<64, 64>>>(devMRGStates);
}else if(usePHILOX)
{
setup_kernel<<<64, 64>>>(devPHILOXStates);
}else {
setup_kernel<<<64, 64>>>(devStates);
}
/* 生成并使用伪随机参数 */
for(i = 0; i < 50; i++) {
if (useMRG) {
generate_kernel<<<64, 64>>>(devMRGStates, sampleCount, devResults);
}else if (usePHILOX){
generate_kernel<<<64, 64>>>(devPHILOXStates, sampleCount, devResults);
}else {
generate_kernel<<<64, 64>>>(devStates, sampleCount, devResults);
}
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy(hostResults, devResults, totalThreads *
sizeof(unsigned int), mcMemcpyDeviceToHost));
/* 显示结果 */
total = 0;
for(i = 0; i < totalThreads; i++) {
total += hostResults[i];
}
printf("Fraction with low bit set was %10.13f\n",
(float)total / (totalThreads * sampleCount * 50.0f));
/* 将结果设置为0*/
MC_CALL(mcMemset(devResults, 0, totalThreads *
sizeof(unsigned int)));
/* 生成和使用均匀的伪随机参数 */
for(i = 0; i < 50; i++) {
if (useMRG) {
generate_uniform_kernel<<<64, 64>>>(devMRGStates, sampleCount, devResults);
}else if(usePHILOX) {
generate_uniform_kernel<<<64, 64>>>(devPHILOXStates, sampleCount, devResults);
}else {
generate_uniform_kernel<<<64, 64>>>(devStates, sampleCount, devResults);
}
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy(hostResults, devResults, totalThreads *
sizeof(unsigned int), mcMemcpyDeviceToHost));
/* 显示结果 */
total = 0;
for(i = 0; i < totalThreads; i++) {
total += hostResults[i];
}
printf("Fraction of uniforms > 0.5 was %10.13f\n",
(float)total / (totalThreads * sampleCount * 50.0f));
/* 将结果设置为0*/
MC_CALL(mcMemset(devResults, 0, totalThreads *
sizeof(unsigned int)));
/* 生成和使用正态分布的伪随机参数 */
for(i = 0; i < 50; i++) {
if (useMRG) {
generate_normal_kernel<<<64, 64>>>(devMRGStates, sampleCount, devResults);
}else if(usePHILOX) {
generate_normal_kernel<<<64, 64>>>(devPHILOXStates, sampleCount, devResults);
}else {
generate_normal_kernel<<<64, 64>>>(devStates, sampleCount, devResults);
}
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy(hostResults, devResults, totalThreads *
sizeof(unsigned int), mcMemcpyDeviceToHost));
/* 显示结果 */
total = 0;
for(i = 0; i < totalThreads; i++) {
total += hostResults[i];
}
printf("Fraction of normals within 1 standard deviation was %10.13f\n",
(float)total / (totalThreads * sampleCount * 50.0f));
/* 清理 */
if (useMRG) {
MC_CALL(mcFree(devMRGStates));
}else if(usePHILOX)
{
MC_CALL(mcFree(devPHILOXStates));
}else {
MC_CALL(mcFree(devStates));
}
MC_CALL(mcFree(devResults));
free(hostResults);
printf("^^^^ kernel_example PASSED\n");
return EXIT_SUCCESS;
}
下面的示例使用mcRAND的主机MTGP的设置API和mcRAND的设备API，借助MTGP32生成器生成整数，并计算具有低位设置的比例。
/*
* 该程序使用设备MCRAND API来计算伪随机整数具有低位设置的比例。
*/
#include <stdio.h>
#include <stdlib.h>
#include <mc_runtime.h>
#include <mcrand_kernel.h>
/* 调用MTGP主机的helper函数 */
#include <mcrand_mtgp32_host.h>
/* 导入MTGP预先计算的参数集 */
#include <mcrand_mtgp32dc_p_11213.h>
#define MC_CALL(x) do { if((x) != mcSuccess) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
#define MCRAND_CALL(x) do { if((x) != MCRAND_STATUS_SUCCESS) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
__global__ void generate_kernel(mcrandStateMtgp32 *state,
int n,
int *result)
{
int id = threadIdx.x + blockIdx.x * 256;
int count = 0;
unsigned int x;
/* 生成伪随机的无符号整数 */
for(int i = 0; i < n; i++) {
x = mcrand(&state[blockIdx.x]);
/* Check if low bit set */
if(x & 1) {
count++;
}
}
/* 存储结果 */
result[id] += count;
}
int main(int argc, char *argv[])
{
int i;
long long total;
mcrandStateMtgp32 *devMTGPStates;
mtgp32_kernel_params *devKernelParams;
int *devResults, *hostResults;
int sampleCount = 10000;
/* 允许覆盖样本计数 */
if (argc == 2) {
sscanf(argv[1],"%d",&sampleCount);
}
/* 为主机上的结果分配空间 */
hostResults = (int *)calloc(64 * 256, sizeof(int));
/* 为设备上的结果分配空间 */
MC_CALL(mcMalloc((void **)&devResults, 64 * 256 *
sizeof(int)));
/* 将结果设置为0*/
MC_CALL(mcMemset(devResults, 0, 64 * 256 *
sizeof(int)));
/* 为设备上的prng状态分配空间 */
MC_CALL(mcMalloc((void **)&devMTGPStates, 64 *
sizeof(mcrandStateMtgp32)));
/* 设置MTGP的prng状态 */
/* 为MTGP内核参数分配空间 */
MC_CALL(mcMalloc((void**)&devKernelParams, sizeof(mtgp32_kernel_params)));
/* 将预定义的参数集重新格式化为内核格式, */
/* 并将内核参数复制到设备内存中 */
MCRAND_CALL(mcrandMakeMTGP32Constants(mtgp32dc_params_fast_11213, devKernelParams));
/* 为每个线程块初始化一个状态 */
MCRAND_CALL(mcrandMakeMTGP32KernelState(devMTGPStates,
mtgp32dc_params_fast_11213, devKernelParams, 64, 1234));
/* 状态设置已完成 */
/* 生成并使用伪随机参数 */
for(i = 0; i < 10; i++) {
generate_kernel<<<64, 256>>>(devMTGPStates, sampleCount, devResults);
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy(hostResults, devResults, 64 * 256 *
sizeof(int), mcMemcpyDeviceToHost));
/* 显示结果 */
total = 0;
for(i = 0; i < 64 * 256; i++) {
total += hostResults[i];
}
printf("Fraction with low bit set was %10.13g\n",
(double)total / (64.0f * 256.0f * sampleCount * 10.0f));
/* 清理 */
MC_CALL(mcFree(devKernelParams));
MC_CALL(mcFree(devMTGPStates));
MC_CALL(mcFree(devResults));
free(hostResults);
printf("^^^^ kernel_mtgp_example PASSED\n");
return EXIT_SUCCESS;
}
下面的示例使用 mcRAND 设备 API，使用64位Scrambled Sobol生成器生成均匀的双精度数。它利用生成的结果来推导出一个球体体积的近似值。
/*
*该程序使用设备 mcRAND API 来计算
* 准随机 3D 点落在球面内的比例
* 半径为1的球内，并且从中推导出球的体积。
*
* 本程序特地使用由
* mcrandGetDirectionVectors64 返回的64位
* Scrambled Sobol方向向量来生成双精度均匀样本。
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <mc_runtime.h>
#include <mcrand_kernel.h>
#include <mcrand.h>
#define THREADS_PER_BLOCK 64
#define BLOCK_COUNT 64
#define TOTAL_THREADS (THREADS_PER_BLOCK * BLOCK_COUNT)
/* 每个维度的64位向量数量 */
#define VECTOR_SIZE 64
#define MC_CALL(x) do { if((x) != mcSuccess) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
#define MCRAND_CALL(x) do { if((x) != MCRAND_STATUS_SUCCESS) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
/* 该内核为每个线程的x、y和z维度初始化状态 */
__global__ void setup_kernel(unsigned long long * sobolDirectionVectors,
unsigned long long *sobolScrambleConstants,
mcrandStateScrambledSobol64 *state)
{
int id = threadIdx.x + blockIdx.x * THREADS_PER_BLOCK;
int dim = 3*id;
/* 每个线程使用3个不同维度 */
mcrand_init(sobolDirectionVectors + VECTOR_SIZE*dim,
sobolScrambleConstants[dim],
1234,
&state[dim]);
mcrand_init(sobolDirectionVectors + VECTOR_SIZE*(dim + 1),
sobolScrambleConstants[dim + 1],
1234,
&state[dim + 1]);
mcrand_init(sobolDirectionVectors + VECTOR_SIZE*(dim + 2),
sobolScrambleConstants[dim + 2],
1234,
&state[dim + 2]);
}
/* 这个内核会生成随机的3D点，并且如果一个点在单位球内，
* 它会将计数器+1
*/
__global__ void generate_kernel(mcrandStateScrambledSobol64 *state,
int n,
long long int *result)
{
int id = threadIdx.x + blockIdx.x * THREADS_PER_BLOCK;
int baseDim = 3 * id;
long long int count = 0;
double x, y, z;
/* 生成准随机的双精度坐标 */
for(int i = 0; i < n; i++) {
x = mcrand_uniform_double(&state[baseDim]);
y = mcrand_uniform_double(&state[baseDim + 1]);
z = mcrand_uniform_double(&state[baseDim + 2]);
/* 检查是否在半径为1的球内 */
if( (x*x + y*y + z*z) < 1.0) {
count++;
}
}
/* 存储结果 */
result[id] += count;
}
int main(int argc, char *argv[])
{
int i;
long long total;
mcrandStateScrambledSobol64 *devSobol64States;
mcrandDirectionVectors64_t *hostVectors64;
unsigned long long int * hostScrambleConstants64;
unsigned long long int * devDirectionVectors64;
unsigned long long int * devScrambleConstants64;
long long int *devResults, *hostResults;
int sampleCount = 10000;
int iterations = 100;
double fraction;
double pi = 3.1415926535897932;
/* 允许覆盖样本数量 */
if (argc == 2) {
sscanf(argv[1],"%d",&sampleCount);
}
/* 为主机上的结果分配空间 */
hostResults = (long long int*)calloc(TOTAL_THREADS,
sizeof(long long int));
/* 为设备上的结果分配空间 */
MC_CALL(mcMalloc((void **)&devResults,
TOTAL_THREADS * sizeof(long long int)));
/* 把结果设置为0 */
MC_CALL(mcMemset(devResults, 0,
TOTAL_THREADS * sizeof(long long int)));
/* 获取指向64位Scrambled方向向量和常数的指针 */
MCRAND_CALL(mcrandGetDirectionVectors64( &hostVectors64,
MCRAND_SCRAMBLED_DIRECTION_VECTORS_64_JOEKUO6));
MCRAND_CALL(mcrandGetScrambleConstants64( &hostScrambleConstants64));
/* 为每个线程的3个状态(x,y,z)分配内存，每个状态用于获取特定的维度 */
MC_CALL(mcMalloc((void **)&devSobol64States,
TOTAL_THREADS * 3 * sizeof(mcrandStateScrambledSobol64)));
/* 分配内存并将每个线程的3组向量复制到设备 */
MC_CALL(mcMalloc((void **)&(devDirectionVectors64),
3 * TOTAL_THREADS * VECTOR_SIZE * sizeof(long long int)));
MC_CALL(mcMemcpy(devDirectionVectors64, hostVectors64,
3 * TOTAL_THREADS * VECTOR_SIZE * sizeof(long long int),
mcMemcpyHostToDevice));
/* 分配内存并将每个线程的3个Scramble常数（每个维度一个常数）
复制到设备 */
MC_CALL(mcMalloc((void **)&(devScrambleConstants64),
3 * TOTAL_THREADS * sizeof(long long int)));
MC_CALL(mcMemcpy(devScrambleConstants64, hostScrambleConstants64,
3 * TOTAL_THREADS * sizeof(long long int),
mcMemcpyHostToDevice));
/* 初始化状态 */
setup_kernel<<<BLOCK_COUNT, THREADS_PER_BLOCK>>>(devDirectionVectors64,
devScrambleConstants64,
devSobol64States);
/* 生成和计数准随机点 */
for(i = 0; i < iterations; i++) {
generate_kernel<<<BLOCK_COUNT, THREADS_PER_BLOCK>>>(devSobol64States, sampleCount, devResults);
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy(hostResults,
devResults,
TOTAL_THREADS * sizeof(long long int),
mcMemcpyDeviceToHost));
/* 统计并显示结果 */
total = 0;
for(i = 0; i < TOTAL_THREADS; i++) {
total += hostResults[i];
}
fraction = (double)total / ((double)TOTAL_THREADS * (double)sampleCount * (double)iterations);
printf("Fraction inside sphere was %g\n", fraction);
printf("(4/3) pi = %g, sampled volume = %g\n",(4.0*pi/3.0),8.0 * fraction);
/* 清理 */
MC_CALL(mcFree(devSobol64States));
MC_CALL(mcFree(devDirectionVectors64));
MC_CALL(mcFree(devScrambleConstants64));
MC_CALL(mcFree(devResults));
free(hostResults);
printf("^^^^ kernel_sobol_example PASSED\n");
return EXIT_SUCCESS;
}
3.7. Thrust 和 mcRAND 示例
下面的示例演示了如何混合使用 mcRAND 和 Thrust。 这是标准 Thrust 示例之一，monte_carlo.cpp 的最小修改版本。 该示例通过随机选择单位正方形中的点，并计算其到原点的距离，来判断其是否位于四分之一单位圆内，以估算 \(\pi\) 。
#include <thrust/iterator/counting_iterator.h>
#include <thrust/functional.h>
#include <thrust/transform_reduce.h>
#include <mcrand_kernel.h>
#include <iostream>
#include <iomanip>
// 我们可以调整 M 和 N 的值来找到性能最佳点
struct estimate_pi :
public thrust::unary_function<unsigned int, float>
{
__device__
float operator()(unsigned int thread_id)
{
float sum = 0;
unsigned int N = 10000; // 每个线程的样本数
unsigned int seed = thread_id;
mcrandState s;
// 给一个随机数生成器设定种子
mcrand_init(seed, 0, 0, &s);
// 在四分之一圆内取 N 个样本
for(unsigned int i = 0; i < N; ++i)
{
// 从单位正方形中获取样本
float x = mcrand_uniform(&s);
float y = mcrand_uniform(&s);
// 计算到原点的距离
float dist = sqrtf(x*x + y*y);
// 若 (u0,u1) 落在四分之一圆内，则将 1.0f 加入 sum
if(dist <= 1.0f)
sum += 1.0f;
}
// 乘以 4 得到整个圆的面积
sum *= 4.0f;
// 除以 N
return sum / N;
}
};
int main(void)
{
// 使用 30,000 个独立种子
int M = 30000;
float estimate = thrust::transform_reduce(
thrust::counting_iterator<int>(0),
thrust::counting_iterator<int>(M),
estimate_pi(),
0.0f,
thrust::plus<float>());
estimate /= M;
std::cout << std::setprecision(3);
std::cout << "pi is approximately ";
std::cout << estimate << std::endl;
return 0;
}
3.8. 泊松分布 API 示例
这个示例展示了泊松分布的三种API类型之间的区别。它是一个模拟商店排队的例子。主机 API (host API)对于生成泊松分布随机数的大向量是最稳健的（也就是它在lambda值的整个范围内具有最好的统计性质）。离散设备API(discrete Device API)几乎和主机API一样健壮，并允许在内核中生成泊松分布的随机数。简单设备API(simple Device API)是最不稳健的，但在为许多不同的lambda值生成泊松分布的随机数时更高效。
/*
* 这个程序使用MCRAND库的泊松分布
* 来模拟商店的队列，模拟时间为16小时。
* 它展示了使用3种不同的API的区别：
* - HOST API -顾客到达的频率由Poisson(4)来描述
* - SIMPLE DEVICE API -顾客到达的频率
* 由Poisson(4*(sin(x/100)+1))来描述，
* 其中x是从开店时间开始计算的分钟数。
* - ROBUST DEVICE API -顾客到达的频率如下：
* - 前3小时内由Poisson(2)来描述。
* - 接下来的3小时内由Poisson(1)来描述。
* - 6小时后由Poisson(3)来描述。
*/
#include <stdio.h>
#include <stdlib.h>
#include <mc_runtime.h>
#include <mcrand_kernel.h>
#include <mcrand.h>
#define MC_CALL(x) do { if((x) != mcSuccess) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)
#define MCRAND_CALL(x) do { if((x)!=MCRAND_STATUS_SUCCESS) { \
printf("Error at %s:%d\n",__FILE__,__LINE__);\
return EXIT_FAILURE;}} while(0)
#define HOURS 16
#define OPENING_HOUR 7
#define CLOSING_HOUR (OPENING_HOUR + HOURS)
#define access_2D(type, ptr, row, column, pitch)\
*((type*)((char*)ptr + (row) * pitch) + column)
enum API_TYPE {
HOST_API = 0,
SIMPLE_DEVICE_API = 1,
ROBUST_DEVICE_API = 2,
};
/* 全局变量 */
API_TYPE api;
int report_break;
int cashiers_load_h[HOURS];
__constant__ int cashiers_load[HOURS];
__global__ void setup_kernel(mcrandState *state)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
/* 每个线程使用相同的种子，不同的
序列号，无偏移 */
mcrand_init(1234, id, 0, &state[id]);
}
__inline__ __device__
void update_queue(int id, int min, unsigned int new_customers,
unsigned int &queue_length,
unsigned int *queue_lengths, size_t pitch)
{
int balance;
balance = new_customers - 2 * cashiers_load[(min-1)/60];
if (balance + (int)queue_length <= 0){
queue_length = 0;
}else{
queue_length += balance;
}
/* 存储结果 */
access_2D(unsigned int, queue_lengths, min-1, id, pitch)
= queue_length;
}
__global__ void simple_device_API_kernel(mcrandState *state,
unsigned int *queue_lengths, size_t pitch)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int new_customers;
unsigned int queue_length = 0;
/* 将状态复制到本地内存，以提高效率 */
mcrandState localState = state[id];
/* 在时间上模拟排队 */
for(int min = 1; min <= 60 * HOURS; min++) {
/* 根据API得出新顾客数量 */
new_customers = mcrand_poisson(&localState,
4*(sin((float)min/100.0)+1));
/* 更新队列 */
update_queue(id, min, new_customers, queue_length,
queue_lengths, pitch);
}
/* 将状态复制回全局内存 */
state[id] = localState;
}
__global__ void host_API_kernel(unsigned int *poisson_numbers,
unsigned int *queue_lengths, size_t pitch)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int new_customers;
unsigned int queue_length = 0;
/* 按时间模拟排队 */
for(int min = 1; min <= 60 * HOURS; min++) {
/* 从全局内存中获取随机数 */
new_customers = poisson_numbers
[blockDim.x * gridDim.x * (min -1) + id];
/* 更新队列 */
update_queue(id, min, new_customers, queue_length,
queue_lengths, pitch);
}
}
__global__ void robust_device_API_kernel(mcrandState *state,
mcrandDiscreteDistribution_t poisson_1,
mcrandDiscreteDistribution_t poisson_2,
mcrandDiscreteDistribution_t poisson_3,
unsigned int *queue_lengths, size_t pitch)
{
int id = threadIdx.x + blockIdx.x * 64;
unsigned int new_customers;
unsigned int queue_length = 0;
/* 将状态从复制到局部内存，以提高效率 */
mcrandState localState = state[id];
/* 按时间模拟排队 */
/* 前3小时 */
for(int min = 1; min <= 60 * 3; min++) {
/* 根据API得出新顾客数量 */
new_customers =
mcrand_discrete(&localState, poisson_2);
/* 更新队列 */
update_queue(id, min, new_customers, queue_length,
queue_lengths, pitch);
}
/* 接下来的3小时 */
for(int min = 60 * 3 + 1; min <= 60 * 6; min++) {
/* 根据API得出顾客数量 */
new_customers =
mcrand_discrete(&localState, poisson_1);
/* 更新队列 */
update_queue(id, min, new_customers, queue_length,
queue_lengths, pitch);
}
/* 6小时后 */
for(int min = 60 * 6 + 1; min <= 60 * HOURS; min++) {
/* 根据API得出新顾客数量 */
new_customers =
mcrand_discrete(&localState, poisson_3);
/* 更新队列 */
update_queue(id, min, new_customers, queue_length,
queue_lengths, pitch);
}
/* 将状态复制回全局内存 */
state[id] = localState;
}
/* 设置报告时间间隔 */
void report_settings()
{
do{
printf("Set time intervals between queue reports");
printf("(in minutes > 0)\n");
if (scanf("%d", &report_break) == 0) continue;
}while(report_break <= 0);
}
/* 设置每小时的收银员数量 */
void add_cachiers(int *cashiers_load)
{
int i, min, max, begin, end;
printf("Cashier serves 2 customers per minute...\n");
for (i = 0; i < HOURS; i++){
cashiers_load_h[i] = 0;
}
while (true){
printf("Adding cashier...\n");
min = OPENING_HOUR;
max = CLOSING_HOUR-1;
do{
printf("Set hour that cahier comes (%d-%d)",
min, max);
printf(" [type 0 to finish adding cashiers]\n");
if (scanf("%d", &begin) == 0) continue;
}while (begin > max || (begin < min && begin != 0));
if (begin == 0) break;
min = begin+1;
max = CLOSING_HOUR;
do{
printf("Set hour that cahier leaves (%d-%d)",
min, max);
printf(" [type 0 to finish adding cashiers]\n");
if (scanf("%d", &end) == 0) continue;
}while (end > max || (end < min && end != 0));
if (end == 0) break;
for (i = begin - OPENING_HOUR;
i < end - OPENING_HOUR; i++){
cashiers_load_h[i]++;
}
}
for (i = OPENING_HOUR; i < CLOSING_HOUR; i++){
printf("\n%2d:00 - %2d:00 %d cashier",
i, i+1, cashiers_load_h[i-OPENING_HOUR]);
if (cashiers_load[i-OPENING_HOUR] != 1) printf("s");
}
printf("\n");
}
/* 设置API类型 */
API_TYPE set_API_type()
{
printf("Choose API type:\n");
int choose;
do{
printf("type 1 for HOST API\n");
printf("type 2 for SIMPLE DEVICE API\n");
printf("type 3 for ROBUST DEVICE API\n");
if (scanf("%d", &choose) == 0) continue;
}while( choose < 1 || choose > 3);
switch(choose){
case 1: return HOST_API;
case 2: return SIMPLE_DEVICE_API;
case 3: return ROBUST_DEVICE_API;
default:
fprintf(stderr, "wrong API\n");
return HOST_API;
}
}
void settings()
{
add_cachiers(cashiers_load);
mcMemcpyToSymbol("cashiers_load", cashiers_load_h,
HOURS * sizeof(int), 0, mcMemcpyHostToDevice);
report_settings();
api = set_API_type();
}
void print_statistics(unsigned int *hostResults, size_t pitch)
{
int min, i, hour, minute;
unsigned int sum;
for(min = report_break; min <= 60 * HOURS;
min += report_break) {
sum = 0;
for(i = 0; i < 64 * 64; i++) {
sum += access_2D(unsigned int, hostResults,
min-1, i, pitch);
}
hour = OPENING_HOUR + min/60;
minute = min%60;
printf("%2d:%02d # of waiting customers = %10.4g |",
hour, minute, (float)sum/(64.0 * 64.0));
printf(" # of cashiers = %d | ",
cashiers_load_h[(min-1)/60]);
printf("# of new customers/min ~= ");
switch (api){
case HOST_API:
printf("%2.2f\n", 4.0);
break;
case SIMPLE_DEVICE_API:
printf("%2.2f\n",
4*(sin((float)min/100.0)+1));
break;
case ROBUST_DEVICE_API:
if (min <= 3 * 60){
printf("%2.2f\n", 2.0);
}else{
if (min <= 6 * 60){
printf("%2.2f\n", 1.0);
}else{
printf("%2.2f\n", 3.0);
}
}
break;
default:
fprintf(stderr, "Wrong API\n");
}
}
}
int main(int argc, char *argv[])
{
int n;
size_t pitch;
mcrandState *devStates;
unsigned int *devResults, *hostResults;
unsigned int *poisson_numbers_d;
mcrandDiscreteDistribution_t poisson_1, poisson_2;
mcrandDiscreteDistribution_t poisson_3;
mcrandGenerator_t gen;
/* 设置收银员、汇报频率(report)和API */
settings();
/* 为设备上的结果分配空间 */
MC_CALL(mcMallocPitch((void **)&devResults, &pitch,
64 * 64 * sizeof(unsigned int), 60 * HOURS));
/* 为主机上的结果分配空间 */
hostResults = (unsigned int *)calloc(pitch * 60 * HOURS,
sizeof(unsigned int));
/* 为设备上的prng状态分配空间 */
MC_CALL(mcMalloc((void **)&devStates, 64 * 64 *
sizeof(mcrandState)));
/* 设置prng状态 */
if (api != HOST_API){
setup_kernel<<<64, 64>>>(devStates);
}
/* 模拟排队 */
switch (api){
case HOST_API:
/* 创建伪随机数生成器 */
MCRAND_CALL(mcrandCreateGenerator(&gen,
MCRAND_RNG_PSEUDO_DEFAULT));
/* 设置种子 */
MCRAND_CALL(mcrandSetPseudoRandomGeneratorSeed(
gen, 1234ULL));
/* 计算n */
n = 64 * 64 * HOURS * 60;
/* 在设备上分配n个无符号整数 */
MC_CALL(mcMalloc((void **)&poisson_numbers_d,
n * sizeof(unsigned int)));
/* 在设备上生成n个无符号整数 */
MCRAND_CALL(mcrandGeneratePoisson(gen,
poisson_numbers_d, n, 4.0));
host_API_kernel<<<64, 64>>>(poisson_numbers_d,
devResults, pitch);
/* 清理 */
MCRAND_CALL(mcrandDestroyGenerator(gen));
break;
case SIMPLE_DEVICE_API:
simple_device_API_kernel<<<64, 64>>>(devStates,
devResults, pitch);
break;
case ROBUST_DEVICE_API:
/* 创建Poisson(1)的直方图 */
MCRAND_CALL(mcrandCreatePoissonDistribution(1.0,
&poisson_1));
/* 创建Poisson(2)的直方图 */
MCRAND_CALL(mcrandCreatePoissonDistribution(2.0,
&poisson_2));
/* 创建Poisson(3)的直方图 */
MCRAND_CALL(mcrandCreatePoissonDistribution(3.0,
&poisson_3));
robust_device_API_kernel<<<64, 64>>>(devStates,
poisson_1, poisson_2, poisson_3,
devResults, pitch);
/* 清理 */
MCRAND_CALL(mcrandDestroyDistribution(poisson_1));
MCRAND_CALL(mcrandDestroyDistribution(poisson_2));
MCRAND_CALL(mcrandDestroyDistribution(poisson_3));
break;
default:
fprintf(stderr, "Wrong API\n");
}
/* 将设备内存复制到主机 */
MC_CALL(mcMemcpy2D(hostResults, pitch, devResults,
pitch, 64 * 64 * sizeof(unsigned int),
60 * HOURS, mcMemcpyDeviceToHost));
/* 显示结果 */
print_statistics(hostResults, pitch);
/* 清理 */
MC_CALL(mcFree(devStates));
MC_CALL(mcFree(devResults));
free(hostResults);
return EXIT_SUCCESS;
}