CUDA GEMM: Naive Kernel to Coalesced

Write, profile, and optimize CUDA matrix multiplication kernels from scratch.

Introduction

30 Second Summary

Every time an AI generates an image, translates a sentence, or recommends a video, the GPU is doing one thing under the hood: multiplying massive matrices together. The speed of that single operation determines whether inference takes milliseconds or seconds.

In this project, you will write a matrix multiplication kernel in CUDA from scratch, profile it with NVIDIA's Nsight Compute, and optimize it based on real performance data.

What You'll Build

A complete CUDA project that implements single-precision GEMM (General Matrix Multiplication), benchmarks it against NVIDIA's cuBLAS library, and uses profiling evidence to guide optimization.

By the end of this project, you'll have:

  • A CMake-based CUDA project with two GEMM kernel implementations (naive and memory-coalesced).
  • A benchmarking harness that measures GFLOPS and compares your kernels against cuBLAS.
  • Nsight Compute profiling reports showing memory access patterns and proving your optimization works.
  • Secret Mission: Add a shared memory tiled kernel that further closes the performance gap to cuBLAS.

Do I need prior GPU programming experience?

No. This project assumes familiarity with C/C++ and a basic understanding of how matrix multiplication works (rows times columns). You will learn CUDA concepts as you go.

Before We Start

Before We Start...

Install CUDA Toolkit and Set Up the Project

To write and run code on your GPU, you need the CUDA Toolkit installed on your machine. The toolkit includes the nvcc compiler (which compiles your CUDA kernels), the cuBLAS library (NVIDIA's optimized math routines), and Nsight Compute (the profiler you will use to find performance bottlenecks).

In this step, you will install CUDA Toolkit 13.3, verify that your RTX 4060 is detected, install CMake for building the project, and create the project folder structure.

In this step, get ready to:
  • Install CUDA Toolkit 13.3 and verify your GPU is detected.
  • Install CMake and confirm it works from the command line.
  • Create the project directory structure for your GEMM kernels.
Download and install CUDA Toolkit 13.3

The CUDA Toolkit is a single installer that bundles everything you need to compile and profile GPU code on Windows.

  • Open your browser and go to the CUDA Toolkit download page:
  • https://developer.nvidia.com/cuda-downloads
  • Select your operating system as Windows, architecture as x86_64, and installer type as the Local installer.
  • Click **Download** and wait for the installer to finish downloading.
  • Run the downloaded installer and follow the on-screen prompts, accepting the default options.

What does the installer include?

The default installation gives you nvcc (the CUDA compiler), cuBLAS (the optimized linear algebra library), Nsight Compute (the GPU profiler), and the CUDA runtime libraries. You do not need to customize the installation.

The installer takes a few minutes

The Local installer is a large download (around 3 GB). Once downloaded, the installation itself takes a few minutes as it registers components with Visual Studio and updates your system PATH.

Verify the CUDA installation and GPU detection

After installation finishes, you need to open a fresh Command Prompt so that it picks up the new PATH entries the installer added.

  • Press the Windows key to open your search bar.
  • Type cmd and press Enter to open a new Command Prompt window.
  • Confirm the CUDA compiler is available by running this command:
nvcc -V

Check your output against the scenarios below.

✔️ I see release 13.3

Your output includes text like release 13.3. The CUDA compiler is installed and ready.

ⓧ I see an older version

If you see a version number lower than 13.3, the installer may not have overwritten your previous CUDA installation. Uninstall the older toolkit from Settings > Apps, then re-run the CUDA 13.3 installer you just downloaded.

After reinstalling, close your Command Prompt and open a brand new one before running nvcc -V again.

ⓧ Command not found

If you see a message saying the command is not recognized, the CUDA bin directory was not added to your system PATH. Try these fixes:

  • Close the current Command Prompt and open a brand new one. The PATH only updates in new sessions.
  • If a new Command Prompt still does not find nvcc, re-run the CUDA Toolkit installer and make sure you accept the default options.
  • Confirm your GPU is detected by running:
nvidia-smi

You should see a table listing your NVIDIA GeForce RTX 4060 along with driver and CUDA version information.

Not seeing your GPU listed?

Make sure your NVIDIA drivers are up to date. You can update them through GeForce Experience or by downloading the latest driver from the NVIDIA website. If nvidia-smi is not recognized at all, your GPU driver may not be installed correctly.

help me troubleshoot why nvidia-smi is not finding my RTX 4060 on Windows

Install CMake and create the project directory

CMake is the build system you will use to compile your CUDA project. It automatically finds the CUDA compiler and links the cuBLAS library for you.

  • Open your browser and go to the CMake download page:
  • https://cmake.org/download/
  • Download the Windows x64 installer (the .msi file).
  • Run the installer. When prompted, select the option to add CMake to the system PATH for all users.
  • Complete the installer with the default settings.

Why add CMake to PATH?

Adding CMake to the system PATH means you can run cmake from any Command Prompt window without needing to type the full installation path. If you skip this option, the cmake command will not be recognized.

  • Close your current Command Prompt and open a new one (press the Windows key, type cmd, press Enter).
  • Verify CMake is installed by running:
cmake --version

Check your output against the scenarios below.

✔️ I see a version number

You see output showing a CMake version (for example, 3.31 or newer). CMake is ready to use.

ⓧ Command not found

If the command is not recognized, the PATH was not updated during installation. Try these fixes:

  • Make sure you opened a brand new Command Prompt after installing CMake.
  • Re-run the CMake installer and select the option to add CMake to the system PATH for all users.
  • Move to your Desktop by running this command:
cd %USERPROFILE%\Desktop
  • Create the project folder and its subdirectories by running these commands:
mkdir cuda-gemm
cd cuda-gemm
mkdir src
mkdir build
  • Confirm the folders were created by running:
dir

You should see src and build listed as directories in the output.

Don't see both folders?

Make sure you ran each mkdir command after the cd cuda-gemm line. If you accidentally created them on the Desktop instead, delete them with rmdir src and rmdir build, then cd cuda-gemm and create them again inside the project folder.

help me fix my project directory structure for the cuda-gemm project

Your development environment is ready. Next up, you will write your first CUDA kernel that performs matrix multiplication on the GPU.

Write and Run a Naive GEMM Kernel

Matrix multiplication is the core operation behind nearly all GPU workloads. To understand how the GPU executes this work, you need to write the kernel yourself and see the raw performance numbers.

In this step, you will write a naive GEMM kernel where each GPU thread computes a single element of the output matrix. You will also build a benchmarking harness to measure performance in GFLOPS and validate correctness against a CPU reference.

In this step, get ready to:
  • Create the CMake build configuration and CUDA kernel source files.
  • Create a benchmarking harness that times the kernel and validates correctness.
  • Build and run the project, observing GFLOPS output for matrix multiplication.
Create the CMake configuration and kernel files

Every CUDA project needs a build system to compile .cu files with nvcc and link GPU libraries. CMake handles this automatically when you declare CUDA as a project language.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type VS Code and press Enter to open it.
  • In VS Code, click File in the top menu bar.
  • Click Open Folder.
  • Select the cuda-gemm folder on your Desktop, then click Open.
  • Create a new file called CMakeLists.txt in the root of your cuda-gemm folder and paste the following code:
cmake_minimum_required(VERSION 3.18)
project(cuda_gemm LANGUAGES CXX CUDA)

# Target your RTX 4060 (compute capability 8.9)
set(CMAKE_CUDA_ARCHITECTURES 89)

# Find the CUDA Toolkit (provides cuBLAS, cudart)
find_package(CUDAToolkit REQUIRED)

# Build the gemm executable from our source files
add_executable(gemm src/main.cu src/kernels.cu)

# Link against CUDA runtime and cuBLAS
target_link_libraries(gemm PRIVATE CUDA::cudart CUDA::cublas)

# Use C++17 and CUDA 17 standards
set_target_properties(gemm PROPERTIES
    CUDA_STANDARD 17
    CXX_STANDARD 17
)

What does this CMake configuration do?

This file tells CMake how to compile your CUDA project.

  • CMAKE_CUDA_ARCHITECTURES 89 generates GPU machine code specifically for your RTX 4060 (Ada Lovelace architecture).
  • find_package(CUDAToolkit REQUIRED) locates the CUDA Toolkit you installed, giving CMake access to cuBLAS and other libraries.
  • target_link_libraries links the CUDA runtime and cuBLAS library so you can use GPU memory management and NVIDIA's optimized matrix multiplication.
  • Save CMakeLists.txt.
  • Create a new file called kernels.cuh inside the src/ folder and paste the following code:
#ifndef KERNELS_CUH
#define KERNELS_CUH

// Naive GEMM: each thread computes one element of C
__global__ void sgemm_naive(int M, int N, int K,
                            const float *A, const float *B, float *C);

// Coalesced GEMM: fix memory access pattern for coalescing
__global__ void sgemm_coalesced(int M, int N, int K,
                                const float *A, const float *B, float *C);

#endif

What is a .cuh header file?

A .cuh file is a CUDA header. It declares kernel function signatures so other files can call them. The __global__ keyword marks a function as a GPU kernel that runs on the device but can be launched from the host (CPU).

You are declaring two kernels here. You will implement sgemm_naive now and use sgemm_coalesced in a later optimization step.

  • Save src/kernels.cuh.

Now you will write the actual kernel implementation. The naive GEMM kernel assigns one thread to compute one element of the output matrix C.

  • Create a new file called kernels.cu inside the src/ folder and paste the naive kernel implementation:
#include "kernels.cuh"

__global__ void sgemm_naive(int M, int N, int K,
                            const float *A, const float *B, float *C) {
    // Each thread computes one element of C
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < M && col < N) {
        float sum = 0.0f;
        // Dot product of row of A and column of B
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

How does the naive GEMM kernel work?

The GPU launches thousands of threads in a 2D grid. Each thread figures out which row and column of C it is responsible for using blockIdx and threadIdx.

That thread then computes a dot product: it walks across one row of matrix A and one column of matrix B, accumulating the sum. The result goes into one element of C. This is correct but not fast, because each thread independently reads from global memory without regard for how neighboring threads access memory.

Next, add the coalesced kernel below the naive kernel in the same file. You will use this kernel in a later step.

  • Add the following code below the sgemm_naive function in src/kernels.cu:
__global__ void sgemm_coalesced(int M, int N, int K,
                                const float *A, const float *B, float *C) {
    // Swap row/col mapping for coalesced memory access
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

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

Why add this kernel now?

The coalesced kernel swaps the order of row and col variable declarations. This changes which dimension maps to threadIdx.x. You will measure the performance difference in Step 5 after profiling with Nsight Compute reveals why memory access order matters.

  • Save src/kernels.cu.

VS Code showing red squiggles on CUDA keywords?

VS Code may not recognize __global__ or blockIdx by default. These are CUDA-specific keywords that nvcc understands. The red squiggles are cosmetic and do not affect compilation. You can ignore them or install the NVIDIA CUDA extension for VS Code.

help me fix CUDA syntax highlighting in VS Code

✔️ Awesome, I've got everything!

Great. Double check you have saved all three files: CMakeLists.txt, src/kernels.cuh, and src/kernels.cu.

ⓧ I'd like to double check the full code

Here is the complete src/kernels.cu file:

#include "kernels.cuh"

__global__ void sgemm_naive(int M, int N, int K,
                            const float *A, const float *B, float *C) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

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

__global__ void sgemm_coalesced(int M, int N, int K,
                                const float *A, const float *B, float *C) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}
Create the benchmarking harness

A kernel is useless without a way to measure it. The benchmarking harness allocates matrices on the GPU, times kernel execution with CUDA events, and validates the output against a known-correct CPU implementation.

  • Create a new file called main.cu inside the src/ folder. Start by adding the includes and utility functions:
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include "cublas_v2.h"
#include "kernels.cuh"

// Fill a matrix with random values between 0 and 1
void fill_matrix(float *mat, int size) {
    for (int i = 0; i < size; i++) {
        mat[i] = (float)(rand() % 100) / 100.0f;
    }
}

// CPU reference implementation for correctness validation
void cpu_gemm(int M, int N, int K, const float *A, const float *B, float *C) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

What do these functions do?

fill_matrix generates random test data so you have matrices to multiply.

cpu_gemm is a triple-nested loop that computes matrix multiplication on the CPU. It is slow but guaranteed correct. You will compare the GPU output against this to catch bugs in your kernel.

  • Add the timing harness function below cpu_gemm in src/main.cu:
// Time a kernel over multiple runs using CUDA events
float run_kernel(void (*launcher)(int, int, int, const float*, const float*, float*, dim3, dim3),
                 int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 dim3 grid, dim3 block, int warmup, int repeats) {
    // Warmup runs to fill caches and stabilize clocks
    for (int i = 0; i < warmup; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaDeviceSynchronize();

    // Timed runs using CUDA events for precise GPU timing
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

Why use CUDA events instead of wall-clock time?

CUDA events measure time on the GPU itself, not the CPU. This gives you precise kernel execution time without including CPU overhead or driver latency. The warmup runs ensure the GPU clocks have ramped up before you start measuring.

  • Add the kernel launcher wrappers and validation function below run_kernel:
// Wrapper that launches the naive kernel with grid/block config
void launch_naive(int M, int N, int K,
                  const float *A, const float *B, float *C,
                  dim3 grid, dim3 block) {
    sgemm_naive<<<grid, block>>>(M, N, K, A, B, C);
}

// Wrapper that launches the coalesced kernel
void launch_coalesced(int M, int N, int K,
                      const float *A, const float *B, float *C,
                      dim3 grid, dim3 block) {
    sgemm_coalesced<<<grid, block>>>(M, N, K, A, B, C);
}

// Compare GPU output against CPU reference, element by element
bool validate(const float *ref, const float *test, int size, float tol) {
    for (int i = 0; i < size; i++) {
        if (fabs(ref[i] - test[i]) > tol) {
            printf("Mismatch at index %d: ref=%.6f test=%.6f\n", i, ref[i], test[i]);
            return false;
        }
    }
    return true;
}

What are the launcher wrappers for?

The <<<grid, block>>> syntax is CUDA's kernel launch syntax. Each launcher wraps one kernel so that run_kernel can accept any launcher as a function pointer. This pattern lets you benchmark different kernels with the same timing code.

validate compares every element of the GPU result against the CPU reference. A tolerance of 1e-3 accounts for floating-point rounding differences between CPU and GPU arithmetic.

  • Add the main function below validate to complete the harness:
int main() {
    int sizes[] = {512, 1024, 2048};
    int num_sizes = 3;
    int warmup = 3;
    int repeats = 10;

    // Naive-only benchmark for this step
    for (int s = 0; s < num_sizes; s++) {
        int M = sizes[s], N = sizes[s], K = sizes[s];
        size_t bytes = M * N * sizeof(float);
        double flops = 2.0 * M * N * K;

        // Allocate host memory
        float *h_A = (float *)malloc(M * K * sizeof(float));
        float *h_B = (float *)malloc(K * N * sizeof(float));
        float *h_C = (float *)malloc(bytes);

        fill_matrix(h_A, M * K);
        fill_matrix(h_B, K * N);

        // Allocate device (GPU) memory
        float *d_A, *d_B, *d_C;
        cudaMalloc(&d_A, M * K * sizeof(float));
        cudaMalloc(&d_B, K * N * sizeof(float));
        cudaMalloc(&d_C, bytes);

        // Copy matrices from host to device
        cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);

        // Configure thread block (32x32 = 1024 threads per block)
        dim3 block(32, 32);
        dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

        // Run and time the naive kernel
        float naive_ms = run_kernel(launch_naive, M, N, K, d_A, d_B, d_C,
                                    grid, block, warmup, repeats);
        double naive_gflops = (flops / (naive_ms / 1000.0)) / 1e9;
        printf("Size %d: Naive = %.1f GFLOPS\n", M, naive_gflops);

        // Validate against CPU for smallest size
        if (s == 0) {
            float *h_ref = (float *)malloc(bytes);
            cpu_gemm(M, N, K, h_A, h_B, h_ref);
            cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
            if (validate(h_ref, h_C, M * N, 1e-3f)) {
                printf("Correctness check PASSED\n");
            } else {
                printf("Correctness check FAILED\n");
            }
            free(h_ref);
        }

        free(h_A); free(h_B); free(h_C);
        cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    }
    return 0;
}

What does the main function do?

  • It loops over three matrix sizes (512, 1024, 2048) to show how performance scales.
  • For each size, it allocates matrices on both CPU and GPU, copies data to the GPU, and launches the naive kernel.
  • dim3 block(32, 32) creates blocks of 1024 threads (the maximum per block on your GPU). The grid dimensions ensure full coverage of the output matrix.
  • GFLOPS is computed as (2 * M * N * K) / (time in seconds) / 1e9. The factor of 2 counts one multiply and one add per inner-loop iteration.
  • Correctness validation runs only on the 512x512 case to keep things fast.
  • Save src/main.cu.

Seeing errors about cublas_v2.h not found?

This header ships with the CUDA Toolkit. If your editor shows it as missing, the code will still compile correctly with CMake because CMake finds the toolkit include paths automatically. You can ignore editor warnings about this include.

help me fix include errors for cublas_v2.h in my CUDA project

✔️ Awesome, I've got everything!

Great. Make sure src/main.cu is saved and contains all four sections: includes/utilities, run_kernel, launchers/validate, and main.

ⓧ I'd like to double check the full code

Here is the complete src/main.cu file:

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include "cublas_v2.h"
#include "kernels.cuh"

void fill_matrix(float *mat, int size) {
    for (int i = 0; i < size; i++) {
        mat[i] = (float)(rand() % 100) / 100.0f;
    }
}

void cpu_gemm(int M, int N, int K, const float *A, const float *B, float *C) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

float run_kernel(void (*launcher)(int, int, int, const float*, const float*, float*, dim3, dim3),
                 int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 dim3 grid, dim3 block, int warmup, int repeats) {
    for (int i = 0; i < warmup; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaDeviceSynchronize();

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

void launch_naive(int M, int N, int K,
                  const float *A, const float *B, float *C,
                  dim3 grid, dim3 block) {
    sgemm_naive<<<grid, block>>>(M, N, K, A, B, C);
}

void launch_coalesced(int M, int N, int K,
                      const float *A, const float *B, float *C,
                      dim3 grid, dim3 block) {
    sgemm_coalesced<<<grid, block>>>(M, N, K, A, B, C);
}

bool validate(const float *ref, const float *test, int size, float tol) {
    for (int i = 0; i < size; i++) {
        if (fabs(ref[i] - test[i]) > tol) {
            printf("Mismatch at index %d: ref=%.6f test=%.6f\n", i, ref[i], test[i]);
            return false;
        }
    }
    return true;
}

int main() {
    int sizes[] = {512, 1024, 2048};
    int num_sizes = 3;
    int warmup = 3;
    int repeats = 10;

    // Naive-only benchmark for this step
    for (int s = 0; s < num_sizes; s++) {
        int M = sizes[s], N = sizes[s], K = sizes[s];
        size_t bytes = M * N * sizeof(float);
        double flops = 2.0 * M * N * K;

        float *h_A = (float *)malloc(M * K * sizeof(float));
        float *h_B = (float *)malloc(K * N * sizeof(float));
        float *h_C = (float *)malloc(bytes);

        fill_matrix(h_A, M * K);
        fill_matrix(h_B, K * N);

        float *d_A, *d_B, *d_C;
        cudaMalloc(&d_A, M * K * sizeof(float));
        cudaMalloc(&d_B, K * N * sizeof(float));
        cudaMalloc(&d_C, bytes);

        cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);

        dim3 block(32, 32);
        dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

        float naive_ms = run_kernel(launch_naive, M, N, K, d_A, d_B, d_C,
                                    grid, block, warmup, repeats);
        double naive_gflops = (flops / (naive_ms / 1000.0)) / 1e9;
        printf("Size %d: Naive = %.1f GFLOPS\n", M, naive_gflops);

        if (s == 0) {
            float *h_ref = (float *)malloc(bytes);
            cpu_gemm(M, N, K, h_A, h_B, h_ref);
            cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
            if (validate(h_ref, h_C, M * N, 1e-3f)) {
                printf("Correctness check PASSED\n");
            } else {
                printf("Correctness check FAILED\n");
            }
            free(h_ref);
        }

        free(h_A); free(h_B); free(h_C);
        cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    }
    return 0;
}
Build and run the project

With all source files in place, you can now use CMake to configure and build your project into an executable.

  • Open a terminal in VS Code by pressing Ctrl+` (backtick key).
  • Make sure you are in the cuda-gemm folder (the prompt should show this path).
  • Configure the CMake build system by running:
cmake -S . -B build

What does this command do?

cmake -S . -B build tells CMake to read your CMakeLists.txt from the current directory (-S .) and generate build files into the build/ directory (-B build). CMake will detect your Visual Studio compiler and the CUDA Toolkit automatically.

  • CMake will print output showing it found the CUDA compiler and the CUDAToolkit. Look for lines mentioning nvcc and Build files have been written to.
  • Build the project in Release mode by running:
cmake --build build --config Release
  • Wait for compilation to finish. You should see a line indicating the build succeeded with no errors.
  • Run the executable by running:
build\Release\gemm.exe

You should see output like this (your exact GFLOPS numbers will vary):

Size 512: Naive = 45.2 GFLOPS
Correctness check PASSED
Size 1024: Naive = 52.8 GFLOPS
Size 2048: Naive = 55.1 GFLOPS

What do these numbers mean?

Your naive kernel likely achieves somewhere between 30-80 GFLOPS on an RTX 4060. The theoretical peak for this GPU is much higher. In the next step, you will compare against cuBLAS (NVIDIA's optimized implementation) to see exactly how much performance you are leaving on the table.

The correctness check passing confirms your kernel computes the same result as the CPU reference. This is critical before you start optimizing.

Build failed with errors?

  • If CMake cannot find the CUDA compiler, make sure nvcc is on your system PATH. Open a new Command Prompt and run nvcc -V to verify.
  • If you see linker errors about cuBLAS, double check that your CMakeLists.txt has find_package(CUDAToolkit REQUIRED) and target_link_libraries(gemm PRIVATE CUDA::cudart CUDA::cublas).
  • If the correctness check fails, make sure your kernel indexing matches the code above exactly. A swapped row/column index will produce wrong results.

help me debug my CUDA GEMM build errors

Your kernel is running correctly and you have a baseline measurement. Next up, you will benchmark against cuBLAS to see how far your naive kernel is from NVIDIA's hand-tuned implementation.

Benchmark Against cuBLAS

Your naive kernel runs and prints GFLOPS numbers. But how do you know if those numbers are good or terrible?

You need a reference point. cuBLAS is NVIDIA's hand-tuned linear algebra library that ships with the CUDA Toolkit. It represents the practical ceiling for GEMM performance on your GPU. In this step, you will benchmark your kernels against cuBLAS so you can see exactly how much performance you are leaving on the table.

In this step, get ready to:
  • Add a cuBLAS SGEMM benchmark function using the cuBLAS v2 API.
  • Print a comparison table showing naive, coalesced, and cuBLAS GFLOPS across multiple matrix sizes.
Add the cuBLAS benchmark function

cuBLAS expects matrices in column-major (Fortran-style) storage. Your kernels use row-major (C-style) storage. To get the correct result without transposing your data, you swap the order of arguments and the dimensions in the cublasSgemm call.

  • In src/main.cu, find the launch_coalesced function. Add the following run_cublas function directly below it:
float run_cublas(cublasHandle_t handle, int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 int warmup, int repeats) {
    float alpha = 1.0f;
    float beta = 0.0f;

    // Warmup runs to fill caches and trigger JIT compilation
    for (int i = 0; i < warmup; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaDeviceSynchronize();

    // Timed runs using CUDA events
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

Why are the arguments swapped?

cuBLAS uses column-major storage by default. Your matrices are stored row-major. The mathematical trick is: computing C = A * B in row-major is equivalent to computing C_T = B_T * A_T in column-major.

So you pass d_B first, then d_A, with swapped dimensions (N, M instead of M, N). The CUBLAS_OP_N flags mean "no transpose" since you are handling the layout swap through argument order.

You won't see a change yet. The next chunk updates main to call this function.

Update main to print a comparison table

Now you will replace the entire main function to run all three implementations (naive, coalesced, cuBLAS) across sizes 512, 1024, and 2048, printing a formatted table.

  • In src/main.cu, find the existing int main() function and replace it entirely with the following code:
int main() {
    int sizes[] = {512, 1024, 2048};
    int num_sizes = 3;
    int warmup = 3;
    int repeats = 10;

    // Initialize cuBLAS library handle
    cublasHandle_t cublas_handle;
    cublasCreate(&cublas_handle);

    // Print table header
    printf("%-6s %-12s %-12s %-12s %-10s %-10s\n",
           "Size", "Naive(GFLOPS)", "Coal.(GFLOPS)", "cuBLAS(GFLOPS)",
           "Naive%", "Coal.%");
    printf("--------------------------------------------------------------\n");

    for (int s = 0; s < num_sizes; s++) {
        int M = sizes[s], N = sizes[s], K = sizes[s];
        size_t bytes = M * N * sizeof(float);
        double flops = 2.0 * M * N * K;

        float *h_A = (float *)malloc(M * K * sizeof(float));
        float *h_B = (float *)malloc(K * N * sizeof(float));
        float *h_C = (float *)malloc(bytes);

        fill_matrix(h_A, M * K);
        fill_matrix(h_B, K * N);

        float *d_A, *d_B, *d_C;
        cudaMalloc(&d_A, M * K * sizeof(float));
        cudaMalloc(&d_B, K * N * sizeof(float));
        cudaMalloc(&d_C, bytes);

        cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);

        dim3 block(32, 32);
        dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

        // Benchmark naive kernel
        float naive_ms = run_kernel(launch_naive, M, N, K, d_A, d_B, d_C,
                                    grid, block, warmup, repeats);
        double naive_gflops = (flops / (naive_ms / 1000.0)) / 1e9;

        // Benchmark coalesced kernel
        float coal_ms = run_kernel(launch_coalesced, M, N, K, d_A, d_B, d_C,
                                   grid, block, warmup, repeats);
        double coal_gflops = (flops / (coal_ms / 1000.0)) / 1e9;

        // Benchmark cuBLAS
        float cublas_ms = run_cublas(cublas_handle, M, N, K, d_A, d_B, d_C,
                                     warmup, repeats);
        double cublas_gflops = (flops / (cublas_ms / 1000.0)) / 1e9;

        // Print row with percentages of cuBLAS performance
        printf("%-6d %-12.1f %-12.1f %-12.1f %-9.1f%% %-9.1f%%\n",
               M, naive_gflops, coal_gflops, cublas_gflops,
               (naive_gflops / cublas_gflops) * 100.0,
               (coal_gflops / cublas_gflops) * 100.0);

        // Validate against CPU for smallest size
        if (s == 0) {
            float *h_ref = (float *)malloc(bytes);
            cpu_gemm(M, N, K, h_A, h_B, h_ref);
            cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
            if (validate(h_ref, h_C, M * N, 1e-3f)) {
                printf("Correctness check PASSED\n");
            } else {
                printf("Correctness check FAILED\n");
            }
            free(h_ref);
        }

        free(h_A); free(h_B); free(h_C);
        cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    }

    // Clean up cuBLAS handle
    cublasDestroy(cublas_handle);
    return 0;
}

What does this code do?

  • cublasCreate(&cublas_handle) initializes the cuBLAS library and returns a handle you pass to every cuBLAS call.
  • The loop now runs naive, coalesced, AND cuBLAS for each matrix size (512, 1024, 2048).
  • The percentage columns show how much of cuBLAS performance each kernel achieves. This is the standard way GPU engineers report optimization progress.
  • cublasDestroy(cublas_handle) frees the library resources at the end.
  • Save src/main.cu.
  • Rebuild and run the project by running these commands:
cd build
cmake --build . --config Release
Release\gemm.exe

You should see a formatted table with columns for Size, Naive GFLOPS, Coalesced GFLOPS, cuBLAS GFLOPS, and the percentage of cuBLAS that each kernel achieves. The naive kernel typically reaches 1-5% of cuBLAS performance.

Seeing a linker error about cublasSgemm?

Make sure your CMakeLists.txt has CUDA::cublas in the target_link_libraries line. Also verify you have find_package(CUDAToolkit REQUIRED) before the executable definition.

If you see a runtime error instead of a table, check that #include "cublas_v2.h" is at the top of src/main.cu.

Help me debug my cuBLAS linker error

Why is the gap so large?

cuBLAS has been hand-optimized over 15+ years. It uses shared memory tiling, register-level blocking, warp-level matrix operations, and hardware-specific instruction scheduling. Your naive kernel makes a new global memory request for every single element in the inner loop. The performance gap tells you exactly how much optimization headroom exists.

✔️ Awesome, I've got everything!

Great. Double check you have saved src/main.cu and your program prints the comparison table with all three implementations.

ⓧ I'd like to double check the full code

Here is the complete src/main.cu at this point:

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include "cublas_v2.h"
#include "kernels.cuh"

void fill_matrix(float *mat, int size) {
    for (int i = 0; i < size; i++) {
        mat[i] = (float)(rand() % 100) / 100.0f;
    }
}

void cpu_gemm(int M, int N, int K, const float *A, const float *B, float *C) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

float run_kernel(void (*launcher)(int, int, int, const float*, const float*, float*,
                                  dim3, dim3),
                 int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 dim3 grid, dim3 block, int warmup, int repeats) {
    for (int i = 0; i < warmup; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaDeviceSynchronize();

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

void launch_naive(int M, int N, int K,
                  const float *A, const float *B, float *C,
                  dim3 grid, dim3 block) {
    sgemm_naive<<<grid, block>>>(M, N, K, A, B, C);
}

void launch_coalesced(int M, int N, int K,
                      const float *A, const float *B, float *C,
                      dim3 grid, dim3 block) {
    sgemm_coalesced<<<grid, block>>>(M, N, K, A, B, C);
}

float run_cublas(cublasHandle_t handle, int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 int warmup, int repeats) {
    float alpha = 1.0f;
    float beta = 0.0f;

    for (int i = 0; i < warmup; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaDeviceSynchronize();

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

bool validate(const float *ref, const float *test, int size, float tol) {
    for (int i = 0; i < size; i++) {
        if (fabs(ref[i] - test[i]) > tol) {
            printf("Mismatch at index %d: ref=%.6f test=%.6f\n", i, ref[i], test[i]);
            return false;
        }
    }
    return true;
}

int main() {
    int sizes[] = {512, 1024, 2048};
    int num_sizes = 3;
    int warmup = 3;
    int repeats = 10;

    cublasHandle_t cublas_handle;
    cublasCreate(&cublas_handle);

    printf("%-6s %-12s %-12s %-12s %-10s %-10s\n",
           "Size", "Naive(GFLOPS)", "Coal.(GFLOPS)", "cuBLAS(GFLOPS)",
           "Naive%", "Coal.%");
    printf("--------------------------------------------------------------\n");

    for (int s = 0; s < num_sizes; s++) {
        int M = sizes[s], N = sizes[s], K = sizes[s];
        size_t bytes = M * N * sizeof(float);
        double flops = 2.0 * M * N * K;

        float *h_A = (float *)malloc(M * K * sizeof(float));
        float *h_B = (float *)malloc(K * N * sizeof(float));
        float *h_C = (float *)malloc(bytes);

        fill_matrix(h_A, M * K);
        fill_matrix(h_B, K * N);

        float *d_A, *d_B, *d_C;
        cudaMalloc(&d_A, M * K * sizeof(float));
        cudaMalloc(&d_B, K * N * sizeof(float));
        cudaMalloc(&d_C, bytes);

        cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);

        dim3 block(32, 32);
        dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

        float naive_ms = run_kernel(launch_naive, M, N, K, d_A, d_B, d_C,
                                    grid, block, warmup, repeats);
        double naive_gflops = (flops / (naive_ms / 1000.0)) / 1e9;

        float coal_ms = run_kernel(launch_coalesced, M, N, K, d_A, d_B, d_C,
                                   grid, block, warmup, repeats);
        double coal_gflops = (flops / (coal_ms / 1000.0)) / 1e9;

        float cublas_ms = run_cublas(cublas_handle, M, N, K, d_A, d_B, d_C,
                                     warmup, repeats);
        double cublas_gflops = (flops / (cublas_ms / 1000.0)) / 1e9;

        printf("%-6d %-12.1f %-12.1f %-12.1f %-9.1f%% %-9.1f%%\n",
               M, naive_gflops, coal_gflops, cublas_gflops,
               (naive_gflops / cublas_gflops) * 100.0,
               (coal_gflops / cublas_gflops) * 100.0);

        if (s == 0) {
            float *h_ref = (float *)malloc(bytes);
            cpu_gemm(M, N, K, h_A, h_B, h_ref);
            cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
            if (validate(h_ref, h_C, M * N, 1e-3f)) {
                printf("Correctness check PASSED\n");
            } else {
                printf("Correctness check FAILED\n");
            }
            free(h_ref);
        }

        free(h_A); free(h_B); free(h_C);
        cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    }

    cublasDestroy(cublas_handle);
    return 0;
}

You now have hard evidence that your kernels are far below the performance ceiling. The next step is to use Nsight Compute to understand exactly why. The profiler will reveal the memory access pattern problem that accounts for most of this gap.

Profile with Nsight Compute

You now have GFLOPS numbers for all three implementations. The naive kernel is dramatically slower than cuBLAS, but the numbers alone don't explain why.

Nsight Compute is NVIDIA's GPU kernel profiler. It captures hardware-level counters that reveal exactly where your kernel's time goes. In this step, you'll profile the naive kernel and find the specific memory access problem that's killing performance.

In this step, get ready to:
  • Profile the naive kernel with the Nsight Compute CLI.
  • Open the profiling report and identify the memory coalescing bottleneck.
Profile the naive kernel with Nsight Compute CLI

The ncu CLI instruments your kernel at the hardware level, replaying it multiple times to collect different sets of performance counters. The result is a detailed report file you can explore in the Nsight Compute UI.

Why replay the kernel multiple times?

GPUs have a limited number of hardware counters, so not all metrics can be captured in a single pass. Nsight Compute automatically re-runs the kernel across multiple passes to gather a complete picture. This takes longer than a normal run, but produces far richer data.

  • Right-click on Command Prompt in your Start menu and select Run as administrator.
  • Navigate to your project folder by running:
cd %USERPROFILE%\Desktop\cuda-gemm
  • Profile the naive kernel with a full metric collection by running this command:
ncu --set full -o naive_profile build\gemm.exe

What does this command do?

  • ncu launches the Nsight Compute CLI profiler.
  • --set full collects the full set of profiling sections, including memory workload analysis.
  • -o naive_profile saves the report to a file called naive_profile.ncu-rep in the current directory.
  • build\gemm.exe is the compiled binary to profile.

The profiler will run your program multiple times. This can take a few minutes since the full set requires many replay passes.

  • Wait for the profiler to finish. You should see a final line indicating the report file was written.
  • Confirm the report exists by running:
dir naive_profile.ncu-rep

You should see naive_profile.ncu-rep listed in the output.

ncu is not recognized?

If ncu is not found, the Nsight Compute binary directory is not in your system PATH. The default install location is C:\Program Files\NVIDIA Corporation\Nsight Compute <version>.

As a quick workaround, search your file system for ncu.exe inside the NVIDIA Corporation folder, then run the command using the full path to that executable.

Also make sure you opened the Command Prompt as Administrator. On Windows with WDDM drivers, full profiling access requires elevated privileges.

help me fix the ncu not recognized error on my Windows machine

Open the report and identify the coalescing problem

The raw report file contains hundreds of metrics. The Nsight Compute UI organizes them into navigable sections by performance bottleneck category.

  • Open the report in the Nsight Compute UI by running:
ncu-ui naive_profile.ncu-rep

Alternative way to open the report

You can also launch ncu-ui on its own, close the Start Activity dialog, then open the report through File > Open.

  • In the Nsight Compute UI, you will see a list of all profiled kernel launches. Click on the sgemm_naive kernel entry to open its detailed results.
  • Scroll down to the Memory Workload Analysis section.

What is Memory Workload Analysis?

This section shows whether memory is a limiting factor for your kernel. It measures hardware unit utilization, communication bandwidth, and memory instruction throughput.

The key metric to look for is sectors per request for global load operations. A sector is an aligned 32-byte chunk of memory. For a perfectly coalesced access (32 threads in a warp each reading consecutive 4-byte floats), a single request loads exactly 4 sectors (128 bytes = 32 threads x 4 bytes).

  • In the Memory Workload Analysis, look at the Global Memory subsection. Find the sectors-per-request value for global loads.

A high sectors-per-request ratio means threads in the same warp are accessing scattered memory addresses. The GPU fetches many more sectors than necessary for each request.

Why is the naive kernel uncoalesced?

Look at the inner loop in sgemm_naive:

sum += A[row * K + k] * B[k * N + col];

In this kernel, threadIdx.x maps to col via blockIdx.x * blockDim.x + threadIdx.x. Consecutive threads (incrementing threadIdx.x) have different col values but share the same row. When reading A[row * K + k], every thread in the warp reads the exact same element of A.

This broadcast pattern generates excessive sector traffic that Nsight Compute flags as inefficient global memory access. The sectors-per-request metric quantifies this waste and gives you a concrete number to improve against in the next step.

  • Write down the sectors-per-request value you see for global loads in the naive kernel report. This is your baseline measurement that proves the bottleneck.

Report looks empty or sections are missing?

Make sure you ran ncu with --set full. Without this flag, only the basic set of sections is collected, which may not include the full Memory Workload Analysis.

If the report shows zero kernel launches, the profiler may have exited before your kernel ran. Ensure the Command Prompt has administrator privileges and re-run the profiling command.

my Nsight Compute report is missing memory metrics or looks empty

You now have profiling evidence that pinpoints the exact bottleneck. Next up, you will fix this access pattern with a coalesced kernel and measure the improvement.

Optimize with Memory Coalescing

In the last step, you profiled the naive kernel and pinpointed uncoalesced memory access as the primary bottleneck.

Now it's time to confirm that your coalesced kernel actually fixes the problem. You'll rebuild the project, run the full benchmark to measure the speedup, profile the coalesced kernel with Nsight Compute to verify improved memory patterns, and then document your findings in a README.

In this step, get ready to:
  • Run the benchmark and compare naive, coalesced, and cuBLAS GFLOPS.
  • Profile the coalesced kernel and confirm lower sectors-per-request.
  • Create a README and push the project to GitHub.
Rebuild and run the full benchmark

Your coalesced kernel swaps the row and column index assignments so that consecutive threads (incrementing threadIdx.x) access consecutive memory addresses. This single change aligns memory requests within a warp to the same cache lines, dramatically reducing the number of memory transactions the GPU needs.

  • Rebuild the project by running these commands in your Command Prompt from inside the cuda-gemm folder:
cd build
cmake --build . --config Release

Why rebuild?

Even though your code hasn't changed since last step, rebuilding ensures the binary is fresh and matches your current source files exactly. CMake's incremental build will detect nothing has changed and finish in seconds.

  • Run the benchmark by executing:
Release\gemm.exe
  • You should see a table with six columns: Size, Naive(GFLOPS), Coal.(GFLOPS), cuBLAS(GFLOPS), Naive%, and Coal.%.
  • The coalesced kernel should show 2-10x higher GFLOPS than the naive kernel, depending on matrix size.
  • You should also see Correctness check PASSED confirming the results are mathematically correct.

What do these numbers mean?

The Naive% and Coal.% columns show what fraction of cuBLAS performance each kernel achieves. The naive kernel typically hits 1-5% of cuBLAS. The coalesced kernel should reach a higher percentage because it eliminates wasted memory bandwidth from strided access patterns.

cuBLAS itself uses shared memory tiling, register blocking, and warp-level matrix operations that go far beyond what either of your kernels do yet.

Not seeing a speedup?

Make sure you rebuilt successfully and that the executable you are running is Release\gemm.exe (not a Debug build). Debug mode disables compiler optimizations that affect GPU code generation.

If both kernels show nearly identical numbers, confirm that your sgemm_naive and sgemm_coalesced functions have different index mappings in src/kernels.cu.

help me debug why my naive and coalesced kernels show the same GFLOPS

Profile the coalesced kernel with Nsight Compute

Numbers alone tell you there is a speedup, but profiling proves why. By comparing the sectors-per-request metric between your naive and coalesced kernels, you get hardware-level evidence that the coalescing fix reduced memory transactions.

  • Profile just the coalesced kernel by running this command (you may need to run Command Prompt as Administrator):
ncu --set full -k regex:coalesced -o coalesced_profile Release\gemm.exe

What does this command do?

  • ncu is the Nsight Compute CLI profiler.
  • --set full collects the full set of profiling sections (memory, compute, occupancy).
  • -k regex:coalesced filters to only profile kernels whose name contains "coalesced". This skips the naive and cuBLAS kernels.
  • -o coalesced_profile saves the report as coalesced_profile.ncu-rep.
  • Open the profiling report in the Nsight Compute UI:
ncu-ui coalesced_profile.ncu-rep
  • In the report, navigate to the Memory Workload Analysis section.
  • Look at the sectors-per-request metric for global loads.
  • Compare this value to what you documented from the naive kernel's profile in the previous step. The coalesced kernel should show a lower number, confirming that consecutive threads are now hitting the same cache lines.

What is a sector?

A sector is a 32-byte chunk of a cache line. When a warp of 32 threads requests memory, the GPU groups those requests by which sectors they hit. Fewer sectors-per-request means the warp's memory addresses are tightly packed. A perfectly coalesced 128-byte load (4 bytes per thread times 32 threads) hits exactly 4 sectors in one transaction.

Nsight Compute won't profile?

On Windows with WDDM display drivers, you may need to run Command Prompt as Administrator. Right-click Command Prompt and select Run as administrator, then re-run the ncu command.

help me troubleshoot Nsight Compute profiling access on Windows

Create a README and push to GitHub

A performance-focused GPU project isn't complete without documentation showing the optimization journey. Your README should tell the story: what you measured, what you found, and what you improved.

  • Create a file called README.md in the root of your cuda-gemm folder.
  • Add the following content to your README, replacing the placeholder numbers with your actual benchmark results:
# CUDA GEMM From Scratch

Single-precision matrix multiplication (SGEMM) implemented in CUDA,
benchmarked against cuBLAS, and profiled with Nsight Compute.

## Performance Results

| Size | Naive (GFLOPS) | Coalesced (GFLOPS) | cuBLAS (GFLOPS) | Naive % | Coalesced % |
|------|----------------|--------------------|-----------------|---------|-------------|
| 512  | [[NAIVE_512="your naive 512 GFLOPS"]]  | [[COAL_512="your coalesced 512 GFLOPS"]]  | [[CUBLAS_512="your cuBLAS 512 GFLOPS"]]  | [[NAIVE_PCT_512="naive %"]]  | [[COAL_PCT_512="coalesced %"]]  |
| 1024 | [[NAIVE_1024="your naive 1024 GFLOPS"]] | [[COAL_1024="your coalesced 1024 GFLOPS"]] | [[CUBLAS_1024="your cuBLAS 1024 GFLOPS"]] | [[NAIVE_PCT_1024="naive %"]] | [[COAL_PCT_1024="coalesced %"]] |
| 2048 | [[NAIVE_2048="your naive 2048 GFLOPS"]] | [[COAL_2048="your coalesced 2048 GFLOPS"]] | [[CUBLAS_2048="your cuBLAS 2048 GFLOPS"]] | [[NAIVE_PCT_2048="naive %"]] | [[COAL_PCT_2048="coalesced %"]] |

## Optimization: Memory Coalescing

The naive kernel maps `threadIdx.x` to the row index, causing threads in
the same warp to access memory addresses that are N elements apart (strided
access). This forces the GPU to issue many separate memory transactions.

The coalesced kernel swaps the mapping so `threadIdx.x` maps to the column
index. Since matrices are stored in row-major order, consecutive columns
are consecutive in memory. Threads in the same warp now access adjacent
addresses, and the GPU can service the entire warp with fewer transactions.

### Nsight Compute Evidence

- **Naive kernel**: High sectors-per-request (uncoalesced global loads)
- **Coalesced kernel**: Lower sectors-per-request (coalesced global loads)

## Next Steps

- Shared memory tiling to reduce redundant global memory reads
- Register-level optimizations for instruction-level parallelism
- Warp-level matrix operations (WMMA/Tensor Cores)

## Build

```bash
mkdir build && cd build
cmake ..
cmake --build . --config Release
Release\gemm.exe
```

Why document the optimization rationale?

GPU engineering hiring managers look for candidates who can explain WHY a change improved performance, not just that it did. Your README demonstrates you understand the hardware memory model and can make data-driven optimization decisions.

  • Save README.md.
  • Initialize a Git repository and push to GitHub by running these commands:
cd ..
git init
git add .
git commit -m "CUDA GEMM: naive + coalesced kernels with benchmarks"
git branch -M main
git remote add origin [[YOUR_REPO_URL="https://github.com/YOUR_USERNAME/cuda-gemm.git"]]
git push -u origin main

What should I push?

Push your source code (src/, CMakeLists.txt, README.md) but not your build/ folder or .ncu-rep profiling reports. Add a .gitignore file with build/ and *.ncu-rep if you want to keep the repository clean.

Git push rejected?

Make sure you created the repository on GitHub first (an empty repo with no README). If you initialized with a README on GitHub, use git pull --rebase origin main before pushing.

help me push my CUDA project to GitHub

✔️ Awesome, I've got everything!

Great. Double check you've saved all files and that your GitHub push completed successfully.

ⓧ I'd like to double check the full code

Here is the complete src/kernels.cu:

#include "kernels.cuh"

__global__ void sgemm_naive(int M, int N, int K,
                            const float *A, const float *B, float *C) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

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

__global__ void sgemm_coalesced(int M, int N, int K,
                                const float *A, const float *B, float *C) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

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

Here is the complete src/kernels.cuh:

#ifndef KERNELS_CUH
#define KERNELS_CUH

// Naive GEMM: each thread computes one element of C
__global__ void sgemm_naive(int M, int N, int K,
                            const float *A, const float *B, float *C);

// Coalesced GEMM: fix memory access pattern for coalescing
__global__ void sgemm_coalesced(int M, int N, int K,
                                const float *A, const float *B, float *C);

#endif

Here is the complete src/main.cu:

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include "cublas_v2.h"
#include "kernels.cuh"

void fill_matrix(float *mat, int size) {
    for (int i = 0; i < size; i++) {
        mat[i] = (float)(rand() % 100) / 100.0f;
    }
}

void cpu_gemm(int M, int N, int K, const float *A, const float *B, float *C) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

float run_kernel(void (*launcher)(int, int, int, const float*, const float*, float*,
                                  dim3, dim3),
                 int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 dim3 grid, dim3 block, int warmup, int repeats) {
    for (int i = 0; i < warmup; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaDeviceSynchronize();

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        launcher(M, N, K, d_A, d_B, d_C, grid, block);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

void launch_naive(int M, int N, int K,
                  const float *A, const float *B, float *C,
                  dim3 grid, dim3 block) {
    sgemm_naive<<<grid, block>>>(M, N, K, A, B, C);
}

void launch_coalesced(int M, int N, int K,
                      const float *A, const float *B, float *C,
                      dim3 grid, dim3 block) {
    sgemm_coalesced<<<grid, block>>>(M, N, K, A, B, C);
}

float run_cublas(cublasHandle_t handle, int M, int N, int K,
                 const float *d_A, const float *d_B, float *d_C,
                 int warmup, int repeats) {
    float alpha = 1.0f;
    float beta = 0.0f;

    for (int i = 0; i < warmup; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaDeviceSynchronize();

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);
    for (int i = 0; i < repeats; i++) {
        cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                    N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    return ms / repeats;
}

bool validate(const float *ref, const float *test, int size, float tol) {
    for (int i = 0; i < size; i++) {
        if (fabs(ref[i] - test[i]) > tol) {
            printf("Mismatch at index %d: ref=%.6f test=%.6f\n", i, ref[i], test[i]);
            return false;
        }
    }
    return true;
}

int main() {
    int sizes[] = {512, 1024, 2048};
    int num_sizes = 3;
    int warmup = 3;
    int repeats = 10;

    cublasHandle_t cublas_handle;
    cublasCreate(&cublas_handle);

    printf("%-6s %-12s %-12s %-12s %-10s %-10s\n",
           "Size", "Naive(GFLOPS)", "Coal.(GFLOPS)", "cuBLAS(GFLOPS)",
           "Naive%", "Coal.%");
    printf("--------------------------------------------------------------\n");

    for (int s = 0; s < num_sizes; s++) {
        int M = sizes[s], N = sizes[s], K = sizes[s];
        size_t bytes = M * N * sizeof(float);
        double flops = 2.0 * M * N * K;

        float *h_A = (float *)malloc(M * K * sizeof(float));
        float *h_B = (float *)malloc(K * N * sizeof(float));
        float *h_C = (float *)malloc(bytes);

        fill_matrix(h_A, M * K);
        fill_matrix(h_B, K * N);

        float *d_A, *d_B, *d_C;
        cudaMalloc(&d_A, M * K * sizeof(float));
        cudaMalloc(&d_B, K * N * sizeof(float));
        cudaMalloc(&d_C, bytes);

        cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice);

        dim3 block(32, 32);
        dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);

        float naive_ms = run_kernel(launch_naive, M, N, K, d_A, d_B, d_C,
                                    grid, block, warmup, repeats);
        double naive_gflops = (flops / (naive_ms / 1000.0)) / 1e9;

        float coal_ms = run_kernel(launch_coalesced, M, N, K, d_A, d_B, d_C,
                                   grid, block, warmup, repeats);
        double coal_gflops = (flops / (coal_ms / 1000.0)) / 1e9;

        float cublas_ms = run_cublas(cublas_handle, M, N, K, d_A, d_B, d_C,
                                     warmup, repeats);
        double cublas_gflops = (flops / (cublas_ms / 1000.0)) / 1e9;

        printf("%-6d %-12.1f %-12.1f %-12.1f %-9.1f%% %-9.1f%%\n",
               M, naive_gflops, coal_gflops, cublas_gflops,
               (naive_gflops / cublas_gflops) * 100.0,
               (coal_gflops / cublas_gflops) * 100.0);

        // Validate against CPU for smallest size
        if (s == 0) {
            float *h_ref = (float *)malloc(bytes);
            cpu_gemm(M, N, K, h_A, h_B, h_ref);
            cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
            if (validate(h_ref, h_C, M * N, 1e-3f)) {
                printf("Correctness check PASSED\n");
            } else {
                printf("Correctness check FAILED\n");
            }
            free(h_ref);
        }

        free(h_A);
        free(h_B);
        free(h_C);
        cudaFree(d_A);
        cudaFree(d_B);
        cudaFree(d_C);
    }

    cublasDestroy(cublas_handle);
    return 0;
}

Here is CMakeLists.txt:

cmake_minimum_required(VERSION 3.18)
project(cuda_gemm LANGUAGES CXX CUDA)

set(CMAKE_CUDA_ARCHITECTURES 89)

find_package(CUDAToolkit REQUIRED)

add_executable(gemm src/main.cu src/kernels.cu)

target_link_libraries(gemm PRIVATE CUDA::cudart CUDA::cublas)

set_target_properties(gemm PROPERTIES
    CUDA_STANDARD 17
    CXX_STANDARD 17
)

Your project is now complete with measured evidence of your optimization. Next up, the Secret Mission challenges you to go even further with shared memory tiling.

Secret mission

Add a Shared Memory Tiled Kernel

Your coalesced kernel improved memory access patterns, but each thread still reads every element it needs directly from global DRAM. What if a block of threads could cooperatively load a chunk of data into fast on-chip shared memory, then everyone reads from that local copy? That's tiling. Implement a shared memory tiled GEMM kernel and close the gap to cuBLAS even further.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely on your local machine, so there are no ongoing cloud costs. Decide whether to keep your resources for future CUDA experiments, pause them, or delete them entirely.

Resources you used:

  • The cuda-gemm/ project folder on your Desktop (source code, build artifacts, profiling reports).
  • CUDA Toolkit 13.3 (includes nvcc, cuBLAS, and Nsight Compute).
  • CMake (build system generator).
  • GitHub repository (your cuda-gemm remote repository).

Keep everything running

No action needed. Choose this if you plan to continue experimenting with CUDA kernels or want to add more optimization passes to your GEMM project.

  • Your cuda-gemm/ folder stays on your Desktop ready for further development.
  • The CUDA Toolkit, CMake, and Nsight Compute remain installed for future GPU projects.
  • Your GitHub repository remains public for your portfolio.

Pause - I'll come back to this later

Free up disk space from build artifacts while keeping your source code intact.

  • Delete the build/ folder inside cuda-gemm/ to reclaim space from compiled binaries. Open a Command Prompt and run:
  • Keep the CUDA Toolkit and CMake installed. They take disk space but cost nothing to leave idle.
  • Keep your GitHub repository as-is for when you return.
cd %USERPROFILE%\Desktop\cuda-gemm
rmdir /s /q build

When you come back, rebuild by running cmake -B build and cmake --build build --config Release from the project folder.

Delete - I don't want to use this again

Remove all project files, installed tools, and the remote repository.

  • Delete the project folder by opening a Command Prompt and running:
cd %USERPROFILE%\Desktop
rmdir /s /q cuda-gemm
  • Uninstall the CUDA Toolkit. Open Windows Settings, select Apps, then Installed apps. Search for entries starting with "NVIDIA CUDA" and uninstall each one.
  • Uninstall CMake. In the same Installed apps list, search for "CMake", click the three-dot menu, and select Uninstall.
  • Delete the GitHub repository. Open your repository page on GitHub, click Settings, scroll down to the Danger Zone section, and click Delete this repository. Confirm by typing the repository name when prompted.

That's a Wrap!

That's a Wrap!

Nice work! You just built a GPU engineering project from scratch, writing CUDA kernels, profiling them with real tools, and optimizing based on hard data.

You've learned how to:

  • Write CUDA GEMM kernels from scratch and compile them with nvcc and CMake, benchmarking against cuBLAS to quantify the performance gap.
  • Profile kernels with Nsight Compute and identify memory coalescing as the primary bottleneck using sectors-per-request metrics.
  • Apply data-driven optimization by fixing memory access patterns, measuring the improvement with profiling evidence, and documenting the journey in a GitHub-ready repository.
  • Secret Mission: Implement shared memory tiling with __shared__ memory and __syncthreads() barriers, demonstrating understanding of the GPU memory hierarchy.

Ready to quiz yourself?