1D Min-convolution in CUDA

Question

I have two arrays, a and b, and I would like to compute the "min convolution" to produce result c. Simple pseudo code looks like the following:

for i = 0 to size(a)+size(b)
    c[i] = inf
    for j = 0 to size(a)
        if (i - j >= 0) and (i - j < size(b))
            c[i] = min(c[i], a[j] + b[i-j])

(edit: changed loops to start at 0 instead of 1)

If the min were instead a sum, we could use a Fast Fourier Transform (FFT), but in the min case, there is no such analog. Instead, I'd like to make this simple algorithm as fast as possible by using a GPU (CUDA). I'd be happy to find existing code that does this (or code that implements the sum case without FFTs, so that I could adapt it for my purposes), but my search so far hasn't turned up any good results. My use case will involve a's and b's that are of size between 1,000 and 100,000.

Questions:

• Does code to do this efficiently already exist?

• If I am going to implement this myself, structurally, how should the CUDA kernel look so as to maximize efficiency? I've tried a simple solution where each c[i] is computed by a separate thread, but this doesn't seem like the best way. Any tips in terms of how to set up thread block structure and memory access patterns?

Answer 1

I have used you algorithm. I think it'll help you.

const int Length=1000;

__global__ void OneD(float *Ad,float *Bd,float *Cd){
    int i=blockIdx.x;
    int j=threadIdx.x;
    Cd[i]=99999.99;
    for(int k=0;k<Length/500;k++){
        while(((i-j)>=0)&&(i-j<Length)&&Cd[i+k*Length]>Ad[j+k*Length]+Bd[i-j]){
            Cd[i+k*Length]=Ad[j+k*Length]+Bd[i-j];
    }}}

I have taken 500 Threads per block. And, 500 blocks per Grid. As, the number of threads per block in my device is restricted to 512, I used 500 threads. I have taken the size of all the arrays as Length (=1000).

Working:

1. i stores the Block Index and j stores the Thread Index.

2. The for loop is used as the number of threads are less than the size of the arrays.

3. The while loop is used for iterating Cd[n].

4. I have not used Shared Memory because, I have taken lots of blocks and threads. So, the amount of Shared Memory required for each block is low.

PS: If your device supports more Threads and Blocks, replace k<Length/500 with k<Length/(supported number of threads)

Answer 2

An alternative which might be useful for large a and b would be to use a block per output entry in c. Using a block allows for memory coalescing, which will be important in what is a memory bandwidth limited operation, and a fairly efficient shared memory reduction can be used to combine per thread partial results into a final per block result. Probably the best strategy is to launch as many blocks per MP as will run concurrently and have each block emit multiple output points. This eliminates some of the scheduling overheads associated with launching and retiring many blocks with relatively low total instruction counts.

An example of how this might be done:

#include <math.h>

template<int bsz>
__global__ __launch_bounds__(512)
void minconv(const float *a, int sizea, const float *b, int sizeb, float *c)
{
    __shared__ volatile float buff[bsz];
    for(int i = blockIdx.x; i<(sizea + sizeb); i+=(gridDim.x*blockDim.x)) {
        float cval = INFINITY;
        for(int j=threadIdx.x; j<sizea; j+= blockDim.x) {
            int t = i - j;
            if ((t>=0) && (t<sizeb))
                cval = min(cval, a[j] + b[t]);
        }
        buff[threadIdx.x] = cval; __syncthreads();
        if (bsz > 256) {
            if (threadIdx.x < 256) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+256]);
            __syncthreads();
        }
        if (bsz > 128) {
            if (threadIdx.x < 128) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+128]); 
            __syncthreads();
        }
        if (bsz > 64) {
            if (threadIdx.x < 64) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+64]);
            __syncthreads();
        }
        if (threadIdx.x < 32) {
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+32]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+16]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+8]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+4]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+2]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+1]);
            if (threadIdx.x == 0) c[i] = buff[0];
        }
    }
}

// Instances for all valid block sizes.
template __global__ void minconv<64>(const float *, int, const float *, int, float *);
template __global__ void minconv<128>(const float *, int, const float *, int, float *);
template __global__ void minconv<256>(const float *, int, const float *, int, float *);
template __global__ void minconv<512>(const float *, int, const float *, int, float *);

[disclaimer: not tested or benchmarked, use at own risk]

This is single precision floating point, but the same idea should work for double precision floating point. For integer, you would need to replace the C99 INFINITY macro with something like INT_MAX or LONG_MAX, but the principle remains the same otherwise.