Is one CUDA block dimension faster than the other?

Question

I have a simple CUDA code that assigns the values of an NxN matrix A to matrix B. In one case, I declare block sizes block(1,32) and have each thread loop over the entries in the first matrix dimension. In the second case, I declare block sizes block(32,1) and have each thread loop over entries in the second matrix dimension.

Is there some really obvious reason why, in my code below, threads that loop over the stride 1 memory are significantly slower than those that the loop over stride N memory? I would have thought it was the other way around (if there is any difference at all).

Am I missing something really obvious (a bug, perhaps)?

The complete code is below.

#include <stdio.h>
#include <sys/time.h>

__global__ void addmat_x(int m, int n, int* A, int *B) 
{    
    int idx, ix;
    int iy = threadIdx.y + blockIdx.y*blockDim.y;
    if (iy < n) 
        for(ix = 0; ix < m; ix++) {
            idx  = iy*m + ix;    /* iy*m is constant */
            B[idx]   = A[idx];
        }
}

__global__ void addmat_y(int m, int n, int* A, int *B) 
{    
    int ix = threadIdx.x + blockIdx.x*blockDim.x;
    int idx, iy;
    if (ix < m)
        for(iy = 0; iy < n; iy++) {
            idx  = iy*m + ix; 
            B[idx]   = A[idx];        
        }
}

double cpuSecond()
{
    struct timeval tp;
    gettimeofday(&tp,NULL);
    return (double) tp.tv_sec + (double)tp.tv_usec*1e-6;
}

int main(int argc, char** argv) 
{
    int *A, *B;
    int *dev_A, *dev_B;
    size_t m, n, nbytes;
    double etime, start;

    m = 1 << 14;  
    n = 1 << 14;  
    nbytes = m*n*sizeof(int);

    A = (int*) malloc(nbytes);
    B = (int*) malloc(nbytes);

    memset(A,0,nbytes);

    cudaMalloc((void**) &dev_A, nbytes);
    cudaMalloc((void**) &dev_B, nbytes);
    cudaMemcpy(dev_A, A, nbytes, cudaMemcpyHostToDevice);

#if 1
    /* One thread per row */
    dim3 block(1,32);  
    dim3 grid(1,(n+block.y-1)/block.y);
    start = cpuSecond();
    addmat_x<<<grid,block>>>(m,n,dev_A, dev_B);
#else
    /* One thread per column */
    dim3 block(32,1);  
    dim3 grid((m+block.x-1)/block.x,1);
    start = cpuSecond();
    addmat_y<<<grid,block>>>(m,n,dev_A, dev_B);
#endif
    cudaDeviceSynchronize();
    etime = cpuSecond() - start;
    printf("GPU Kernel %10.3g (s)\n",etime);

    cudaFree(dev_A);
    cudaFree(dev_B);
    free(A);
    free(B);

    cudaDeviceReset();
}

32 threads grouped together in x ie. (32,1) will coalesce on loads and stores. This is a more efficient use of the GPU memory subsystem. — Robert Crovella
@Robert Crovella Thanks. In other words, the thread grouping is more important than the layout of global memory access? — Donna
They are related. In order to do an efficient global load or store, you want adjacent threads in a warp to read adjacent locations in memory, roughly speaking. This is a very basic GPU programming concept covered in many questions and many tutorials. In the "good" case, a warp read is behaving something like what is depicted on slide 45 here. In the bad case, a warp read is behaving something like what is depicted on slide 53. — Robert Crovella
@Robert Crovella The slides are great. At the time the warp is launched, how can it be known what global memory addresses will be needed? A second question is, why does the dimension of 1d threadblocks matter, if they are all treated as 1d arrays by the SM? — Donna
At the time the warp is launched, how can it be known what global memory addresses will be needed? You don't know this at the time the warp is launched. You have to inspect the code. You can look at a global memory read operation, and compute the indices that each thread in the warp will generate, just by looking at your code. why does the dimension of 1d threadblocks matter, if they are all treated as 1d arrays by the SM? It doesn't matter. What matters is the relationship between the threads and the indexing created by your code. Study the indices your code generates, in each case. — Robert Crovella

Robert Crovella Robert Crovella · Accepted Answer · 2018-04-14T01:45:55

Lets compare the global memory indexing generated by each thread, in each case.

addmat_x:

Your block dimension is (1,32). This means 1 thread wide in x, 32 threads "long" in y. The threadId.x value for each thread will be 0. The threadIdx.y value for the threads in the warp will range from 0 to 31, as you move from thread to thread in the warp. With that, let's inspect your creation of idx in that kernel:

m = 1 << 14;  
...
int iy = threadIdx.y + blockIdx.y*blockDim.y;
idx  = iy*m + ix;

let's choose the first block, whose blockIdx.y is 0. Then:

idx  = threadIdx.y*(1<<14) + ix;

For the first loop iteration, ix is 0. The idx values generated by each thread will be:

threadIdx.y:   | idx:
   0              0
   1                (1<<14)
   2              2*(1<<14)
  ...
   31            31*(1<<14)

For a given loop iteration, the distance from the load or store index from one thread to the next will be 1<<14. i.e. not adjacent. Scattered.

addmat_y:

Your block dimension is (32,1). This means 32 threads wide in x, 1 thread "long" in y. The threadIdx.y value for each thread will be 0. The threadIdx.x value for the threads in the warp will range from 0 to 31, as you move from thread to thread. Now let's inspect your creation of idx in that kernel:

m = 1 << 14;  
...
int ix = threadIdx.x + blockIdx.x*blockDim.x;
idx  = iy*m + ix;

Let's choose the first block, whose blockIdx.x is 0. Then:

idx  = iy*m + threadIdx.x;

For the first loop iteration, iy is 0, so we simply have:

idx  = threadIdx.x;

This generates the following index pattern across the warp:

threadIdx.x:   | idx:
   0              0
   1              1
   2              2
  ...
   31            31

These indices are adjacent, it is not a scattered load or store, the addresses will coalesce nicely, and this represents "efficient" use of global memory. It will perform faster than the first case.

Is one CUDA block dimension faster than the other?

1 Answers