Asked by Hao Zhang
on 13 Dec 2018

Hi, I am developing a SPH (smoothed particle hydrodynamic) code of solving fluid equations in Matlab. I have successfully vectorized and implemented the code on GPU. The speed up is astonishing, ~10 times faster than CPU code. In the process I got invaluable suggestions from Matt J and Joss Knight for speeding up the code. Thank them for their wonderful suggestions.

In SPH, one has to find the neighbor particles in a given radius for every particles in the domain. I use the matlab function knnsearch for this purpose. Now the other part of the code (except neighbor search but solving the fluid equations) is so fast that the limiting part now is the knnsearch, which uses kdtree algorithm runing on CPU. It takes 85% percent of the runtime (see the following code profiler results)

The function 'knnCPU_kdtree_func' uses the matlab built-in function knnsearch with kdtree algorithm runing on CPU. The other functions are solving the real fluid equations runing on GPU only consumes 10% of the total time.

I wonder is there any GPU implementation of k-nearest neighbor search that I can free download and using as a function call in my matlab code? Many thanks.

Answer by Matt J
on 14 Dec 2018

Edited by Matt J
on 14 Dec 2018

Accepted Answer

A quick Google search turned up this,

They have Matlab bindings in their legacy version.

Hao Zhang
on 14 Dec 2018

Many thanks! This seems to be a very relevant code!

I have downloaded and tried to use nvcc to compile their knncuda.cu file into a ptx file. This step succeed.

In the next step I use the matlab GPUcoder to create CUDAKernel object and trying to run it on GPU, for the entry I choose '_Z17compute_distancesPfiiS_iiiS_'.

clear;close all;clc;

gd=gpuDevice();

reset(gd);

N=int32(10);

rx=rand(N,1);

ry=rand(N,1);

ref=[gpuArray(single(rx)) gpuArray(single(ry))];

ref_width=N;

ref_pitch=int32(2);

query=ref;

query_width=N;

query_pitch=int32(2);

height=int32(2);

dist=zeros(N,1,'single','gpuArray');

%%% create CUDAKernel object %%%

knncuda_kernel=parallel.gpu.CUDAKernel('knncuda.ptx','knncuda.cu','_Z17compute_distancesPfiiS_iiiS_');

%%% create CUDAKernel object %%%

%%% run CUDAKernel %%%

result=feval(knncuda_kernel,ref,ref_width,ref_pitch,query,query_width,query_pitch,height,dist)

%%% run kernel %%%

No error is reported, however I always get the same function return (result) as one of the input argument (ref). Why is that? I have very limited knowledge of cuda C... Could you please help? Thanks!

Matt J
on 14 Dec 2018

It's all 3rd party black boxes to me, I'm afraid.

Sign in to comment.

Answer by Matt J
on 13 Dec 2018

Edited by Matt J
on 13 Dec 2018

In SPH, one has to find the neighbor particles in a given radius for every particles in the domain. I use the matlab function knnsearch for this purpose.

Perhaps instead, you should use rangesearch with a pre-constructed KDTreeSearcher object. You could then use parfor to process different chunks of query data in parallel.

Hao Zhang
on 13 Dec 2018

Hi Matt, I found using knnsearch followed by finding which neighbors are within radius is faster than directly call rangesearch. Just have to use an enough large K so that all neighbors within radius is included in the K-nearest neighbors.

I tried to use parfor to process chunks of query data in parallel, but the execution time is even longer than call knnsearch once incuding all query points. The following is a test code:

clear;clc;close all

n_parallel=4;

if isempty(gcp('nocreate'))==1

pooljob=parpool('local',n_parallel);

end

N=1e5;

knn_K=50;

rx=rand(N,1);

ry=rand(N,1);

BucketSize_kdtree=50;

Mdl=KDTreeSearcher([rx ry],'Distance','euclidean','BucketSize',BucketSize_kdtree); %%% kdtree search object

idx_Neighbor_cell=cell(n_parallel,1);

d_Neighbor_cell=cell(n_parallel,1);

logic_beyor_cell=cell(n_parallel,1);

N_par=round(N/n_parallel);

tic;

[idx_Neighbor1,d_Neighbor1]=knnsearch(Mdl,[rx ry],'K',knn_K);

toc

tic;

rx_query_cell=mat2cell(rx,[N_par*ones(1,n_parallel-1) N-N_par*(n_parallel-1)],1);

ry_query_cell=mat2cell(ry,[N_par*ones(1,n_parallel-1) N-N_par*(n_parallel-1)],1);

parfor i=1:n_parallel

[idx_Neighbor_cell{i},d_Neighbor_cell{i}]=knnsearch(Mdl,[rx_query_cell{i} ry_query_cell{i}],'K',knn_K);

end

idx_Neighbor=cell2mat(idx_Neighbor_cell);

d_Neighbor=cell2mat(d_Neighbor_cell);

toc;

The runtime is:

Elapsed time is 0.389322 seconds.

Elapsed time is 1.199775 seconds.

The parfor version takes much longer time despite using 4 cores.

Does anyone know why this is and how to optimize it further? Thanks!

Sign in to comment.

Answer by Joss Knight
on 13 Dec 2018

knnsearch is supported on the GPU, so just use it!

Hao Zhang
on 14 Dec 2018

Thanks for the explanations. I tried pdist2 with gputimeit and got nearly the same results. The argument on allocation of large matrix is a good point. It is indeed not clear how much the time is expended on distance calculations. But if I run this

clear;clc;close all

N=2.5e4;

knn_K=50;

rx=rand(N,1);

ry=rand(N,1);

BucketSize_kdtree=50;

tic;

Mdl=KDTreeSearcher([rx ry],'Distance','euclidean','BucketSize',BucketSize_kdtree); %%% kdtree search object

[idx_Neighbor1,d_Neighbor1]=knnsearch(Mdl,[rx ry],'K',knn_K);

toc

gd=gpuDevice();

reset(gd);

rx=gpuArray(single(rx));

ry=gpuArray(single(ry));

for i=1:10

tic;

D=pdist2([rx ry],[rx ry],'euclidean');

wait(gd);

toc;

end

Elapsed time is 0.097022 seconds.

Elapsed time is 0.175050 seconds.

Elapsed time is 0.199612 seconds.

Elapsed time is 0.187805 seconds.

Elapsed time is 0.191508 seconds.

Elapsed time is 0.190413 seconds.

Elapsed time is 0.191164 seconds.

Elapsed time is 0.190247 seconds.

Elapsed time is 0.191383 seconds.

Elapsed time is 0.180551 seconds.

Elapsed time is 0.181040 seconds.

Does it mean that not much time is spent on allocation since D is already there? Or each time when pdist2 is called, the memory is recallocated again? That sounds a bit stupid to me, is there any way to let matlab just allocate the memory once? Since I have to run neighbor finding at every time step of my simulation, it is ok that only the first time takes longer. Thanks

Matt J
on 14 Dec 2018

Hao Zhang
on 14 Dec 2018

Ok, I am willing to write a cuda C code, but I am not a computer scientist, I would imagine cuda C will be a steep learning curve for me...

I wonder is there any standard cuda libaray that I can free download to do this?

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.