Naive Thread-Per-Element CUDA GEMM

The Naive Strategy: One Thread Per Output Element

The standard matrix multiplication computes $C=AB$, where $A$ is an $M\times K$ matrix, $B$ is a $K\times N$ matrix, and the resulting matrix $C$ is $M\times N$. The value of each element $C_{ij}$ is the dot product of the $i$-th row of $A$ and the $j$-th column of $B$:

Since the calculation for each $C_{ij}$ is independent of all other elements in $C$, this problem is a natural fit for parallelization. The most straightforward strategy is to assign exactly one thread to compute each element of the output matrix $C$.

To do this, we will launch a 2D grid of threads that has the same dimensions as our output matrix $C$. Each thread will use its unique 2D index to identify which row of $A$ and which column of $B$ it needs to process.

Mapping Thread Indices to Matrix Coordinates

As you learned in Module 4, CUDA provides built-in variables—threadIdx, blockIdx, blockDim, and gridDim—that a thread can use to determine its identity and location within the grid. For our 2D grid, a thread can compute its global row and col index in the output matrix $C$ like this:

// Get the global row and column index for this thread
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

Here, blockIdx.y * blockDim.y calculates the starting row index of the current thread block, and threadIdx.y provides the offset within that block. The same logic applies to the column index using the .x dimension.

From 2D Coordinates to 1D Memory

While we think of matrices in 2D, they are stored in the GPU's global memory as a flat, 1D array. The standard memory layout in C/C++ is row-major order, where all elements of the first row are stored contiguously, followed by the second row, and so on.

To access an element at (row, col) in a matrix with W columns, we use the formula: index = row * W + col. Applying this to our GEMM kernel:

• To access $A_{ik}$ (element at (row, i) in matrix $A$), the index is row * K + i.
• To access $B_{kj}$ (element at (i, col) in matrix $B$), the index is i * N + col.

With this indexing math, we have all the components needed to write our kernel.

Implementing the Naive Kernel

The following video provides a complete, from-scratch walkthrough of creating a naive matrix multiplication kernel. It covers the host-side setup (allocating memory, setting grid/block dimensions) and the device-side kernel implementation.

After watching the video segments, the resulting kernel should look conceptually like this:

__global__ void sgemm_naive(int M, int N, int K, const float* A, const float* B, float* C) {
    // Calculate the global row and column that this thread is responsible for.
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    // Boundary check to ensure the thread is within the bounds of the output matrix C.
    // This is necessary when M or N are not perfect multiples of the block dimensions.
    if (row < M && col < N) {
        float value = 0.0f;
        // Compute the dot product of the i-th row of A and the j-th column of B.
        for (int i = 0; i < K; ++i) {
            // A has K columns, B has N columns.
            value += A[row * K + i] * B[i * N + col];
        }
        // Write the final result to the corresponding element in C.
        C[row * N + col] = value;
    }
}

This kernel correctly implements matrix multiplication. If you were to run it and verify the results against a standard CPU implementation, it would produce the right answer. However, its performance is far from optimal.

The Villain: Why This Kernel is Slow

Our naive kernel is what we might call a "villain"—it's correct on the surface, but it hides a fatal performance flaw. The issue lies not in the computation, but in the memory access pattern.

Let's analyze the work done by threads within a single block. Consider the threads responsible for the first row of a block (i.e., threadIdx.y = 0). All of these threads (with threadIdx.x from 0 to 31, for example) are working on the same output row. To compute their respective dot products, they all need to read the exact same row from matrix $A$.

The following resource explains this problem with an excellent analogy and clear diagrams.

As the video explains, each thread makes its own independent trip to the slow, off-chip global memory (the "pantry") to fetch data. When multiple threads in a block need the same data, they all perform this slow walk independently, resulting in a catastrophic waste of memory bandwidth. For a $1024\times 1024$ matrix, each element of $A$ is re-read from global memory 1024 times!

This kernel is overwhelmingly memory-bound. The powerful CUDA cores spend most of their time idle, waiting for data to arrive. This is the central problem we must solve to achieve high performance.

Conclusion

In this lesson, you have successfully implemented your first CUDA kernel for matrix multiplication. You've translated the mathematical definition of GEMM into a parallel algorithm by assigning one thread to each element of the output matrix, and you've mastered the indexing required to map 2D matrix coordinates to 1D physical memory.

Key Takeaways:

• The naive GEMM strategy maps a 2D grid of threads directly to the output matrix $C$, with each thread computing one dot product.
• Global thread indices (row, col) are computed using a combination of blockIdx, blockDim, and threadIdx.
• Accessing elements in a row-major 1D array requires converting (row, col) coordinates to a linear index: row * width + col.
• While functionally correct, the naive kernel is extremely inefficient due to massive redundant global memory access, making it severely memory-bound.

You have now built and diagnosed your performance "villain." In our next lesson, we will introduce the hero that will defeat it: shared memory. You will learn to implement a tiled matrix multiplication algorithm, a fundamental technique that caches data in fast on-chip memory to dramatically reduce global memory traffic and begin our journey toward a truly high-performance kernel.