Hi Guys,
I’m writing a matmul kernel, with different approaches to measure the performance.
Input and kernel config:
Both input matrix A and B are 2048 x 2048 matrices. I have a (32, 32, 0) as grid dimension and (8,8,0) as block dimension.
My kernel setup is as follows:
Each block is responsible for computing 64 x 64 sub matrix of the resultant matrix. Since we have (8 x 8) block, each thread will compute 8 x 8 sub matrix of the resultant matrix.
Kernel also uses shared memory for two 64 x 64 matrices for both matrix A and matrix B.
For instance, since tile size is 8 each thread will load (8 x 2048) sub matrix from the matrix A and (2048 x 8) sub matrix from the matrix B to the shared memory and perform matmul to obtain (8 x 8) sub matrix . Due to the SM size (48 KB) limitation per thread block. I’m splitting the SM loading step into multiple time steps, where in each time step a single thread will load strides of size 64.
__global__ void matmul_2048(const int *A, const int *B, int nrows,
int ncols, int tile_sz, int *result) {
int global_tidx = blockIdx.x * blockDim.x + threadIdx.x;
int global_tidy = blockIdx.y * blockDim.y + threadIdx.y;
int thr_mat_range;
int startRowIdx2D = global_tidx * tile_sz;
int startColIdx2D = global_tidy * tile_sz;
const int chunk_size = 64;
__shared__ int matA[chunk_size][chunk_size], matB[chunk_size][chunk_size];
int chunks = nrows / chunk_size;
/*
0. Split rows into chunk_size chunks (aka Timesteps)
1. Load 64 x 64 submatrix of A into SM matA
2. Load 64 x 64 submatrix of B into SM matB
3. Perform 64 x 64 matmul operation and write to global memory.
*/
for (int i=0; i < chunks; i++) {
//corresponding start row index of 2D matrix A, the
//thread is responsible for.
startRowIdx2D = global_tidx * tile_sz;
for (int t=0; t < tile_sz; t++) {
int rowIdx1D = startRowIdx2D * ncols + i * chunk_size;
thr_mat_range = threadIdx.x * tile_sz + t;
for (int j=0; j < chunk_size; j++) {
matA[thr_mat_range][j] = A[rowIdx1D + j];
}
startRowIdx2D++;
}
//corresponding start column index of 2D matrix B, the
//thread is responsible for.
startColIdx2D = global_tidy * tile_sz;
int startRow = i * chunk_size;
for (int t=0; t < tile_sz; t++) {
thr_mat_range = threadIdx.y * tile_sz + t;
for (int j=0; j < chunk_size; j++) {
int colIdx1D = (startRow + j) * ncols + startColIdx2D;
matB[j][thr_mat_range] = B[colIdx1D];
}
startColIdx2D++;
}
//why we need syncthreads for correctness??
__syncthreads();
//start target row of the output 2D matrix
int outrow = global_tidx * tile_sz;
//perform matmul of 64x64 chunks from both matrix A and B.
for (int rowt=0; rowt<tile_sz; rowt++) {
int row = threadIdx.x * tile_sz + rowt;
//start target column of the output 2D matrix.
int outcol = global_tidy * tile_sz;
for (int colt=0; colt<tile_sz; colt++) {
int col = threadIdx.y * tile_sz + colt;
int temp = 0;
for (int k=0; k < chunk_size; k++) {
temp += matA[row][k] * matB[k][col];
}
result[outrow * ncols + outcol] += temp;
outcol++;
}
outrow++;
}
}
return;
}
Each thread in the entire gird access independent addresses without any overlapping between any two threads. I’m not sure why I need to perform __syncthreads() after loading the sub matrices into shared memory before performing the matmul operation. I think since each thread’s address access is different, I don’t think there is a datarace. I have also run the kernel under compute-sanitizer racecheck tool, but the tool didn’t give any errors.
Interestingly, If I replace the temp with a direct write to the result array, I get more mismatches between expected (CPU) and produced (GPU) values. Delaying the write with an additional register (temp), produces fewer mismatched values.
Could you please help me with the reasoning for the need of __sync_threads for correctness?
Please let me know if I’m making any mistakes here.
Thank you very much,