Matrix Layout & Coalesced Access

Picture a bookcase full of vectors

A matrix is a stack of vectors. Each row (or column) can be sliced out as its own vector, but the whole bookcase has to occupy one contiguous chunk of memory. The choice of which dimension is contiguous is the layout.

GPU mantra: coalesce or die. When 32 neighboring threads in a warp read 32 contiguous addresses, the hardware satisfies them in one (or close to one) memory transaction. Break that pattern and bandwidth tanks.

Layout	In-memory order	Best when
Row-major	Row elements adjacent (stride along columns is 1)	Threads stride across a row
Column-major	Column elements adjacent (stride along rows is 1)	Threads stride down a column (Fortran/cuBLAS)

Reality check on RTX 4090, 4096×4096 FP32 row-sum:

Row-major + threads scan rows (coalesced): 610 GB/s
Column-major + threads scan rows (strided): 105 GB/s

That's a 6× swing caused entirely by access pattern.

Leading dimension (ld)

The 'leading dimension' is the stride between consecutive rows (or columns). For tightly packed matrices, ld == cols for row-major or ld == rows for column-major. But you'll often pad — ld > cols — to avoid memory bank conflicts. Every cuBLAS call takes lda, ldb, ldc for exactly this reason.

Code

Row-major matrix sum (coalesced) — fast·cuda

__global__ void row_sum(const float* M, float* out, int rows, int cols, int ld) {
    int row = blockIdx.x;                  // one block per row
    int tid = threadIdx.x;
    if (row >= rows) return;

    // Threads in the warp read consecutive M[row][0..cols] elements
    // — perfectly coalesced. Each thread accumulates a partial sum,
    // then we'd reduce them in shared memory (next lesson).
    float s = 0.f;
    for (int c = tid; c < cols; c += blockDim.x) {
        s += M[row * ld + c];
    }
    // ... reduction across threads in this block (lesson 2) ...
}

Column-major matrix sum read row-wise (strided) — slow·cuda

// Same logical operation but on column-major storage:
//   M[row, col] = M_buf[col * ld + row]
// Threads in the warp now read M_buf[col*ld+row] for row in 0..warp_size,
// striding by `ld` each step → uncoalesced.
__global__ void row_sum_colmajor(const float* M_buf, float* out,
                                  int rows, int cols, int ld) {
    int row = blockIdx.x;
    int tid = threadIdx.x;
    if (row >= rows) return;
    float s = 0.f;
    for (int c = tid; c < cols; c += blockDim.x) {
        s += M_buf[c * ld + row];   // each access strides by `ld`
    }
    // Same math, ~6× slower on bandwidth-bound kernels.
}

Exercise

Take any 4096×4096 row-major matrix on a GPU you have access to. Time two kernels: (a) sum each row using one block per row, threads stride along columns; (b) sum each column using one block per column, threads stride down rows of a row-major buffer. The difference (often 5–10×) is purely the access pattern, not the math. Fix (b) by transposing the matrix first, or restructuring so each block scans its column contiguously after a transpose.