Calculate squared Euclidean distance matrix on GPU

Question

Let p be a matrix of first set of locations where each row gives the coordinates of a particular point. Similarly, let q be a matrix of second set of locations where each row gives the coordinates of a particular point.

Then formula for pairwise squared Euclidean distance is:

k(i,j) = (p(i,:) - q(j,:))*(p(i,:) - q(j,:))',

where p(i,:) denotes i-th row of matrix p, and p' denotes the transpose of p.

I would like to compute matrix k on CUDA-enabled GPU (NVidia Tesla) in C++. I have OpenCV v.2.4.1 with GPU support but I'm open to other alternatives, like Thrust library. However, I'm not too familiar with GPU programming. Can you suggest an efficient way to accomplish this task? What C++ libraries should I use?

By the way, you're computing the square of the Euclidean distance... — Pedro, Jun 25 '12 at 21:49
@Pedro: I guess that is why it was described as "pairwise *squared* Euclidean distance" (emphasis mine)...... — talonmies, Jun 26 '12 at 10:00
@talonmies: Yup, by the looks of the timestamp of the edit, it just overlapped with my comment. — Pedro, Jun 26 '12 at 11:51

score 4 · Accepted Answer · edited Feb 02 '19 at 10:34

The problem looks simple enough to make a library overkill.

Without knowing the range of i and j, I'd suggest you partition k into blocks of a multiple of 32 threads each and in each block, compute

float sum, myp[d];
int i = blockIdx.x*blockDim.x + threadIdx.x;
for ( int kk = 0 ; kk < d ; kk++ )
    myp[kk] = p(i,kk);
for ( j = blockIdx.y*blockDim.y ; j < (blockIdx.y+1)*blockDim ; j++ ) {
    sum = 0.0f;
    #pragma unroll
    for ( int kk = 0 ; kk < d ; kk++ ) {
        temp = myp[kk] - q(j,kk);
        sum += temp*temp;
        }
    k(i,j) = sum;
    }

where I am assuming that your data has d dimensions and writing p(i,k), q(j,k) and k(i,j) to mean an access to a two-dimensional array. I also took the liberty in assuming that your data is of type float.

Note that depending on how k is stored, e.g. row-major or column-major, you may want to loop over i per thread instead to get coalesced writes to k.

thanks for the cuda c code! I have an OpenCV library (that I use anyway) that supports GPU matrix operations, so might try it first. — Alexey, Jun 25 '12 at 21:35

Calculate squared Euclidean distance matrix on GPU

1 Answers1

Linked