Naive Matrix Multiplication

Welcome back to 2D Grids and Matrix Math. We are now in lesson 2, which means we already have an important tool from the course: mapping a 2D thread layout to matrix rows and columns. In the last lesson, we used that idea to place values into their correct positions. In this lesson, we will use the same mapping for something much more useful: naive matrix multiplication. By the end, we will have a CUDA kernel where each thread computes one element of an output matrix.

Before we look at the code, let's build a mental model. If matrix A has shape rowsA x colsA and matrix B has shape colsA x colsB, then the result matrix C has shape rowsA x colsB.

Each output entry is a dot product:

This fits the GPU very naturally:

1. A single thread chooses a single (row, col) in C;
2. That thread loops over k;
3. It stores the final sum into one output position.

So, the big idea is simple: one thread per output cell.

We create a small test case on the CPU to understand the math before scaling up. Matrix $A$ is $2 \times 3$ and matrix $B$ is $3 \times 2$, which results in a $2 \times 2$ matrix $C$:

In our code, we represent these matrices as flat 1D arrays using row-major order. We also prepare an h_expected array so the program can verify the kernel's output:

int main() {
    const int rowsA = 2, colsA = 3, colsB = 2;
    const size_t bytesA = rowsA * colsA * sizeof(float);
    const size_t bytesB = colsA * colsB * sizeof(float);
    const size_t bytesC = rowsA * colsB * sizeof(float);

    // Row-major representation of the matrices above
    float h_A[] = {1.0f, 2.0f, 3.0f,
                   4.0f, 5.0f, 6.0f};
    float h_B[] = {7.0f,  8.0f,
                   9.0f, 10.0f,
                   11.0f, 12.0f};
    
    float h_C[rowsA * colsB] = {};
    float h_expected[] = {58.0f, 64.0f, 139.0f, 154.0f};

Notice that we only define three variables for the dimensions: rowsA, colsA, and colsB. Because matrix multiplication is only valid when the number of columns in $A$ matches the number of rows in $B$, a separate rowsB variable would be redundant. Using three variables (often referred to as $M, K$, and $N$ in mathematical libraries) reinforces this requirement and keeps the logic cleaner.

The first expected value is easy to check by hand. For the top-left result ($C_{0,0}$), we compute the dot product of the first row of $A$ and the first column of $B$:

$1 \cdot 7 + 2 \cdot 9 + 3 \cdot 11 = 58$

Using a small, hand-verifiable test case is a best practice when developing GPU kernels, as it makes debugging indexing errors much simpler.

Now, we define the kernel. As you may recall from the previous lesson, 2D indexing provides each thread with a global row and col. The key difference now is what those coordinates mean. They are no longer just labels for memory positions; they tell the thread which element of the output matrix C it owns.

/**
 * Naive matrix multiplication kernel.
 * A is rowsA x colsA
 * B is colsA x colsB
 * C is rowsA x colsB
 */
__global__ void matrixMultiply(const float* A, const float* B, float* C,
                               int rowsA, int colsA, int colsB) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

A few design choices are worth noting. The input pointers are const because the kernel reads from A and B but does not change them. Also, the y-direction maps to rows, and the x-direction maps to columns, which matches the pattern we built in lesson 1.

Once a thread knows its output coordinate, it can compute the dot product for that specific cell. A bounds check is still necessary because the grid may be rounded up, launching a few extra threads. Inside the loop, k walks across one row of A and one column of B. This version is called naive because every thread reads directly from global memory, with no shared memory optimization yet.

    if (row < rowsA && col < colsB) {
        float sum = 0.0f;
        for (int k = 0; k < colsA; k++) {
            sum += A[row * colsA + k] * B[k * colsB + col];
        }
        C[row * colsB + col] = sum;
    }
}

The three index formulas are the heart of the kernel:

1. A[row * colsA + k]: element from row row of A;
2. B[k * colsB + col]: element from column col of B;
3. C[row * colsB + col]: output location for the finished sum.

This is exactly where 2D matrix thinking meets the 1D memory layout we covered in the previous lesson.

With the kernel ready, we can prepare the device-side data. We allocate three buffers on the GPU: one for d_A, one for d_B, and one for d_C. Then, we copy only the input matrices from the host to the device. The output matrix does not require an initial copy because the kernel will write all its values.

    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, bytesA));
    CUDA_CHECK(cudaMalloc(&d_B, bytesB));
    CUDA_CHECK(cudaMalloc(&d_C, bytesC));

    CUDA_CHECK(cudaMemcpy(d_A, h_A, bytesA, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, bytesB, cudaMemcpyHostToDevice));

The launch configuration follows the same 2D structure from the previous lesson, but now the grid is shaped around the output matrix. Since C has rowsA rows and colsB columns, we need enough threads to cover those dimensions. The x-direction spans output columns, and the y-direction spans output rows.

    dim3 threadsPerBlock(16, 16);
    dim3 numBlocks((colsB + threadsPerBlock.x - 1) / threadsPerBlock.x,
                   (rowsA + threadsPerBlock.y - 1) / threadsPerBlock.y);

The round-up formula is important:

• (size + blockSize - 1) / blockSize ensures we launch enough blocks, even when the matrix size is not a perfect multiple of the block size.

For this tiny 2 x 2 output, one block is enough. Still, this exact formula also works for much larger rectangular matrices.

At this point, everything is ready. We launch the matrixMultiply kernel, check for launch errors using cudaGetLastError, wait for the device to finish with cudaDeviceSynchronize, and copy the result back to the host. Then, we compare each output value against h_expected. Notice that the check uses std::fabs from the <csmath> library with a small tolerance of 1e-5f, which is a good habit for floating-point results. Finally, we print a success message, show the top-left value h_C[0], free device memory with cudaFree, and return a status code.

    matrixMultiply<<<numBlocks, threadsPerBlock>>>(d_A, d_B, d_C, rowsA, colsA, colsB);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    CUDA_CHECK(cudaMemcpy(h_C, d_C, bytesC, cudaMemcpyDeviceToHost));

    bool success = true;
    for (int i = 0; i < rowsA * colsB; ++i) {
        if (std::fabs(h_C[i] - h_expected[i]) > 1e-5f) {
            success = false;
            break;
        }
    }

    std::cout << "Matrix multiplication verification: "
              << (success ? "SUCCESS" : "FAILURE") << std::endl;
    std::cout << "Top-left result: " << h_C[0] << " (Expected: 58.0)" << std::endl;

    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    return success ? 0 : 1;
}

The program produces:

Matrix multiplication verification: SUCCESS
Top-left result: 58 (Expected: 58.0)

This output tells us two useful things at once: the full verification passed and the first element matches the hand-computed value 58. Together, they provide strong evidence that the indexing and dot product logic are correct.

In this lesson, we turned 2D thread indexing into a real matrix operation. We defined matrix shapes for general rectangular multiplication, mapped each thread to one output element, used a loop over k to compute a dot product, launched a 2D grid based on the output matrix, and verified the result on the host.

This kernel is called naive for a reason: it is correct and clear, but not yet optimized. This is a great place to be while learning because we now understand the full flow from math to memory to execution. In the practice section ahead, we will reinforce this pattern so we can write the kernel confidently on our own.