先看一个C递归函数:
uint32_t c_recursion(uint32_t t, uint32_t depth) {
uint32_t res = 0;
if (depth == 0) { //底层返回
res = t;
}else {
res += c_recursion(hash(t * 10086u), depth - 1); //向下传递本层处理过的参数
res += c_recursion(hash(t * 5201314u), depth - 1);
}
return hash(res * 0x47a92fu); //本层返回 合并的、经过本层处理的值
}
给定一个初始值,从这个初始值递推分叉出去,最后逐层返回到最上层,在C递归中很简单的函数
在CUDA中实现的思路:
递推过程中,本层参数与下层第一分叉参数,共享同一个tid,第二分叉tid加本层线程总数
完整代码如下:
using namespace std;
//#include <windows.h>
struct Timer {
LARGE_INTEGER tc; //最低频率
LARGE_INTEGER start;//起始时间
LARGE_INTEGER t1;//上次时间
void init() {
QueryPerformanceFrequency(&tc);//系统精确计时的最小单位值
QueryPerformanceCounter(&start); //时间计数开始
t1 = start;
}
//与上次的计时
double TimingButton_ms() {
LARGE_INTEGER now;
QueryPerformanceCounter(&now);
double res = (double)((now.QuadPart - t1.QuadPart) * 1000.0 / tc.QuadPart);
t1 = now;
return res;
}
//与最初的计时
double TotalTime_ms() {
LARGE_INTEGER now;
QueryPerformanceCounter(&now);
return (double)((now.QuadPart - start.QuadPart) * 1000.0 / tc.QuadPart);
}
};
// hash作者:https://nullprogram.com/blog/2018/07/31/
//bias: 0.17353355999581582 ( very probably the best of its kind )
__device__ __host__
uint32_t lowbias32(uint32_t x) {
x ^= x >> 16;
x *= 0x7feb352dU;
x ^= x >> 15;
x *= 0x846ca68bU;
x ^= x >> 16;
return x;
}
#define hash(x) lowbias32((x))
uint32_t c_recursion(uint32_t t, uint32_t depth) {
uint32_t res = 0;
if (depth == 0) { //底层返回
res = t;
}else {
res += c_recursion(hash(t * 10086u), depth - 1); //向下传递本层处理过的参数
res += c_recursion(hash(t * 5201314u), depth - 1);
}
return hash(res * 0x47a92fu); //本层返回 合并的、经过本层处理的值
}
__global__
void first_push(uint32_t* now_params) {
const uint32_t tid = (blockIdx.y * gridDim.x * blockDim.x) + (blockIdx.x * blockDim.x) + threadIdx.x;
now_params[tid] = 5;
}
__global__
void push(uint32_t* now_params, uint32_t* next_params) {
const uint32_t tid = (blockIdx.y * gridDim.x * blockDim.x) + (blockIdx.x * blockDim.x) + threadIdx.x;
const uint32_t now_total = gridDim.y * gridDim.x * blockDim.x;
uint32_t t = now_params[tid];
next_params[tid] = hash(t * 10086u);
next_params[tid + now_total] = hash(t * 5201314u);
}
__global__
void first_back(uint32_t* now_retuenvalue, uint32_t* now_params) {
const uint32_t tid = (blockIdx.y * gridDim.x * blockDim.x) + (blockIdx.x * blockDim.x) + threadIdx.x;
now_retuenvalue[tid] = hash(now_params[tid]* 0x47a92fu);
}
__global__
void back(uint32_t* now_retuenvalue, uint32_t* now_params, uint32_t* last_retuenvalue) {
const uint32_t tid = (blockIdx.y * gridDim.x * blockDim.x) + (blockIdx.x * blockDim.x) + threadIdx.x;
const uint32_t now_total = gridDim.y * gridDim.x * blockDim.x;
uint32_t res = last_retuenvalue[tid]+ last_retuenvalue[tid + now_total];
now_retuenvalue[tid] = hash(res * 0x47a92fu);
}
std::map<std::string, void*> devMemory;
std::map<std::string, void*> hostMemory;
void cleanup() {
for (auto iter = devMemory.begin(); iter != devMemory.end(); ++iter) {
cudaFree(iter->second);
cout << "cleanup devMemory:" << iter->first << endl;
}
for (auto iter = hostMemory.begin(); iter != hostMemory.end(); ++iter) {
cudaFree(iter->second);
cout << "cleanup hostMemory:" << iter->first << endl;
}
}
#define ERRORINSP(status,str) if (status != cudaSuccess) {cout << str <<"\n" <<cudaGetErrorString(status) << endl; cleanup(); return 0;}
int main(void) {
const uint32_t max_depth = 20;//30层使用内存8GB 可能超过显存
const uint32_t value_count = (1u << max_depth) - 1;
Timer timer={};
timer.init();
uint32_t f = c_recursion(5, max_depth-1);
cout << "C递归计时:" << timer.TotalTime_ms() << "ms" << endl;
std::printf("C递归结果: 0x%X \n", f);
timer.init();
cudaError_t cudaStatus;
cudaStatus = cudaMalloc((void**)&devMemory["params"], value_count * sizeof(uint32_t));
ERRORINSP(cudaStatus, "cudaMalloc failed!");
cudaStatus = cudaMalloc((void**)&devMemory["retValues"], value_count * sizeof(uint32_t));
ERRORINSP(cudaStatus, "cudaMalloc failed!");
size_t memsize = value_count * sizeof(uint32_t) * 2;
cout << "已创建的设备内存大小:" << memsize/(1ULL<<20) <<"MB" << endl;
first_push << < 1, 1 >> > ((uint32_t*)devMemory["params"]); //初始值
cudaStatus = cudaDeviceSynchronize();
ERRORINSP(cudaStatus, "launching Kernel:")
cudaStatus = cudaGetLastError();
ERRORINSP(cudaStatus, "Kernel launch failed:")
uint32_t p_offset = 0;
for (int i = 0; i < max_depth - 1; i++) {
uint32_t n = 0x1U << i;
uint32_t* now_params = (uint32_t*)devMemory["params"] + p_offset;
uint32_t* next_params = (uint32_t*)devMemory["params"] + p_offset + n;
dim3 bolck = { 1,1,1 };
dim3 thread = { n,1,1 };
if (i > 10) {
thread.x = 1u << 10u;
bolck.x = 1u << (i-10);
}
push << < bolck, thread >> > (now_params, next_params); //参数递推
cudaStatus = cudaDeviceSynchronize();
ERRORINSP(cudaStatus,"launching Kernel:")
cudaStatus = cudaGetLastError();
ERRORINSP(cudaStatus, "Kernel launch failed:")
p_offset += n;
}
p_offset = value_count;
for (int i = max_depth-1; i >= 0; i--) {
uint32_t n = 0x1U << i;
p_offset -= n;
uint32_t* now_params = (uint32_t*)devMemory["params"] + p_offset;
uint32_t* now_retuenvalue = (uint32_t*)devMemory["retValues"] + p_offset;
dim3 bolck = { 1,1,1 };
dim3 thread = { n,1,1 };
if (i > 10) {
thread.x = 1u << 10u;
bolck.x = 1u << (i - 10);
}
if (i == max_depth - 1) {
first_back << <bolck, thread >> > (now_retuenvalue, now_params); //底层返回
}else{
uint32_t* last_retuenvalue = (uint32_t*)devMemory["retValues"] + p_offset+ n;
back << <bolck, thread >> > (now_retuenvalue, now_params, last_retuenvalue); //返回值
}
cudaStatus = cudaDeviceSynchronize();
ERRORINSP(cudaStatus, "launching Kernel:")
cudaStatus = cudaGetLastError();
ERRORINSP(cudaStatus, "Kernel launch failed:")
}
cout << "CUDA计时:" << timer.TotalTime_ms() << "ms" << endl;
uint32_t data = 0;
cudaStatus = cudaMemcpy(&data, devMemory["retValues"], 1 * sizeof(uint32_t), cudaMemcpyDeviceToHost);
std::printf("CUDA递归结果: 0x%X \n", data);
// system("pause");
cleanup();
return 0;
}
运行结果如下
C递归计时:8.3675ms
C递归结果: 0x96AA2994
已创建的设备内存大小:7MB
CUDA计时:301.847ms
CUDA递归结果: 0x96AA2994
cleanup devMemory:params
cleanup devMemory:retValues
结果是正确的,说明效果达到了。
如果提高一个维度,把初始的0看做一个像素
我使用它来改造用C语言画光(七):比尔-朗伯定律 - 知乎 (zhihu.com)
1024*1024像素,255次采样,4层追踪,:3分钟 ,5层:8分钟,6层:13分钟
场景更复杂的情况下比原版快十倍
cuda还可以用2048*2048模拟4次采样,这些是同时发生的