我刚在库达,我正在尝试实现一个内核来计算我的大都会蒙特卡罗模拟的能量。
我将把这个函数的线性版本放在这里:
float calc_energy(struct frame frm, float L, float rc){
int i,j;
float E=0, rij, dx, dy, dz;
for(i=0; i<frm.natm; i++)
{
for(j=i+1; j<frm.natm; j++)
{
dx = fabs(frm.conf[j][0] - frm.conf[i][0]);
dy = fabs(frm.conf[j][1] - frm.conf[i][1]);
dz = fabs(frm.conf[j][2] - frm.conf[i][2]);
dx = dx - round(dx/L)*L;
dy = dy - round(dy/L)*L;
dz = dz - round(dz/L)*L;
/*rij*/
rij = sqrt(dx*dx + dy*dy + dz*dz);
if (rij <= rc)
{
E = E + (4*((1/pow(rij,12))-(1/pow(rij,6))));
}
}
}
return E;然后我尝试使用Cuda并行化:这是我的想法:
void calc_energy(frame* s, float L, float rc)
{
extern __shared__ float E;
int i = blockDim.x*blockIdx.x + threadIdx.x;
int j = blockDim.y*blockIdx.y + threadIdx.y;
float rij, dx, dy, dz;
dx = fabs(s->conf[j][0] - s->conf[i][0]);
dy = fabs(s->conf[j][1] - s->conf[i][1]);
dz = fabs(s->conf[j][2] - s->conf[i][2]);
dx = dx - round(dx/L)*L;
dy = dy - round(dy/L)*L;
dz = dz - round(dz/L)*L;
rij = sqrt(dx*dx + dy*dy + dz*dz);
if (rij <= rc)
{
E += (4*((1/pow(rij,12))-(1/pow(rij,6)))); //<- here is the big problem
}
} 我的主要问题是如何从每个线程求和变量E并将其返回给主机?我打算使用尽可能多的线程和块。
显然,在计算变量E时,代码的一部分丢失了。
我读过一些关于削减方法的内容,但我想知道这里是否有必要这样做。
我使用以下代码调用内核:
calc_energy<<<dimGrid,dimBlock>>>(d_state, 100, 5);编辑:
我知道我需要使用还原方法。幼崽对我做得很好。
继续代码的实现,我意识到我有一个新的问题,也许是因为我在这方面缺乏知识。
在我的嵌套循环中,变量(frm.natm)可以达到10^5的量级。考虑到我的GPU (GTX750 is ),每个块的线程数是1024,每个网格的块数是1024。如果我正确理解,内核中的最大运行数是1024x1024 = 1048576 (实际上要少)。
因此,如果我需要在嵌套循环中进行10^5 x 10^5 = 10^10的计算,那么考虑该算法的最佳方法是什么?选择一个固定的数字(适合我的GPU)并拆分计算将是一个好主意?
发布于 2018-05-29 07:27:20
我的主要问题是如何从每个线程求和变量E并将其返回给主机?
您需要首先使用某种形式的分块并行缩减来在块级别上对每个线程进行求和(为此我建议使用CUB块wise缩减实现 )。
一旦每个块都有来自其线程的部分和,则需要合并块和。这可以通过每个块中的一个线程在原子上完成,也可以通过第二个内核调用(带一个块),或者在主机上完成。如何和在哪里使用最后的和将确定哪些选项是最适合您的应用程序。
发布于 2019-01-22 05:59:07
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
#include <algorithm>
#include <cstdlib>
int main(void)
{
thrust::host_vector<int> h_vec(100);
std::generate(h_vec.begin(), h_vec.end(), rand);
thrust::device_vector<int> d_vec = h_vec;
int x = thrust::reduce(d_vec.begin(), d_vec.end(), 0, thrust::plus<int>());
std::cout<< x<< std::endl;
return 0;
}https://stackoverflow.com/questions/50574268
复制相似问题