简而言之,在您的问题中,您提到使用网格跨步循环作为折叠所需状态的方法。我认为这种方法或类似的方法(如@talonmies建议的)是最明智的方法。选择一种方法,将线程数减少到保持机器繁忙/充分利用所需的最低限度。然后让每个线程计算多个操作,重新使用提供的随机生成器状态。
我从你的代码 shell 开始,把它变成了一个经典的“hello world”MC 问题:根据正方形面积与内切圆面积之比计算 pi,随机生成点来估计面积。
那么我们将考虑3种方法:
创建一个大型一维网格以及每个线程的状态,以便每个线程计算一个随机点并对其进行测试。
创建一个更小的一维网格以及每个线程的状态,并允许每个线程计算多个随机点/测试,以便生成与情况 1 相同数量的点/测试。
创建一个与方法 2 大小相同的网格,同时创建一个“唯一驻留线程索引”,然后提供足够的状态来覆盖唯一驻留线程。每个线程将有效地使用通过“常驻线程索引”提供的状态,而不是普通的全局唯一线程索引。计算与情况 1 相同数量的点/测试。
“唯一驻留线程索引”的计算并不是一件小事。在主机代码中:
我们必须确定理论上可以驻留的最大块数。我为此使用了一个简单的启发式方法,但可以说还有更好的方法。我只需将每个多处理器的最大驻留线程数除以每个块所选线程的数量。为此我们必须使用整数除法。
然后初始化足够的状态以覆盖最大块数乘以 GPU 上的 SM 数。
在设备代码中:
- 确定我在哪个 SM 上。
- 使用该 SM,执行原子测试以检索每个 SM 的唯一块号。
- 将 1 和 2 的结果相结合,再加上
threadIdx.x
变量来创建“唯一驻留线程索引”。这将成为我们典型的线程索引,而不是通常的计算。
从结果的角度来看,可以得出以下结论:
与其他答案中所述相反,初始化时间是not微不足道。它不应该被忽视。对于这个简单的问题,初始化内核运行时间超过了计算内核运行时间。因此,我们最重要的结论是,我们应该寻求方法来最小化随机生成器状态的创建。我们当然不想不必要地重新运行初始化。因此,对于这个特定的代码/测试,我们可能应该根据这些结果放弃方法 1。
从计算内核运行时的角度来看,一个内核没有什么值得称赞的。
从代码复杂度的角度来看,方法 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