V0. Naive matrix multiplication

1. A Grid can contain multiple Blocks, and a Block can contain multiple Threads.

Grid
 ├── Block
 │    ├── Thread
 │    ├── Thread
 ├──  Block
 │    ├── ...

Grid and Blocks $\color{green}{\text{Up to 3 dimensions }}$

Grid
 ├── blockIdx.x
 ├── blockIdx.y
 └── blockIdx.z

Block
 ├── threadIdx.x
 ├── threadIdx.y
 └── threadIdx.z

- The position of each block and each thread is determined by three coordinates: (x, y, z).

2. Why is GEMM typically 2D?

Because it corresponds to the two dimensions of a matrix.
However, for 3D data—such as medical CT scans, 3D volumes, or video frame sequences—it is natural to use three dimensions.
The GEMM implementation you are looking at uses a 2D block and a 2D grid;
it is simply that the underlying CUDA architecture $\color{red}{\text{still provides three dimensions (x, y, and z).}}$

For example, a 2D block:
 dim3 block(16, 16);
is actually equivalent to:
 dim3 block(16, 16, 1);
Therefore, the thread coordinates are still:
 threadIdx.x
 threadIdx.y
 threadIdx.z
It is just that:
 threadIdx.z is always 0.

3.Define the number of threads.

dim3 block: defines the number of threads in a block. dim3 grid: defines the number of blocks in a grid.

blockDim (block Dimension) It can be understood as: "It represents the number of threads in the x, y, and z directions for the current block."
It and threadIdx form a pair of very important concepts.

4. The difference between 'threadIdx' and 'blockDim'

threadIdx: Position within the current block
blockDim: Dimensions of the current block
blockIdx: Index of the current block within the grid
gridDim: Total number of blocks in the grid

# code in the project:
int row = blockIdx.y * blockDim.y + threadIdx.y; 
int col = blockIdx.x * blockDim.x + threadIdx.x; 

5. dim3 block(16, 16); why?

dim3 block(16, 16) is widely used as a standard baseline for 2D matrix or image processing kernels because it achieves clean warp alignment without overwhelming SM resources.

6.Total number of blocks required.

dim3 grid((N + block.x - 1) / block.x, (M + block.y - 1) / block.y);
- For 512 size matrix: (512+16-1)/16 = 527/16= 32

7. Conclusion:

$\color{red}{\text{Construct a thread matrix corresponding to matrix C.}}$
Each thread corresponds to a single value in matrix C.
The calculation is performed by directly reading the corresponding row from matrix A and column from matrix B from main memory.
matrix data is linear in memory, but we can logically say it is 2D.
CUDA's 2D thread organization matches this logical structure nicely.
So 2D CUDA indexing is a programming abstraction, not a requirement imposed by the physical memory.

8. Hardware level, the GPU have things like

GPU
├── SM
│   ├── CUDA cores / execution units
│   ├── registers
│   ├── shared memory
│   └── ...
├── SM
├── SM
└── ...

9. There's an important hardware detail: warps

NVIDIA hardware executes threads in groups called warps, normally 32 threads.
dim3 block(16, 16);
So 16 × 16 = 256 threads, block becomes: 256/32=8 warps

The important point is that the hardware scheduling unit is the warp, not the 16×16 geometric picture.

$\color{purple}{\text{GPU execution → SM → Warp → Threads}}$

So the matrix is linear in memory; the 2D/3D thread layout is mainly a convenient way to express how different threads correspond to logical elements of multidimensional data.

10. CUDA Programming Structure and Workflow

This can be broadly categorized as follows:
(1) Allocating GPU memory
(2) Transferring data from CPU to GPU
(3) Invoking the CUDA kernel to perform computation
(4) Transferring data back from GPU to CPU
(5) Finally, freeing the GPU memory

11. Why is a warmup needed?

When executing a CUDA program for the first time, there may be some additional overhead, such as:
(a) CUDA runtime initialization.
(b) Overhead associated with the first kernel launch.
(c) GPU context initialization.
(d) Memory-related initialization.
(e) Caches not yet in a stable state.

12. cudaEventRecord(start)

Insert a "marker" into the execution queue of the GPU. When the GPU executes this mark in sequence, the event becomes completed; the CPU can wait for this status to know that the GPU has executed this point.
cudaEventSynchronize causes the CPU to block and wait until the GPU executes the specified Event.

Sample output (Google Colab ,T4)

# ./build/v0_naive 1024
=== V0 Naive GEMM ===
M=1024 N=1024 K=1024
  [verify] N > 512, skipping full CPU verification (too slow); suggest verifying correctness with a smaller scale separately
  [bench] avg_time=5.180 ms, GFLOPS=414.54

========= Next page ==========