Pycuda Blocks and Grids to work with big datas

I need help to know the size of my blocks and grids. I’m building a python app to perform metric calculations based on scipy as: Euclidean distance, Manhattan, Pearson, Cosine, joined other.

The project is PycudaDistances.

It seems to work very well with small arrays. When I perform a more exhaustive test, unfortunately it did not work. I downloaded movielens set (http://www.grouplens.org/node/73).

Using Movielens 100k, I declared an array with shape (943, 1682). That is, users are 943 and 1682 films evaluated. The films not by a classifier user I configured the value to 0.

With a much larger array algorithm no longer works. I face the following error:

pycuda._driver.LogicError: cuFuncSetBlockShape failed: invalid value.

Researching this error, I found an explanation of telling Andrew that supports 512 threads to join and to work with larger blocks it is necessary to work with blocks and grids.

I wanted a help to adapt the algorithm Euclidean distance arrays to work from small to giant arrays.

def euclidean_distances(X, Y=None, inverse=True):
    X, Y = check_pairwise_arrays(X,Y)
    rows = X.shape[0]
    cols = Y.shape[0]
    solution = numpy.zeros((rows, cols))
    solution = solution.astype(numpy.float32)

    kernel_code_template = """
    #include <math.h>
    
    __global__ void euclidean(float *x, float *y, float *solution) {

        int idx = threadIdx.x + blockDim.x * blockIdx.x;
        int idy = threadIdx.y + blockDim.y * blockIdx.y;
        
        float result = 0.0;
        
        for(int iter = 0; iter < %(NDIM)s; iter++) {
            
            float x_e = x[%(NDIM)s * idy + iter];
            float y_e = y[%(NDIM)s * idx + iter];
            result += pow((x_e - y_e), 2);
        }
        int pos = idx + %(NCOLS)s * idy;
        solution[pos] = sqrt(result);
    }
    """
    kernel_code = kernel_code_template % {
        'NCOLS': cols,
        'NDIM': X.shape[1]
    }

    mod = SourceModule(kernel_code)

    func = mod.get_function("euclidean")
    func(drv.In(X), drv.In(Y), drv.Out(solution), block=(cols, rows, 1))

    return numpy.divide(1.0, (1.0 + solution)) if inverse else solution

For more details see: https://github.com/vinigracindo/pycudaDistances/blob/master/distances.py

Answer

To size the execution paramters for your kernel you need to do two things (in this order):

1. Determine the block size

Your block size will mostly be determined by hardware limitations and performance. I recommend reading this answer for more detailed information, but the very short summary is that your GPU has a limit on the total number of threads per block it can run, and it has a finite register file, shared and local memory size. The block dimensions you select must fall inside these limits, otherwise the kernel will not run. The block size can also effect the performance of kernel, and you will find a block size which gives optimal performance. Block size should always be a round multiple of the warp size, which is 32 on all CUDA compatible hardware released to date.

2. Determine the grid size

For the sort of kernel you have shown, the number of blocks you need is directly related to the amount of input data and the dimensions of each block.

If, for example, your input array size was 943×1682, and you had a 16×16 block size, you would need a 59 x 106 grid, which would yield 944×1696 threads in the kernel launch. In this case the input data size is not a round multiple of the block size, you will need to modify your kernel to ensure that it doesn’t read out-of-bounds. One approach could be something like:

__global__ void euclidean(float *x, float *y, float *solution) {
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    int idy = threadIdx.y + blockDim.y * blockIdx.y;

     if ( ( idx < %(NCOLS)s ) && ( idy < %(NDIM)s ) ) {

        .....
     }
}

The python code to launch the kernel could look like something similar to:

bdim = (16, 16, 1)
dx, mx = divmod(cols, bdim[0])
dy, my = divmod(rows, bdim[1])

gdim = ( (dx + (mx>0)) * bdim[0], (dy + (my>0)) * bdim[1]) )
func(drv.In(X), drv.In(Y), drv.Out(solution), block=bdim, grid=gdim)

This question and answer may also help understand how this process works.

Please note that all of the above code was written in the browser and has never been tested. Use it at your own risk.

Also note it was based on a very brief reading of your code and might not be correct because you have not really described anything about how the code is called in your question.