求达人指教一个高效的CUDA并行算法!

[hejian29]

假设我有n个点做为样本， 存放在global memory中。
我要求另外m个点每个点到n个点中的点的最短距离。 m个点也是存放在global memory中
也就是说，对于m个点中的每个点。 我要求出这个点到每个样本点（n)的距离， 然后取最小值。

我现在是这么做的，但效率非常低。

望达人指教一个高效的算法。


//surfaceParticlePos是n个样本点
//int index = __mul24(blockIdx.x,blockDim.x) + threadIdx.x; 是当前第m个被求解的点
//每个线程处理一个被求解的点。

__global__ void dcomputeParticle2SuraceDist(float4* particlePos, float4* surfaceParticlePos, float* dist, uint iNumSurfParticles)
{
int index = __mul24(blockIdx.x,blockDim.x) + threadIdx.x;
float3 pos = make_float3(particlePos[index]);
float _dist = 999.0f;

for(int i = 0; i < iNumSurfParticles; i++)
{
float3 spos = make_float3(surfaceParticlePos[i]);
float l = length(spos - pos);
if( l < _dist )
_dist = l;
}
dist[index] = dist;
}

[springleo]

1.n和m的取值范围大概是多少？这个很重要，

（1）如果n很大，那么可以分成n个线程，m个点可以分别作为函数的参数传入，

__global__ void calc_dist_kernel(float4 pos,float4* surface_pos,float *dist)
{

}
求最小值可以用并行，复杂度从O(n) 变成O(log(n))，
__global__ void get_min_dist_kernel(float *dist,float *min_dist)
{

}


for(i=0;i<m;i++)
{
calc_dist_kernel....
get_min_dist_kernel...
}




2.如果n很小，m很大，那就分成m个线程，把n 个点拷到共享内存中


总之
访问全局内存一定要对齐，如果不能对齐，改变数据结构和算法直到对齐

尽可能利用共享内存和寄存器，以减少全局内存访问次数

[Cyrosly]

虽然有些限制,但你先试试

#define __GPU_ARCH__ 10
#define __GPU_arch_sm_10__ 10
#define __GPU_arch_sm_11__ 11
#define __GPU_arch_sm_12__ 12
#define __GPU_arch_sm_13__ 13
#define architecture(s) __GPU_arch_sm_##s##__

#if __GPU_ARCH__<archtecture(13)
#define IMUL(x,y) __umul24(x,y)
#else
#define IMUL(x,y) ((x)*(y))
#endif

#if __GPU_ARCH__<architecture(13)
#define MAXSAMPLE 512
#else
#define MAXSAMPLE 768
#endif

__device__ float distance(const float4& v0,const float4& v1)
{
    float ds=v0.x-v1.x;
    ds*=ds;
    float t=v0.y-v1.y;
    ds=t*t+ds;
    t=v0.z-v0.z;
    ds=t*t+ds;
    return sqrtf(ds);
}

__global__ void findNearestRefPoints(float* dis,float4* particle,float4* sample)
{
    uint tidx=IMUL(blockDim.x,blockIdx.x)+threadIdx.x;
    uint tidy=IMUL(blockDim.y,blockIdx.y)+threadIdx.y;
    uint gloc=IMUL(IMUL(gridDim.x,blockDim.x),tidy)+tidx;

    __shared__ float smem[MAXSAMPLE];
    float4 src=particle[gloc];
    uint sloc=IMUL(blockDim.x,threadIdx.y)+threadIdx.x;
    float4 ref=sample[sloc];
    smem[sloc]=distance(src,ref);
    
    uint offset=IMUL(blockDim.x,blockDim.y); offset>>=1; //必须是2的幂

    while(offset>0){
        __syncthreads();
        if(sloc<offset){
            if(smem[sloc]>smem[sloc+offset]){
                smem[sloc]=smem[sloc+offset]; //不需要保存原值
            }
        }
        offset>>=1;
    }

    if(!(threadIdx.x&threadIdx.y)){
        gloc=IMUL(gridDim.x,blockIdx.y)+blockIdx.x;
        dis[gloc]=smem[0];
    }
}
    

[Cyrosly]

忘记说了,我现在是0积分

[Cyrosly]

如果sample的数量不大,可以考虑存储在const cache里

[OpenHero]

呵呵，这个需要基础扎实，不然出来的东西还是有问题~，非最优化的算法实现这个已经有人做了；要做必须得作出更好的结果；

需要有几个地方注意，内存的访问，和线程的调度

[Cyrosly]

不行,是每个点都要与所有的参考点一一比较.所以上边的程序不可以,所以改成这样:

#define __GPU_ARCH__ 10
#define __GPU_arch_sm_10__ 10
#define __GPU_arch_sm_11__ 11
#define __GPU_arch_sm_12__ 12
#define __GPU_arch_sm_13__ 13
#define architecture(s) __GPU_arch_sm_##s##__

#if __GPU_ARCH__<archtecture(13)
#define IMUL(x,y) __umul24(x,y)
#else
#define IMUL(x,y) ((x)*(y))
#endif

#if __GPU_ARCH__<architecture(13)
#define MAXSAMPLE 768
#else
#define MAXSAMPLE 1024
#endif

#define I0 0
#define I1 !I0
#define SWAP() (I0=I1)

__global__ void findNearestRefPoints(float* dis,float4* particle,float4* sample)
{
    uint tidx=IMUL(blockDim.x,blockIdx.x)+threadIdx.x;
    uint tidy=IMUL(blockDim.y,blockIdx.y)+threadIdx.y;
    uint gloc=IMUL(IMUL(gridDim.x,blockDim.x),tidy)+tidx;

    __shared__ float smem[MAXSAMPLE];
    
    uint sloc=IMUL(blockDim.x,threadIdx.y)+threadIdx.x;
    float4 ref=sample[sloc];
    float dst[2]={.0f,.0f};
    
    uint barea=IMUL(blockDim.x,blockDim.y);
    uint offset=area>>1;
    gloc=IMUL(gridDim.x,blockIdx.y)+blockIdx.x;
    
    for(uint i=0;i<MAXSAMPLE;i++)
    {    
        float4 src=particle[gloc*barea+i];
        smem[sloc]=distance(src,ref);
        while(offset>0){
            __syncthreads();
            if(sloc<offset){
                if(smem[sloc]>smem[sloc+offset]){
                    smem[sloc]=smem[sloc+offset];
                }
            }
            offset>>=1;
        }
        dst[I0]=smem[0];
        SWAP();
        if(dst[I0]<dst[I1]){
            dst[I0]=dst[I1];
        }
    }

    if(!(threadIdx.x&threadIdx.y)){
        dis[gloc]=dst[I0];
    }
}