Optimizing GEMM in HIP

Matrix multiplication is one of the most fundamental GPU workloads. It underlies the compute-intensive layers of deep neural networks - fully connected layers, convolutional layers expressed as implicit GEMMs, and attention mechanisms - and is central to scientific computing, computer vision, and recommendation systems. GPUs are heavily optimized for matrix multiplication, and understanding how to write an efficient GEMM kernel is an effective way to learn how GPU hardware resources interact.

This tutorial walks through seven progressive optimization steps applied to a general-purpose single-precision (FP32) matrix multiplication kernel (General Matrix Multiply, or GEMM: \(\pmb{C} = \pmb{A} \times \pmb{B}\)). Each step builds on the previous one, introducing a specific technique and explaining how to measure its effect with the ROCm performance analysis stack.

The complete source files for all steps are available at:

• Step 1 - Naive

• Step 2 - LDS tiling

• Step 3 - Register tiling

• Step 4 - Double buffering

• Step 5 - Vectorized loads

• Step 6 - Register pressure

• Step 7 - Generic kernel

Note

All examples target \(4096 \times 4096\) matrices in FP32 row-major layout and are validated against the identity-matrix test (\(B = I \implies C = A\)). They compile with amdclang++ -O3 -std=c++17 and run on any ROCm-supported GPU architecture.

Prerequisites

Before starting this tutorial, ensure the following are in place.

• ROCm installed and amdclang++ available on PATH.

• Familiarity with the HIP execution model (grids, blocks, wavefronts) and its mapping to AMD GPU hardware (dispatches, workgroups, wavefronts).

• Application tracing and profiling using rocprofv3 installed for performance analysis.

Tip

Hardware performance counter availability varies by ROCm version, Linux kernel version, and system permissions. A counter that appears in rocprofv3 --list-avail might still return only zero values on a given system. If that happens, the counter is unavailable. Kernel duration (End_Timestamp - Start_Timestamp from rocprofv3 --kernel-trace) is always available and provides a reliable baseline across all steps and architectures.

GEMM fundamentals

A GEMM multiplies an \(M \times K\) matrix \(\pmb{A}\) by a \(K \times N\) matrix \(\pmb{B}\) to produce an \(M \times N\) output matrix \(\pmb{C}\). Each element \(C_{ij}\) is the inner product of the i-th row of \(\pmb{A}\) and the j-th column of \(\pmb{B}\):

A CPU implementation applies three nested loops over m, n, and k, performing one multiply-accumulate per iteration - 2 × M × N × K scalar operations in total.

A GPU implementation is embarrassingly parallel across the output elements. A HIP kernel assigns one thread (or a small tile of threads) to each output element, eliminating the m and n loops entirely and leaving only the k reduction loop inside each thread. With \(M \times N\) output elements and a modern AMD GPU fielding tens of thousands of concurrent threads, the full output matrix can be computed in a single dispatch - provided data can be supplied fast enough to keep the Compute Units (CUs) busy.

Background: the GEMM arithmetic intensity

For an \(M \times K \times N\) GEMM the arithmetic intensity - floating-point operations per byte of DRAM traffic - is:

\[I = \frac{2 \cdot M \cdot N \cdot K}{\text{sizeof}(\text{float}) \cdot (M \cdot K + K \cdot N + M \cdot N)}\]

For \(M = N = K = n\) this simplifies to \(\frac{n}{6}\). With \(n = 4096\) that gives roughly 683 FLOPs/byte, far above the roofline ridge point of any current AMD GPU. GEMM is therefore compute-bound in principle - but only if data is supplied fast enough to keep the compute units busy. The naive kernel falls well below the roofline because it is memory-bound in practice: global memory latency stalls dominate (see Roofline model for background on roofline analysis and performance bottlenecks).

The optimization steps that follow progressively close the gap between actual and theoretical throughput by improving data reuse and instruction-level efficiency.

Step 1: Naive kernel

The naive kernel assigns one thread per output element. Each thread reads a full row of \(\pmb{A}\) and a full column of \(\pmb{B}\) directly from global memory.

While this maps naturally onto the GPU’s parallel execution model, it produces severe cache thrashing. Consider that any element \(A[i][k]\) is needed by all N threads that compute a different output column in the same row, and \(B[k][j]\) is needed by all M threads that compute a different output row in the same column. With tens of thousands of threads in flight simultaneously, the working set far exceeds the L2 cache, so data that should be reused is evicted before the next thread requests it. The result is that total DRAM traffic is a large multiple of the minimum required bandwidth (sizeof(float) * (M*K + K*N + M*N)), and the kernel is firmly memory-bound despite GEMM’s high arithmetic intensity in principle.

__global__ void matrix_multiply_naive(const float* __restrict__ A,
                                      const float* __restrict__ B,
                                      float* __restrict__ C,
                                      int m,
                                      int n,
                                      int k)
{
    const int col = blockIdx.x * blockDim.x + threadIdx.x;
    const int row = blockIdx.y * blockDim.y + threadIdx.y;

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

Launch configuration:

const dim3 block(BLOCK_SIZE, BLOCK_SIZE);
const dim3 grid((N + BLOCK_SIZE - 1) / BLOCK_SIZE,
                (M + BLOCK_SIZE - 1) / BLOCK_SIZE);

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_naive.hip -o mm_naive
./mm_naive

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_naive

Kernel duration is End_Timestamp - Start_Timestamp (both in nanoseconds). The trace also captures VGPR_Count and Scratch_Size for each dispatch.

The --kernel-trace CSV is sufficient to establish the baseline: the naive kernel’s End_Timestamp - Start_Timestamp will be the slowest of all steps because every global memory access is a cache miss.

Step 2: LDS tiling

The root cause of the naive kernel’s cache thrashing is that all threads share a single transparent L2 cache with no way to guarantee that a loaded value stays resident until every thread that needs it has read it. AMD GPUs expose Local Data Share (LDS) - a low-latency, high-bandwidth on-chip memory that is explicitly managed by the programmer, functioning as a programmable L1 cache. Unlike CPU hardware caches, data placed in LDS stays there until the kernel explicitly overwrites or discards it.

The key insight is that every element of \(\pmb{A}\) is used by \(N\) threads (one per output column) and every element of \(\pmb{B}\) is used by \(M\) threads. Caching a TILE_SIZE * TILE_SIZE strip of \(\pmb{A}\) and \(\pmb{B}\) in LDS lets all TILE_SIZE² threads in a block reuse that data without touching global memory again.

constexpr int TILE_SIZE = 16;

Shared memory allocation:

// Allocate two LDS tiles: one strip of A (row-major) and one of B (row-major).
// Each tile holds TILE_SIZE x TILE_SIZE float elements = 1 KiB.
__shared__ float tile_a[TILE_SIZE][TILE_SIZE];
__shared__ float tile_b[TILE_SIZE][TILE_SIZE];

Load phase (cooperative, one element per thread):

// Cooperative load: thread (tx, ty) loads one element of each tile.
//
// tile_a ← A[row][t*TILE_SIZE + tx]
//   Row comes from the block's output row; column strides across K.
const int a_col = t * TILE_SIZE + tx;
tile_a[ty][tx]  = (row < m && a_col < k) ? A[row * k + a_col] : 0.0f;

// tile_b ← B[t*TILE_SIZE + ty][col]
//   Row strides across K; column comes from the block's output column.
const int b_row = t * TILE_SIZE + ty;
tile_b[ty][tx]  = (b_row < k && col < n) ? B[b_row * n + col] : 0.0f;

// Ensure all threads have written their element before any thread reads.
__syncthreads();

Compute phase (inner product from LDS):

// Accumulate the partial dot product for this tile from LDS.
// Both tile_a and tile_b reside in LDS; no global memory is touched here.
for(int i = 0; i < TILE_SIZE; ++i)
{
    sum += tile_a[ty][i] * tile_b[i][tx];
}

// Ensure all threads have finished reading before the next iteration
// overwrites the tiles.
__syncthreads();

LDS bank conflict analysis

Moving data into LDS is only half the battle - how threads access that data determines whether the LDS delivers its full bandwidth. LDS is divided into independently addressable banks. When threads in the same cycle access different addresses that map to the same bank, the hardware must serialize those accesses. This is called a bank conflict, and it directly reduces the effective LDS bandwidth by the degree of the conflict (a k-way conflict takes k cycles instead of one).

Bank mapping is straightforward: consecutive 4-byte words are assigned to consecutive banks in round-robin order. For a 32-bank LDS, word at byte address a maps to bank (a / 4) % 32. Two threads accessing the same address are not in conflict - the hardware broadcasts the value to both.

The number of LDS banks varies across AMD GPU architectures:

Architecture	Banks	Entries per bank (4 bytes each)
CDNA	32	512
CDNA2	32	512
CDNA3	32	512
CDNA4	64	640
RDNA2	64	512
RDNA3	64	512
RDNA3.5	64	512
RDNA4	64	512

Note

On RDNA GPUs the 64 banks are sub-divided into two sets of 32 banks, each affiliated with a pair of SIMD32 units within the WGP. A wavefront executes on one SIMD32 and maps its accesses to the affiliated 32-bank set. In CU mode (-mcumode), a workgroup runs on a single CU and accesses one set of 32 banks; in WGP mode (the default) the workgroup spans the full WGP.

For this kernel’s compute phase - an inner product over the K-strip:

sum += tile_a[ty][i] * tile_b[i][tx];

• the access pattern with float data is **inherently conflict-free

regardless of tile size**:

• tile_a[ty][i]: all threads in a wavefront that share the same ty read the same address. The hardware broadcasts the value - no conflict.

• tile_b[i][tx]: each thread has a unique tx, and because each float is exactly 4 bytes (= one bank slot), consecutive tx values always map to consecutive banks. No two threads in the same cycle can hit the same bank at a different address, regardless of how large TILE_SIZE is. The i loop iterates sequentially within each thread, so accesses to different rows of tile_b are never concurrent.

In other words, with FP32 data and the simple inner-product pattern of this step, bank conflicts are a non-issue. The bank mechanism is worth understanding now, however, because it becomes a real concern later.

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_lds.hip -o mm_lds
./mm_lds

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_lds

Compare End_Timestamp - Start_Timestamp against the Step 1 baseline. To measure L2-to-High Bandwidth Memory (HBM) read traffic on CDNA GPUs:

# CDNA only
rocprofv3 --pmc TCP_TCC_READ_REQ_sum --output-format csv -- ./mm_lds

What to observe after this step

Compare the following counters against the Step 1 naive kernel baseline.

Counter	What to look for
Kernel duration (kernel-trace CSV)	End_Timestamp - Start_Timestamp should drop significantly versus the naive kernel, confirming that LDS data reuse is reducing global memory traffic.
TCP_TCC_READ_REQ_sum (CDNA)	CDNA only: reduction proportional to TILE_SIZE confirms fewer L2-to-HBM read requests.
LDSBankConflict	Should remain at or near zero for TILE_SIZE=16 with FP32 data (confirms no bank conflicts).

Step 3: Register tiling

In the LDS kernel each thread computes exactly one output element, reading TILE_SIZE values from tile_a and TILE_SIZE values from tile_b for every K-strip. If instead each thread computes a THREAD_TILE_M × THREAD_TILE_N sub-tile in registers, it amortizes the LDS load cost across THREAD_TILE_M × THREAD_TILE_N outputs.

The outer product of a length-THREAD_TILE_M column fragment of A and a length-THREAD_TILE_N row fragment of B produces a full THREAD_TILE_M × THREAD_TILE_N block of C contributions using only THREAD_TILE_M + THREAD_TILE_N LDS reads instead of THREAD_TILE_M × THREAD_TILE_N.

Tile parameters:

constexpr int BLOCK_TILE_M  = 128; // Rows of C computed by one block
constexpr int BLOCK_TILE_N  = 128; // Columns of C computed by one block
constexpr int K_TILE_SIZE   = 16;  // K-strip width loaded into LDS each iteration
constexpr int THREAD_TILE_M = 8;   // Rows of C computed by one thread
constexpr int THREAD_TILE_N = 8;   // Columns of C computed by one thread

Transposed tile_b layout:

tile_b is stored transposed in LDS as tile_b_T[BLOCK_TILE_N][K_TILE_SIZE] (column-major for B). This layout choice is motivated by two concerns:

1. Stride-1 reads during the outer-product compute phase: In the outer-product pattern, each thread loads a fragment of THREAD_TILE_N consecutive elements from tile_b_T along the ki dimension. Because the inner dimension of tile_b_T is K_TILE_SIZE, consecutive ki values are adjacent in memory - stride-1 access. Without transposition (tile_b[K_TILE_SIZE][BLOCK_TILE_N]), the fragment load would stride across the large BLOCK_TILE_N dimension, producing scattered LDS reads.

2. Bank-conflict safety for sub-4-byte types: In Step 2’s simple inner-product loop, each float occupies exactly one 4-byte LDS bank slot, so bank conflicts cannot arise regardless of tile dimensions. This property does not hold for smaller data types. When architecture-specific builtins introduce half-precision (FP16, 2 bytes) or quarter-precision (FP8, 1 byte) data in follow-up sections, multiple elements pack into a single 4-byte bank slot. If the row stride of tile_b equals or is a multiple of the bank count, different threads can address different sub-word elements within the same bank slot, producing conflicts that are impossible with FP32. The transposed layout decouples the fragment access stride (K_TILE_SIZE) from the tile’s outer dimension (BLOCK_TILE_N), avoiding this class of conflict regardless of element size or bank count.

LDS allocation:

// tile_a: row-major strip of A, shape [BLOCK_TILE_M][K_TILE_SIZE]
__shared__ float tile_a[BLOCK_TILE_M][K_TILE_SIZE];
// tile_b_T: column-major strip of B (transposed in LDS), shape [BLOCK_TILE_N][K_TILE_SIZE]
// Storing transposed eliminates the bank-conflict risk that arises when
// TILE_SIZE equals the LDS bank count, and produces stride-1 reads during
// both the load phase and the compute phase.
__shared__ float tile_b_T[BLOCK_TILE_N][K_TILE_SIZE];

Cooperative tile load:

// --- Cooperative load of tile_a (BLOCK_TILE_M x K_TILE_SIZE) ---
// Each thread loads one or more elements by striding through the flat
// tile index space with step = num_threads.
{
    const int tile_elems = BLOCK_TILE_M * K_TILE_SIZE;
    for(int idx = tid; idx < tile_elems; idx += num_threads)
    {
        const int tile_row  = idx / K_TILE_SIZE;
        const int tile_col  = idx % K_TILE_SIZE;
        const int global_row = block_row_start + tile_row;
        const int global_col = t * K_TILE_SIZE + tile_col;
        tile_a[tile_row][tile_col] =
            (global_row < m && global_col < k)
                ? A[global_row * k + global_col]
                : 0.0f;
    }
}

// --- Cooperative load of tile_b_T (BLOCK_TILE_N x K_TILE_SIZE, transposed) ---
// We load B[b_row][b_col] into tile_b_T[b_col - block_col_start][b_row - t*K].
// Equivalently: thread loads the element at (col_within_tile, k_within_tile)
// of the transposed tile, reading B in row-major order.
{
    const int tile_elems = BLOCK_TILE_N * K_TILE_SIZE;
    for(int idx = tid; idx < tile_elems; idx += num_threads)
    {
        const int tile_col_idx = idx / K_TILE_SIZE; // column within block tile
        const int tile_k_idx   = idx % K_TILE_SIZE; // k index within strip
        const int global_row   = t * K_TILE_SIZE + tile_k_idx;
        const int global_col   = block_col_start + tile_col_idx;
        tile_b_T[tile_col_idx][tile_k_idx] =
            (global_row < k && global_col < n)
                ? B[global_row * n + global_col]
                : 0.0f;
    }
}

__syncthreads();

Outer-product accumulation:

// Outer-product accumulation entirely in registers.
// For each k-index within the strip:
//   - Load a column fragment of tile_a into registers (a_frag).
//   - Load a row fragment of tile_b_T into registers (b_frag).
//   - Accumulate the TM x TN outer product into acc.
for(int ki = 0; ki < K_TILE_SIZE; ++ki)
{
    float a_frag[THREAD_TILE_M];
    #pragma unroll
    for(int i = 0; i < THREAD_TILE_M; ++i)
    {
        a_frag[i] = tile_a[thread_row + i][ki];
    }

    float b_frag[THREAD_TILE_N];
    #pragma unroll
    for(int j = 0; j < THREAD_TILE_N; ++j)
    {
        b_frag[j] = tile_b_T[thread_col + j][ki];
    }

    #pragma unroll
    for(int i = 0; i < THREAD_TILE_M; ++i)
    {
        #pragma unroll
        for(int j = 0; j < THREAD_TILE_N; ++j)
        {
            acc[i][j] += a_frag[i] * b_frag[j];
        }
    }
}

__syncthreads();

Write-back:

// Write the register accumulator to global memory.
// Each store is a coalesced write: threads with consecutive tx values write
// to consecutive columns (thread_col strides by THREAD_TILE_N = 8, but the
// inner j loop produces 8 consecutive output columns per thread, so adjacent
// threads in x cover adjacent 8-wide strips → effectively coalesced).
#pragma unroll
for(int i = 0; i < THREAD_TILE_M; ++i)
{
    #pragma unroll
    for(int j = 0; j < THREAD_TILE_N; ++j)
    {
        const int out_row = block_row_start + thread_row + i;
        const int out_col = block_col_start + thread_col + j;
        if(out_row < m && out_col < n)
        {
            C[out_row * n + out_col] = acc[i][j];
        }
    }
}

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_register_tiling.hip -o mm_register_tiling
./mm_register_tiling

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_register_tiling

For LDS and arithmetic instruction counts:

rocprofv3 --pmc SQ_INSTS_LDS --output-format csv -- ./mm_register_tiling

What to observe

Compare the following counters against the LDS tiling kernel from Step 2.

Counter	What to look for
VALUInsts (RDNA, CDNA, CDNA2) / SQ_INSTS_VALU (CDNA3, CDNA4)	Increase in VALU instructions per wave (more Fused Multiply-Accumulate (FMA) operations per LDS read).
SQ_INSTS_LDS	Reduction in LDS instructions per wave (THREAD_TILE_M + THREAD_TILE_N instead of 2 × THREAD_TILE_M × THREAD_TILE_N).
SQ_WAIT_INST_LDS	Reduction in LDS stall cycles (register reuse hides LDS latency).

Tile parameter tuning guidance

The tile dimensions (BLOCK_TILE_M, BLOCK_TILE_N, K_TILE_SIZE, THREAD_TILE_M, THREAD_TILE_N) are not one-size-fits-all. Optimal values depend on several interacting constraints:

Matrix dimensions

Tile sizes should evenly divide the matrix dimensions to avoid boundary handling overhead. Padding the matrices to a multiple of the tile size is common in production GEMM libraries.

LDS capacity

Each workgroup allocates BLOCK_TILE_M × K_TILE_SIZE and K_TILE_SIZE × BLOCK_TILE_N floats in LDS. LDS capacity varies across AMD GPU architectures (typically 64–160 KiB per compute unit). Exceeding the LDS budget reduces occupancy by limiting how many workgroups can be resident simultaneously.

Register file

Each thread holds a THREAD_TILE_M × THREAD_TILE_N accumulator array plus fragment temporaries. Larger thread tiles increase arithmetic intensity but consume more Vector General-Purpose Registers (VGPRs), reducing occupancy. Step 6 addresses this tradeoff directly.

Architecture-to-architecture variation

LDS bank counts, wavefront widths, VGPR file size, and L1/L2 cache line sizes differ across CDNA and RDNA families. Tile sizes that are optimal on a CDNA3 GPU might not be optimal on an RDNA4 GPU. Use rocprofv3 --pmc to measure LDS efficiency, VGPR usage, and occupancy on each target, and re-tune accordingly.

Step 4: Double buffering

Every iteration of the K-strip loop stalls at __syncthreads() waiting for LDS tile loads to complete before compute can begin. Software double buffering hides this latency by maintaining two pairs of LDS buffers - one pair for A tiles and one for B tiles - each with a ping and a pong slot. The next tile is loaded into the background slot while the current slot is being consumed.

Both buffering strategies are hidden behind a TilePolicy interface so that the kernel body is identical regardless of the chosen approach.

Policy interface overview:

Method	Responsibility
prologue	Load tile 0 into buffer 0 and synchronize (double-buffer only; no-op for single).
prefetch	Issue the load for the next tile into the background buffer.
acquire	Synchronize before compute (single-buffer: __syncthreads(); double: no-op).
release	Synchronize after compute (both: __syncthreads()).
buf_idx	Return which buffer to read for the current iteration.

Single-buffer policy (baseline - same logic as Step 3):

struct SingleBufferTilePolicy
{
    static constexpr int num_buffers = 1;

    // LDS holds one pair of tiles.
    struct SharedStorage
    {
        float tile_a[BLOCK_TILE_M][K_TILE_SIZE];
        float tile_b_T[BLOCK_TILE_N][K_TILE_SIZE];
    };

    // No preloading: the first prefetch() call loads tile 0.
    static __device__ void prologue(SharedStorage& /*smem*/,
                                    const float* /*A*/,
                                    const float* /*B*/,
                                    int /*block_row_start*/,
                                    int /*block_col_start*/,
                                    int /*m*/,
                                    int /*n*/,
                                    int /*k*/,
                                    int /*tid*/,
                                    int /*num_threads*/,
                                    int /*num_tiles*/)
    {}

    // Load tile `tile_idx - 1` (the CURRENT tile t) into the single buffer.
    // The kernel calls prefetch(t + 1), so tile_idx = t + 1.  We load
    // tile t = tile_idx - 1 here; acquire() then __syncthreads() to ensure
    // it is visible before compute begins.
    static __device__ void prefetch(SharedStorage&  smem,
                                    const float*    A,
                                    const float*    B,
                                    int             tile_idx,
                                    int             block_row_start,
                                    int             block_col_start,
                                    int             m,
                                    int             n,
                                    int             k,
                                    int             tid,
                                    int             num_threads,
                                    int             /*num_tiles*/)
    {
        // tile_idx = t + 1, so the tile to load is t = tile_idx - 1.
        // The loop runs t in [0, num_tiles), so tile_idx in [1, num_tiles],
        // and load_idx in [0, num_tiles - 1] - always in range.
        const int load_idx = tile_idx - 1;

        const int tile_elems_a = BLOCK_TILE_M * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_a; idx += num_threads)
        {
            const int tile_row   = idx / K_TILE_SIZE;
            const int tile_col   = idx % K_TILE_SIZE;
            const int global_row = block_row_start + tile_row;
            const int global_col = load_idx * K_TILE_SIZE + tile_col;
            smem.tile_a[tile_row][tile_col] =
                (global_row < m && global_col < k)
                    ? A[global_row * k + global_col]
                    : 0.0f;
        }

        const int tile_elems_b = BLOCK_TILE_N * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_b; idx += num_threads)
        {
            const int tile_col_idx = idx / K_TILE_SIZE;
            const int tile_k_idx   = idx % K_TILE_SIZE;
            const int global_row   = load_idx * K_TILE_SIZE + tile_k_idx;
            const int global_col   = block_col_start + tile_col_idx;
            smem.tile_b_T[tile_col_idx][tile_k_idx] =
                (global_row < k && global_col < n)
                    ? B[global_row * n + global_col]
                    : 0.0f;
        }
    }

    // Wait for the prefetch just issued to complete.
    static __device__ void acquire(int /*tile_idx*/)
    {
        __syncthreads();
    }

    // Protect the buffer: no thread may start the next prefetch until all
    // threads have finished reading the current tile.
    static __device__ void release(int /*tile_idx*/)
    {
        __syncthreads();
    }

    // Always return buffer index 0.
    static __device__ int buf_idx(int /*tile_idx*/)
    {
        return 0;
    }

    static __device__ const float* get_tile_a(const SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_a[0][0];
    }

    static __device__ const float* get_tile_b_T(const SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_b_T[0][0];
    }
};

Software double-buffer policy:

struct SoftwareDoubleBufferTilePolicy
{
    static constexpr int num_buffers = 2;

    // LDS holds TWO pairs of tiles, indexed by [buf][...].
    struct SharedStorage
    {
        float tile_a[2][BLOCK_TILE_M][K_TILE_SIZE];
        float tile_b_T[2][BLOCK_TILE_N][K_TILE_SIZE];
    };

    // Load tile 0 into buf[0] and synchronise so that the first iteration
    // can call acquire() as a no-op.
    static __device__ void prologue(SharedStorage&  smem,
                                    const float*    A,
                                    const float*    B,
                                    int             block_row_start,
                                    int             block_col_start,
                                    int             m,
                                    int             n,
                                    int             k,
                                    int             tid,
                                    int             num_threads,
                                    int             /*num_tiles*/)
    {
        const int tile_elems_a = BLOCK_TILE_M * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_a; idx += num_threads)
        {
            const int tile_row   = idx / K_TILE_SIZE;
            const int tile_col   = idx % K_TILE_SIZE;
            const int global_row = block_row_start + tile_row;
            const int global_col = tile_col; // tile 0: K offset = 0
            smem.tile_a[0][tile_row][tile_col] =
                (global_row < m && global_col < k)
                    ? A[global_row * k + global_col]
                    : 0.0f;
        }

        const int tile_elems_b = BLOCK_TILE_N * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_b; idx += num_threads)
        {
            const int tile_col_idx = idx / K_TILE_SIZE;
            const int tile_k_idx   = idx % K_TILE_SIZE;
            const int global_row   = tile_k_idx; // tile 0: K offset = 0
            const int global_col   = block_col_start + tile_col_idx;
            smem.tile_b_T[0][tile_col_idx][tile_k_idx] =
                (global_row < k && global_col < n)
                    ? B[global_row * n + global_col]
                    : 0.0f;
        }

        // Ensure tile 0 is fully in LDS before the main loop begins.
        __syncthreads();
    }

    // Issue the load for tile `tile_idx` into the alternate buffer.
    // This load runs concurrently with the compute on the current buffer.
    static __device__ void prefetch(SharedStorage&  smem,
                                    const float*    A,
                                    const float*    B,
                                    int             tile_idx,
                                    int             block_row_start,
                                    int             block_col_start,
                                    int             m,
                                    int             n,
                                    int             k,
                                    int             tid,
                                    int             num_threads,
                                    int             num_tiles)
    {
        if(tile_idx >= num_tiles)
        {
            return;
        }

        // Write into the buffer that is NOT currently being read.
        const int dst_buf = tile_idx % 2;

        const int tile_elems_a = BLOCK_TILE_M * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_a; idx += num_threads)
        {
            const int tile_row   = idx / K_TILE_SIZE;
            const int tile_col   = idx % K_TILE_SIZE;
            const int global_row = block_row_start + tile_row;
            const int global_col = tile_idx * K_TILE_SIZE + tile_col;
            smem.tile_a[dst_buf][tile_row][tile_col] =
                (global_row < m && global_col < k)
                    ? A[global_row * k + global_col]
                    : 0.0f;
        }

        const int tile_elems_b = BLOCK_TILE_N * K_TILE_SIZE;
        for(int idx = tid; idx < tile_elems_b; idx += num_threads)
        {
            const int tile_col_idx = idx / K_TILE_SIZE;
            const int tile_k_idx   = idx % K_TILE_SIZE;
            const int global_row   = tile_idx * K_TILE_SIZE + tile_k_idx;
            const int global_col   = block_col_start + tile_col_idx;
            smem.tile_b_T[dst_buf][tile_col_idx][tile_k_idx] =
                (global_row < k && global_col < n)
                    ? B[global_row * n + global_col]
                    : 0.0f;
        }
        // No sync here: the load runs concurrently with compute on the other buffer.
    }

    // The current buffer is already filled (by prologue or the previous release).
    // No synchronisation is needed here.
    static __device__ void acquire(int /*tile_idx*/) {}

    // Single __syncthreads() serves dual purpose:
    //   - Ensures every thread has finished reading buf[tile_idx % 2]
    //     (so the next prefetch may safely overwrite it).
    //   - Ensures the prefetch issued at the top of this iteration
    //     (loading into buf[(tile_idx+1) % 2]) has completed.
    static __device__ void release(int /*tile_idx*/)
    {
        __syncthreads();
    }

    static __device__ int buf_idx(int tile_idx)
    {
        return tile_idx % 2;
    }

    static __device__ const float* get_tile_a(const SharedStorage& smem, int buf)
    {
        return &smem.tile_a[buf][0][0];
    }

    static __device__ const float* get_tile_b_T(const SharedStorage& smem, int buf)
    {
        return &smem.tile_b_T[buf][0][0];
    }
};

Compile-time policy validation:

template<typename TilePolicy>
struct TilePolicyTraits
{
    static_assert(TilePolicy::num_buffers == 1 || TilePolicy::num_buffers == 2,
                  "TilePolicy::num_buffers must be 1 (single buffer) or 2 (double buffer)");

    // SharedStorage must be a non-empty, trivially-destructible aggregate
    // that the compiler can place in __shared__ memory.
    static_assert(sizeof(typename TilePolicy::SharedStorage) > 0,
                  "TilePolicy::SharedStorage must have non-zero size");
    static_assert(std::is_trivially_destructible_v<typename TilePolicy::SharedStorage>,
                  "TilePolicy::SharedStorage must be trivially destructible "
                  "(required for __shared__ variables)");
};

Unified kernel template:

template<typename TilePolicy>
__global__ void matrix_multiply_buffered(const float* __restrict__ A,
                                         const float* __restrict__ B,
                                         float* __restrict__ C,
                                         int m,
                                         int n,
                                         int k)
{
    // Instantiate the trait checker; all static_asserts fire at compile time.
    (void)TilePolicyTraits<TilePolicy>{};

    __shared__ typename TilePolicy::SharedStorage smem;

    const int tx          = threadIdx.x;
    const int ty          = threadIdx.y;
    const int tid         = ty * blockDim.x + tx;
    const int num_threads = blockDim.x * blockDim.y;

    const int thread_row      = ty * THREAD_TILE_M;
    const int thread_col      = tx * THREAD_TILE_N;
    const int block_row_start = blockIdx.y * BLOCK_TILE_M;
    const int block_col_start = blockIdx.x * BLOCK_TILE_N;

    float acc[THREAD_TILE_M][THREAD_TILE_N] = {};

    const int num_tiles = (k + K_TILE_SIZE - 1) / K_TILE_SIZE;

    // [Sphinx prologue start]
    // Prologue: SingleBuffer → no-op; DoubleBuffer → preload tile 0 + sync.
    TilePolicy::prologue(smem, A, B, block_row_start, block_col_start,
                         m, n, k, tid, num_threads, num_tiles);
    // [Sphinx prologue end]

    for(int t = 0; t < num_tiles; ++t)
    {
        // [Sphinx prefetch acquire start]
        // Issue the load for the NEXT tile into the alternate buffer (or the
        // single buffer for SingleBufferTilePolicy).
        TilePolicy::prefetch(smem, A, B, t + 1, block_row_start, block_col_start,
                             m, n, k, tid, num_threads, num_tiles);

        // Ensure the CURRENT tile (tile t) is fully in LDS before reading it.
        // SingleBuffer: __syncthreads() after the load just issued above.
        // DoubleBuffer: no-op (tile t was already in LDS from the previous release).
        TilePolicy::acquire(t);
        // [Sphinx prefetch acquire end]

        // Outer-product accumulation from LDS (register tiling, same as Step 3).
        const float* ta  = TilePolicy::get_tile_a(smem, TilePolicy::buf_idx(t));
        const float* tb_T = TilePolicy::get_tile_b_T(smem, TilePolicy::buf_idx(t));

        for(int ki = 0; ki < K_TILE_SIZE; ++ki)
        {
            float a_frag[THREAD_TILE_M];
            #pragma unroll
            for(int i = 0; i < THREAD_TILE_M; ++i)
            {
                // tile_a is laid out as [BLOCK_TILE_M][K_TILE_SIZE]
                a_frag[i] = ta[(thread_row + i) * K_TILE_SIZE + ki];
            }

            float b_frag[THREAD_TILE_N];
            #pragma unroll
            for(int j = 0; j < THREAD_TILE_N; ++j)
            {
                // tile_b_T is laid out as [BLOCK_TILE_N][K_TILE_SIZE]
                b_frag[j] = tb_T[(thread_col + j) * K_TILE_SIZE + ki];
            }

            #pragma unroll
            for(int i = 0; i < THREAD_TILE_M; ++i)
            {
                #pragma unroll
                for(int j = 0; j < THREAD_TILE_N; ++j)
                {
                    acc[i][j] += a_frag[i] * b_frag[j];
                }
            }
        }

        // [Sphinx release start]
        // SingleBuffer: second __syncthreads() – protects LDS before prefetch
        //               overwrites it in the next iteration.
        // DoubleBuffer: single __syncthreads() – ensures:
        //   (a) compute has finished reading buf[t%2], so the next prefetch
        //       may safely write into it.
        //   (b) the prefetch issued above (for tile t+1 into buf[(t+1)%2])
        //       has completed before the next iteration's acquire() returns.
        TilePolicy::release(t);
        // [Sphinx release end]
    }

    // Write accumulator to global C
    #pragma unroll
    for(int i = 0; i < THREAD_TILE_M; ++i)
    {
        #pragma unroll
        for(int j = 0; j < THREAD_TILE_N; ++j)
        {
            const int out_row = block_row_start + thread_row + i;
            const int out_col = block_col_start + thread_col + j;
            if(out_row < m && out_col < n)
            {
                C[out_row * n + out_col] = acc[i][j];
            }
        }
    }
}

Note

The LDS footprint doubles with software double buffering: two copies of the A and B tile buffers are needed instead of one. For the parameters in this example (BLOCK_TILE_M = BLOCK_TILE_N = 128, K_TILE_SIZE = 16):

• Single-buffer LDS: 128 × 16 × 4 × 2 = 16 KiB

• Double-buffer LDS: 16 KiB × 2 = 32 KiB

This is within the 64–160 KiB LDS budget on all supported architectures, but leaves less headroom for occupancy. Use rocprofv3 --pmc MeanOccupancyPerCU to verify occupancy does not drop when switching from single- to double-buffered policy.

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_double_buffer.hip -o mm_double_buffer
./mm_double_buffer

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_double_buffer

To measure LDS stall cycles (RDNA3 and all CDNA):

# RDNA3 and all CDNA
rocprofv3 --pmc SQ_WAIT_INST_LDS --output-format csv -- ./mm_double_buffer

What to observe

Compare the following counters against the register tiling kernel from Step 3.

Counter	What to look for
SQ_WAIT_INST_LDS	Reduction in LDS stall cycles (load latency hidden by prefetch).
Kernel duration (kernel-trace CSV)	Should remain similar to Step 3 (prefetch hides latency but does not reduce total data fetched).
TCP_TCC_READ_REQ_sum	CDNA only: roughly constant vs Step 3 (same number of L2-to-HBM reads; only latency is hidden, not traffic).

Step 5: Vectorized loads

Each global memory load instruction in the tile-loading loop fetches one float per thread. Replacing it with a float2 or float4 load fetches 2 or 4 float values per instruction - the same total data moves through the cache hierarchy, but in fewer instructions. This reduces pressure on the VMEM instruction-issue pipeline and can improve overall throughput when instruction issue is the bottleneck rather than memory bandwidth.

To put this in context, consider how a single coalesced scalar load of a float maps to cache lines on each architecture family:

Architecture	Wavefront width	L1/L0 cache line size	Cache lines per coalesced scalar load
CDNA	64 lanes	64 bytes	64 × 4 B / 64 B = 4
CDNA2	64 lanes	64 bytes	64 × 4 B / 64 B = 4
CDNA3	64 lanes	128 bytes	64 × 4 B / 128 B = 2
CDNA4	64 lanes	128 bytes	64 × 4 B / 128 B = 2
RDNA2	32 lanes	128 bytes	32 × 4 B / 128 B = 1
RDNA3	32 lanes	128 bytes	32 × 4 B / 128 B = 1
RDNA3.5	32 lanes	128 bytes	32 × 4 B / 128 B = 1
RDNA4	32 lanes	128 bytes	32 × 4 B / 128 B = 1

On RDNA GPUs, a coalesced scalar float load already fills exactly one cache line - wider vector loads do not reduce cache traffic. On CDNA and CDNA2 a scalar load spans 4 cache lines. On CDNA3 and CDNA4, it spans 2. In all cases, vector loads do not change the number of cache lines accessed; they reduce the number of instructions the wavefront must issue to move the same amount of data.

Two vector widths are shown alongside the scalar baseline:

Width	HIP type	Alignment	Instruction
1	float	4 bytes	buffer_load_dword (one 4 B element per thread)
2	float2	8 bytes	buffer_load_dwordx2 (two 4 B elements per thread)
4	float4	16 bytes	buffer_load_dwordx4 (four 4 B elements per thread)

Vector type helper:

template<int Width>
struct VectorType;

template<>
struct VectorType<1>
{
    using type                      = float;
    static constexpr int width      = 1;
    static constexpr int byte_align = alignof(float); // 4 bytes
};

template<>
struct VectorType<2>
{
    using type                      = float2;
    static constexpr int width      = 2;
    static constexpr int byte_align = 8; // float2 must be 8-byte aligned
};

template<>
struct VectorType<4>
{
    using type                      = float4;
    static constexpr int width      = 4;
    static constexpr int byte_align = 16; // float4 must be 16-byte aligned
};

Vectorized load function:

template<int Width>
__device__ void vectorized_load(float* __restrict__       dst,
                                const float* __restrict__ src,
                                int                       tile_elems,
                                int                       tid,
                                int                       num_threads)
{
    using Vec = typename VectorType<Width>::type;

    static_assert(sizeof(Vec) == Width * sizeof(float),
                  "Vector type size mismatch");

    // Number of complete vector elements in the tile.
    // The static_asserts in the kernel guarantee tile_elems % Width == 0.
    const int vec_elems = tile_elems / Width;

    const Vec* src_vec = reinterpret_cast<const Vec*>(src);
    Vec*       dst_vec = reinterpret_cast<Vec*>(dst);

    for(int idx = tid; idx < vec_elems; idx += num_threads)
    {
        dst_vec[idx] = src_vec[idx];
    }
}

Alignment requirements:

The reinterpret_cast in vectorized_load is only safe when the source pointer is aligned to sizeof(Vec) bytes. Two alignment guarantees apply:

1. hipMalloc returns a pointer aligned to at least 256 bytes, satisfying float4 (16 bytes) for the base of any matrix.

2. Each row of A begins at offset row × K × sizeof(float). For float4 loads the row length K must be a multiple of 4. A compile-time static_assert enforces this for the tile parameters.

A runtime check is included in the example to catch misaligned user-supplied pointers:

void check_alignment(const void* ptr, int required_bytes, const char* name)
{
    const uintptr_t addr = reinterpret_cast<uintptr_t>(ptr);
    if(addr % static_cast<uintptr_t>(required_bytes) != 0)
    {
        std::cerr << "WARNING: " << name << " pointer " << ptr
                  << " is not " << required_bytes << "-byte aligned. "
                  << "Vectorised loads of width " << required_bytes
                  << " bytes will produce undefined behaviour.\n";
    }
}

Note

Only the A-tile load is vectorized because its elements are contiguous in global memory (row-major, stride 1). The B-tile is loaded column-by-column (stride N in global memory), which is not amenable to simple vector loads. Architecture-specific builtins (MFMA and WMMA), covered in follow-up sections, address this asymmetry.

Vectorized loads and smaller data types

For FP32, vectorized loads are a moderate optimization: they reduce instruction count, but a scalar load already fills cache lines well (see table above). The picture changes significantly for smaller data types. With FP16 (2 bytes) or FP8 (1 byte), a scalar load per thread no longer fills a full cache line:

Element type	Size	RDNA scalar load (32 lanes)	Cache line fill (128 B line)
FP32	4 bytes	32 × 4 = 128 B	100% (1 full line)
FP16	2 bytes	32 × 2 = 64 B	50% (half a line wasted)
FP8	1 byte	32 × 1 = 32 B	25% (three quarters wasted)

The wasted portion of each cache line is fetched from DRAM but never used - this is pure bandwidth overhead. A 2-wide vector load for FP16 or a 4-wide vector load for FP8 restores the full cache line fill. On CDNA3/CDNA4 (128 B cache line, 64-wide wavefronts), FP8 scalar loads similarly fill only half a line.

This is why the TilePolicy parameterizes the vector load width: the optimal width depends on both the element type and the target architecture. For the FP32 case in this tutorial, the benefit is modest, but for the low-precision builtins introduced in follow-up sections, vectorized loads become essential to avoid wasting memory bandwidth.

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_vectorized.hip -o mm_vectorized
./mm_vectorized

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_vectorized

To compare VMEM instruction cycles across scalar and vector variants:

# RDNA (combined VMEM counter)
rocprofv3 --pmc SQ_INST_CYCLES_VMEM --output-format csv -- ./mm_vectorized

# CDNA (read counter only; GEMM tile loads are reads)
rocprofv3 --pmc SQ_INST_CYCLES_VMEM_RD --output-format csv -- ./mm_vectorized

What to observe

Compare the following counters against the double-buffered kernel from Step 4.

Counter	What to look for
SQ_INST_CYCLES_VMEM (all RDNA GPUs) / SQ_INST_CYCLES_VMEM_RD (all CDNA)	Reduction in VMEM instruction cycles (fewer instructions for the same total data). CDNA GPUs split this into SQ_INST_CYCLES_VMEM_RD (reads) and SQ_INST_CYCLES_VMEM_WR (writes); use the read counter for GEMM tile loads.
TCP_TOTAL_CACHE_ACCESSES	CDNA only: total L2 cache accesses should remain roughly constant (cache-line traffic does not change – only the instruction count does).

Step 6: Register pressure and occupancy

GPU occupancy - the ratio of active wavefronts to the hardware maximum - is set by the most constrained resource. For register-tiled GEMM kernels that resource is typically the VGPR file: each thread holds THREAD_TILE_M × THREAD_TILE_N accumulator registers plus fragment arrays, and the compiler might allocate additional temporaries.

Higher register usage means fewer concurrent wavefronts per CU, which reduces the GPU’s ability to hide memory latency through wavefront switching. Conversely, forcibly reducing register usage can increase register spilling to scratch memory (a slow VGPR-to-DRAM path), degrading performance.

Three kernel variants illustrate the tradeoff:

No annotation (compiler decides freely):

__global__ void matrix_multiply_no_hint(const float* __restrict__ A,
                                        const float* __restrict__ B,
                                        float* __restrict__ C,
                                        int m,
                                        int n,
                                        int k)
{
    __shared__ float tile_a[BLOCK_TILE_M][K_TILE_SIZE];
    __shared__ float tile_b_T[BLOCK_TILE_N][K_TILE_SIZE];
    gemm_body(A, B, C, m, n, k, tile_a, tile_b_T);
}

``__launch_bounds__``:

constexpr int MIN_WAVES_PER_EU = 2; // tune this with rocprofv3

__global__
__launch_bounds__(BLOCK_THREADS, MIN_WAVES_PER_EU)
void matrix_multiply_launch_bounds(const float* __restrict__ A,
                                   const float* __restrict__ B,
                                   float* __restrict__ C,
                                   int m,
                                   int n,
                                   int k)
{
    __shared__ float tile_a[BLOCK_TILE_M][K_TILE_SIZE];
    __shared__ float tile_b_T[BLOCK_TILE_N][K_TILE_SIZE];
    gemm_body(A, B, C, m, n, k, tile_a, tile_b_T);
}

The two-argument form __launch_bounds__(max_threads, min_waves_per_eu) tells the compiler to allocate VGPRs such that at least min_waves_per_eu wavefronts can be resident per EU simultaneously, given max_threads / warpSize wavefronts per block.

``[[clang::amdgpu_waves_per_eu]]``:

// tune WAVES_PER_EU_MIN and WAVES_PER_EU_MAX with rocprofv3 --kernel-trace
// (compare VGPR_Count across variants).
constexpr int WAVES_PER_EU_MIN = 2;
constexpr int WAVES_PER_EU_MAX = 4;

__global__ void matrix_multiply_waves_per_eu
[[clang::amdgpu_waves_per_eu(WAVES_PER_EU_MIN, WAVES_PER_EU_MAX)]]
(const float* __restrict__ A,
 const float* __restrict__ B,
 float* __restrict__ C,
 int m,
 int n,
 int k)
{
    __shared__ float tile_a[BLOCK_TILE_M][K_TILE_SIZE];
    __shared__ float tile_b_T[BLOCK_TILE_N][K_TILE_SIZE];
    gemm_body(A, B, C, m, n, k, tile_a, tile_b_T);
}

This Clang attribute directly instructs the backend to target a wavefront occupancy in the range [min, max] per EU.

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_launch_bounds.hip -o mm_launch_bounds
./mm_launch_bounds

Finding optimal values with rocprofv3

1. Collect static kernel metadata for all three variants:

   rocprofv3 --kernel-trace --output-format csv -- ./mm_launch_bounds
   

   In the resulting CSV, compare the VGPR_Count and Scratch_Size columns across the three kernels.

2. Collect dynamic occupancy counters:

   rocprofv3 --pmc SQ_WAVES_sum MeanOccupancyPerCU SQ_WAIT_INST_LDS --output-format csv -- ./mm_launch_bounds
   
   # if SQ_WAIT_INST_LDS not available
   rocprofv3 --pmc SQ_WAVES_sum MeanOccupancyPerCU --output-format csv -- ./mm_launch_bounds
   
   # if MeanOccupancyPerCU and SQ_WAVES_sum not available
   rocprofv3 --pmc SQ_LEVEL_WAVES SQ_WAIT_INST_LDS --output-format csv -- ./mm_launch_bounds
   

3. Open both CSVs and note the values of:

   • VGPR_Count (vector registers allocated per thread, from the kernel-trace CSV)

   • Scratch_Size (>0 means VGPRs are spilling to DRAM - avoid this)

   • MeanOccupancyPerCU or SQ_LEVEL_WAVES: mean active wavefronts per CU, from the PMC CSV

4. Calculate the theoretical maximum occupancy from VGPR_Count using the formula for your architecture (available in the AMD ISA reference).

5. If the compiler allocated more VGPRs than necessary and occupancy is below the target, add __launch_bounds__ with a min_waves_per_eu that reflects the desired occupancy.

6. Re-profile and re-check Scratch_Size - if it increases significantly, the compiler was forced to spill and the constraint is too aggressive.

Note

The numeric values MIN_WAVES_PER_EU = 2 and WAVES_PER_EU_MAX = 4 in the example are illustrative. Optimal values depend on the target GPU and the exact kernel register usage shown in the rocprofv3 --kernel-trace CSV. Always re-profile with rocprofv3 after applying the annotation to confirm that occupancy improves without introducing scratch-memory spilling.

What to observe

Compare the following counters across the three kernel variants (no annotation, __launch_bounds__, and [[clang::amdgpu_waves_per_eu]]).

Counter / column	What to look for
VGPR_Count (kernel-trace CSV)	From --kernel-trace CSV: should decrease after adding the annotation.
Scratch_Size (kernel-trace CSV)	From --kernel-trace CSV: must remain zero (non-zero means register spilling to DRAM).
MeanOccupancyPerCU + SQ_WAVES_sum / SQ_LEVEL_WAVES	From --pmc CSV: mean active wavefronts per CU; should rise with fewer VGPRs.
SQ_WAIT_INST_LDS	LDS stall cycles (fall when more wavefronts are resident).

Step 7: Generic kernel

By this point the kernel is already a strong general-purpose implementation: it uses LDS tiling to exploit data reuse, register tiling to maximize arithmetic intensity, software double buffering to hide load latency, and vectorized loads to reduce memory transaction overhead. For many workloads this is sufficient.

Squeezing out the last few percent of throughput, however, requires architecture-specific matrix-multiply instructions: Matrix Fused Multiply-Accumulate (MFMA) on CDNA GPUs and Wave Matrix Multiply-Accumulate (WMMA) on RDNA3 and RDNA4. These instructions perform a small matrix multiply directly in hardware and deliver substantially higher FLOP/s than an equivalent sequence of scalar FMAs.

Steps 1–6 produced a well-optimized scalar GEMM kernel, but repeating the same work for each architecture-specific instruction set - MFMA on CDNA, WMMA on RDNA3, the relaxed WMMA variant on RDNA4 - would mean maintaining several near-identical copies of the kernel with only the inner computation swapped out. Any future improvement (a new tiling strategy, a wider vector load, a different buffering depth) would have to be applied to every copy independently.

The goal of this step is to factor the kernel so that the optimized orchestration is written once and architecture-specific builtins can be dropped in later as a policy, without touching any kernel code.

The insight from the preceding steps is that all the work decomposes into exactly two independent concerns:

1. Data movement (TilePolicy) - tile shape, LDS layout, vector load width, and buffering strategy. These optimizations are identical regardless of which instruction is used to compute the output; a float4 load into a double-buffered LDS tile is just as beneficial whether the inner loop uses scalar FMAs or MFMA.

2. Arithmetic (ComputePolicy) - how a thread’s register fragment is loaded from LDS and how the output accumulator is updated. This is the only part that differs between scalar code and architecture-specific builtins.

Following the principle of lifting an algorithm into its most general form, these two responsibilities are encapsulated in two orthogonal policy classes: TilePolicy and ComputePolicy. A single kernel template matrix_multiply_generic<TilePolicy, ComputePolicy> orchestrates the common control flow while delegating all architecture-specific details to the policies.

The kernel presented here uses ScalarFMAPolicy as the ComputePolicy. It already incorporates all the optimizations from the preceding steps: LDS tiling, register tiling, software double buffering, and vectorized loads. When an MFMA or WMMA ComputePolicy is provided in one of the architecture-specific builtins chapters, the same data-movement infrastructure and the same kernel orchestration are reused unchanged - only the inner arithmetic changes.

Policy interfaces

Both interfaces are documented in full in matrix_multiply_generic.hip. The key design decisions are:

TilePolicy interface summary:

Requirement	Rationale
num_buffers = 1 | 2	Governs LDS layout (one or two buffer pairs) and sync strategy.
block_tile_m, block_tile_n, k_tile_size	Tile shape needed by the kernel to compute grid dimensions.
SharedStorage (trivially destructible)	Placed in __shared__; must not require a destructor call.
prologue, prefetch, acquire, release	Data movement hooks; see below.

ComputePolicy interface summary:

Requirement	Rationale
thread_tile_m, thread_tile_n	Output sub-tile per thread; kernel derives block dimensions from these.
effective_lanes	Unique output lanes per wavefront.
thread_tile_offset, tid, lane_id, *row, *col	Architecture-aware thread-to-output mapping.
elem_a, elem_b	Element types of the register fragments (for example, float or __half); the kernel loop is fully templated on these.
k_step	Number of k-indices consumed per mma() call (1 for scalar FMA; higher values for builtins that process multiple k-indices per call); the kernel loop advances ki by this amount.
load_a, load_b, mma, store_c	Fragment load, multiply-accumulate, and write-back.

Data type scope: policy coverage

elem_a and elem_b parameterize the register fragment type and are already fully wired through the kernel loop. A future ComputePolicy can set elem_a = __half and the float-to-half conversion happens entirely inside load_a / load_b - the kernel body is untouched.

However, two things are not yet parameterized and are hardcoded to float in this file:

• The LDS storage type - SharedStorage in every TilePolicy holds float arrays, and the cooperative load helpers write float into LDS.

• The global memory pointer type - the kernel signature takes const float* A, const float* B, float* C.

This means two distinct cases arise when introducing builtins in follow-up sections:

Case	Example	What needs to change
FP32 in memory, low-precision fragments	Global float → LDS float → __half fragment for MFMA	Only ComputePolicy::load_a and load_b (conversion in registers). TilePolicy and the kernel signature are unchanged.
Low-precision in memory	Global __half → LDS __half → __half fragment	TilePolicy needs an InputT template parameter so SharedStorage and the cooperative load helpers use InputT instead of float. The kernel signature changes from const float* to const InputT*.

The global memory and LDS remain float throughout this tutorial, so only Case 1 applies. The DirectLoadTilePolicy below is a Case 1 example at the TilePolicy level; a ComputePolicy that converts to FP16 in its load_a / load_b would be a Case 1 example at the compute level. Case 2 is left as an extension for architecture-specific follow-up sections that operate on native FP16 or FP8 input matrices.

Compile-time validation

Both policy interfaces are validated with C++17 static_assert traits:

TilePolicy traits:

template<typename TilePolicy>
struct TilePolicyTraits
{
    // num_buffers must be 1 or 2
    static_assert(TilePolicy::num_buffers == 1 || TilePolicy::num_buffers == 2,
                  "TilePolicy::num_buffers must be 1 (single) or 2 (double-buffer)");

    // Tile dimensions must be positive multiples of k_tile_size
    static_assert(TilePolicy::block_tile_m > 0,
                  "TilePolicy::block_tile_m must be a positive integer");
    static_assert(TilePolicy::block_tile_n > 0,
                  "TilePolicy::block_tile_n must be a positive integer");
    static_assert(TilePolicy::k_tile_size > 0,
                  "TilePolicy::k_tile_size must be a positive integer");

    // SharedStorage must be trivially destructible (no LDS destructor call)
    static_assert(std::is_trivially_destructible_v<typename TilePolicy::SharedStorage>,
                  "TilePolicy::SharedStorage must be trivially destructible");

    // SharedStorage must actually occupy LDS
    static_assert(sizeof(typename TilePolicy::SharedStorage) > 0,
                  "TilePolicy::SharedStorage must have non-zero size");
};

ComputePolicy traits:

template<typename ComputePolicy>
struct ComputePolicyTraits
{
    static_assert(ComputePolicy::thread_tile_m > 0,
                  "ComputePolicy::thread_tile_m must be a positive integer");
    static_assert(ComputePolicy::thread_tile_n > 0,
                  "ComputePolicy::thread_tile_n must be a positive integer");
    static_assert(ComputePolicy::effective_lanes > 0,
                  "ComputePolicy::effective_lanes must be a positive integer");
    static_assert(ComputePolicy::k_step > 0,
                  "ComputePolicy::k_step must be a positive integer");

    // Accumulator must be trivially constructible (no runtime init in registers)
    static_assert(std::is_trivially_constructible_v<typename ComputePolicy::Accumulator>,
                  "ComputePolicy::Accumulator must be trivially constructible");

    // Fragment element types must have a well-defined size (covers both
    // scalar arithmetic types like float and vector types like half2).
    static_assert(sizeof(typename ComputePolicy::elem_a) > 0,
                  "ComputePolicy::elem_a must have non-zero size");
    static_assert(sizeof(typename ComputePolicy::elem_b) > 0,
                  "ComputePolicy::elem_b must have non-zero size");
};

C++17 versus C++20: concept syntax

The static_assert approach above is portable to C++17 but requires the traits struct to be instantiated explicitly, and error messages appear at the trait instantiation site. C++20 requires clauses provide a more ergonomic alternative:

template<typename TilePolicy>
concept TilePolicyConcept = requires
{
    requires TilePolicy::num_buffers == 1 || TilePolicy::num_buffers == 2;
    requires TilePolicy::block_tile_m > 0;
    requires TilePolicy::block_tile_n > 0;
    requires TilePolicy::k_tile_size > 0;
    requires std::is_trivially_destructible_v<typename TilePolicy::SharedStorage>;
    requires sizeof(typename TilePolicy::SharedStorage) > 0;
};

// The kernel declaration then becomes self-documenting:
template<TilePolicyConcept TilePolicy, typename ComputePolicy>
__global__ void matrix_multiply_generic(/* … */);
//
// Constraint violations are reported at the point of the template
// instantiation with a clear "constraint not satisfied" diagnostic rather than
// a cryptic static_assert failure inside a trait struct.

With C++20, the constraint is expressed directly in the function signature and violations are reported at the point of the invalid template instantiation with a clear “constraint not satisfied” diagnostic.

Concrete policies

ScalarFMAPolicy - portable scalar FP32 outer-product (no builtins):

template<int ThreadTileM_, int ThreadTileN_, int BlockDimX_>
struct ScalarFMAPolicy
{
    static constexpr int thread_tile_m  = ThreadTileM_;
    static constexpr int thread_tile_n  = ThreadTileN_;
    // All lanes in a wavefront compute unique output elements (no mirroring).
    static constexpr int effective_lanes = 64; // AMD wavefront size
    static constexpr int k_step         = 1;  // one scalar FMA per ki iteration

    using elem_a = float;
    using elem_b = float;

    struct Accumulator
    {
        float data[ThreadTileM_][ThreadTileN_];
    };

    __device__ static void zero(Accumulator& acc)
    {
        #pragma unroll
        for(int i = 0; i < ThreadTileM_; ++i)
            #pragma unroll
            for(int j = 0; j < ThreadTileN_; ++j)
                acc.data[i][j] = 0.0f;
    }

    // Map thread id to the top-left corner of its output sub-tile.
    // lane_id is accepted but unused (interface forward-compatibility).
    __device__ static void thread_tile_offset(int tid,
                                              int /*lane_id*/,
                                              int* thread_row,
                                              int* thread_col)
    {
        const int tx = tid % BlockDimX_;
        const int ty = tid / BlockDimX_;
        *thread_row  = ty * ThreadTileM_;
        *thread_col  = tx * ThreadTileN_;
    }

    // Load a column fragment of tile_a into registers.
    // tile_a is stored row-major: tile_a[row][ki] = ptr[row * k_tile_size + ki]
    __device__ static void load_a(const float* tile_a_ptr,
                                  int          tile_a_row,
                                  int          ki,
                                  int          k_tile_size,
                                  int          /*lane_id*/,
                                  elem_a (&frag)[ThreadTileM_])
    {
        #pragma unroll
        for(int i = 0; i < ThreadTileM_; ++i)
            frag[i] = tile_a_ptr[(tile_a_row + i) * k_tile_size + ki];
    }

    // Load a row fragment of tile_b_T into registers.
    // tile_b_T is stored column-major: tile_b_T[col][ki] = ptr[col * k_tile_size + ki]
    __device__ static void load_b(const float* tile_b_T_ptr,
                                  int          tile_b_col,
                                  int          ki,
                                  int          k_tile_size,
                                  int          /*lane_id*/,
                                  elem_b (&frag)[ThreadTileN_])
    {
        #pragma unroll
        for(int j = 0; j < ThreadTileN_; ++j)
            frag[j] = tile_b_T_ptr[(tile_b_col + j) * k_tile_size + ki];
    }

    // Outer-product FMA: acc += a_frag ⊗ b_frag
    __device__ static void mma(Accumulator&       acc,
                               const elem_a (&a_frag)[ThreadTileM_],
                               const elem_b (&b_frag)[ThreadTileN_])
    {
        #pragma unroll
        for(int i = 0; i < ThreadTileM_; ++i)
            #pragma unroll
            for(int j = 0; j < ThreadTileN_; ++j)
                acc.data[i][j] += a_frag[i] * b_frag[j];
    }

    // Write accumulator to global memory.
    // lane_id accepted but unused; all lanes write their unique sub-tile.
    __device__ static void store_c(const Accumulator& acc,
                                   float*             C,
                                   int                out_row_base,
                                   int                out_col_base,
                                   int                m,
                                   int                n,
                                   int                /*lane_id*/)
    {
        #pragma unroll
        for(int i = 0; i < ThreadTileM_; ++i)
            #pragma unroll
            for(int j = 0; j < ThreadTileN_; ++j)
            {
                const int r = out_row_base + i;
                const int c = out_col_base + j;
                if(r < m && c < n)
                    C[r * n + c] = acc.data[i][j];
            }
    }
};

SingleBufferTilePolicy - single LDS buffer pair, equivalent to Step 3:

template<int BlockTileM_, int BlockTileN_, int KTileSize_>
struct SingleBufferTilePolicy
{
    static constexpr int num_buffers  = 1;
    static constexpr int block_tile_m = BlockTileM_;
    static constexpr int block_tile_n = BlockTileN_;
    static constexpr int k_tile_size  = KTileSize_;

    struct SharedStorage
    {
        float tile_a[BlockTileM_][KTileSize_];
        float tile_b_T[BlockTileN_][KTileSize_];
    };

    __device__ static void prologue(SharedStorage& /*smem*/,
                                    const float* /*A*/, const float* /*B*/,
                                    int /*block_row_start*/, int /*block_col_start*/,
                                    int /*m*/, int /*n*/, int /*k*/,
                                    int /*tid*/, int /*num_threads*/)
    {
        // No prologue needed for single-buffer; tile is loaded at the top of
        // each loop iteration in prefetch().
    }

    __device__ static void prefetch(SharedStorage& smem,
                                    const float* A, const float* B,
                                    int          t_next,
                                    int block_row_start, int block_col_start,
                                    int m, int n, int k,
                                    int tid, int num_threads)
    {
        // For single-buffer we load the *current* tile (t_next is actually t
        // from the caller's perspective - see kernel loop).
        // Rename for clarity: t_next here is the index we want to load.
        const int t = t_next;
        cooperative_load_tile_a(&smem.tile_a[0][0], A, block_row_start,
                                block_tile_m, t, k_tile_size, m, k, tid, num_threads);
        cooperative_load_tile_b_T(&smem.tile_b_T[0][0], B, block_col_start,
                                   block_tile_n, t, k_tile_size, k, n, tid, num_threads);
    }

    // After prefetch: wait for the load to complete before compute.
    __device__ static void acquire(int /*t*/) { __syncthreads(); }

    // After compute: wait for all threads to finish reading before the next
    // iteration overwrites the tile.
    __device__ static void release(int /*t*/) { __syncthreads(); }

    __device__ static int buf_idx(int /*t*/) { return 0; }

    __device__ static float* get_tile_a(SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_a[0][0];
    }
    __device__ static float* get_tile_b_T(SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_b_T[0][0];
    }
};

SoftwareDoubleBufferTilePolicy - ping-pong LDS buffers, equivalent to Step 4:

template<int BlockTileM_, int BlockTileN_, int KTileSize_>
struct SoftwareDoubleBufferTilePolicy
{
    static constexpr int num_buffers  = 2;
    static constexpr int block_tile_m = BlockTileM_;
    static constexpr int block_tile_n = BlockTileN_;
    static constexpr int k_tile_size  = KTileSize_;

    struct SharedStorage
    {
        float tile_a[2][BlockTileM_][KTileSize_];
        float tile_b_T[2][BlockTileN_][KTileSize_];
    };

    // Prologue: load tile 0 into buffer 0 and synchronise so compute can
    // begin immediately at the start of the first iteration.
    __device__ static void prologue(SharedStorage& smem,
                                    const float* A, const float* B,
                                    int block_row_start, int block_col_start,
                                    int m, int n, int k,
                                    int tid, int num_threads)
    {
        cooperative_load_tile_a(&smem.tile_a[0][0][0], A, block_row_start,
                                block_tile_m, 0, k_tile_size, m, k, tid, num_threads);
        cooperative_load_tile_b_T(&smem.tile_b_T[0][0][0], B, block_col_start,
                                   block_tile_n, 0, k_tile_size, k, n, tid, num_threads);
        __syncthreads(); // tile 0 visible to all threads before loop entry
    }

    // Prefetch tile t_next into the background buffer (1 - current) while
    // the current buffer is in use for compute.  No sync here - the sync that
    // ends the previous compute phase (release) also covers this load because
    // the background buffer is not touched until two iterations later.
    __device__ static void prefetch(SharedStorage& smem,
                                    const float* A, const float* B,
                                    int          t_next,
                                    int block_row_start, int block_col_start,
                                    int m, int n, int k,
                                    int tid, int num_threads)
    {
        if(t_next * k_tile_size < k) // guard: do not load beyond K
        {
            const int buf = t_next % 2; // background buffer index
            cooperative_load_tile_a(&smem.tile_a[buf][0][0], A, block_row_start,
                                    block_tile_m, t_next, k_tile_size, m, k, tid, num_threads);
            cooperative_load_tile_b_T(&smem.tile_b_T[buf][0][0], B, block_col_start,
                                       block_tile_n, t_next, k_tile_size, k, n, tid, num_threads);
        }
    }

    // acquire: the prologue already placed a barrier; nothing needed here.
    // The previous iteration's release() synced both the completed compute
    // and the in-flight prefetch into the background buffer.
    __device__ static void acquire(int /*t*/) {}

    // release: one __syncthreads covers two things:
    //   (a) all threads finished computing from the current buffer, and
    //   (b) all threads finished writing the prefetched tile into the
    //       background buffer (the load issued in prefetch() above).
    // Both must be visible before the next iteration swaps buffers.
    __device__ static void release(int /*t*/) { __syncthreads(); }

    __device__ static int buf_idx(int t) { return t % 2; }

    __device__ static float* get_tile_a(SharedStorage& smem, int buf)
    {
        return &smem.tile_a[buf][0][0];
    }
    __device__ static float* get_tile_b_T(SharedStorage& smem, int buf)
    {
        return &smem.tile_b_T[buf][0][0];
    }
};

Generic kernel template

template<typename TilePolicy, typename ComputePolicy>
__global__ void
    __launch_bounds__(TilePolicy::block_tile_m / ComputePolicy::thread_tile_m
                          * TilePolicy::block_tile_n / ComputePolicy::thread_tile_n)
    matrix_multiply_generic(const float* __restrict__ A,
                            const float* __restrict__ B,
                            float* __restrict__       C,
                            int                       m,
                            int                       n,
                            int                       k)
{
    // Instantiate TilePolicyTraits/ComputePolicyTraits to trigger static_assert
    // checks at the point where this kernel is instantiated.
    (void)TilePolicyTraits<TilePolicy>{};
    (void)ComputePolicyTraits<ComputePolicy>{};

    // [Sphinx gemm kernel smem start]
    __shared__ typename TilePolicy::SharedStorage smem;
    // [Sphinx gemm kernel smem end]

    const int tid         = threadIdx.y * blockDim.x + threadIdx.x;
    const int lane_id     = tid % 64; // wavefront lane index (AMD: 64-wide)
    const int num_threads = blockDim.x * blockDim.y;

    // Top-left corner of the block tile in global C
    const int block_row_start = blockIdx.y * TilePolicy::block_tile_m;
    const int block_col_start = blockIdx.x * TilePolicy::block_tile_n;

    // Thread's output sub-tile offset within the block tile
    int thread_row = 0, thread_col = 0;
    ComputePolicy::thread_tile_offset(tid, lane_id, &thread_row, &thread_col);

    // Initialise accumulator to zero
    typename ComputePolicy::Accumulator acc;
    ComputePolicy::zero(acc);

    const int num_tiles = (k + TilePolicy::k_tile_size - 1) / TilePolicy::k_tile_size;

    // [Sphinx gemm kernel prologue start]
    // Double-buffer prologue: load tile 0 for policies that support it.
    // For SingleBufferTilePolicy this is a no-op.
    TilePolicy::prologue(smem, A, B,
                         block_row_start, block_col_start,
                         m, n, k, tid, num_threads);
    // [Sphinx gemm kernel prologue end]

    // [Sphinx gemm kernel loop start]
    for(int t = 0; t < num_tiles; ++t)
    {
        // --- Prefetch next tile (or load current tile for single-buffer) ---
        // For SoftwareDoubleBufferTilePolicy: issues the load for tile (t+1)
        //   into the background buffer while we compute from the foreground.
        // For SingleBufferTilePolicy: loads tile t (t_next == t here).
        TilePolicy::prefetch(smem, A, B, t + (TilePolicy::num_buffers == 1 ? 0 : 1),
                             block_row_start, block_col_start,
                             m, n, k, tid, num_threads);

        // --- Acquire: synchronise before compute ---
        // Single-buffer:  __syncthreads() after load
        // Double-buffer:  no-op (previous release() already synced)
        TilePolicy::acquire(t);

        // --- Compute: outer-product accumulation from the active buffer ---
        const int buf = TilePolicy::buf_idx(t);
        const float* tile_a_ptr   = TilePolicy::get_tile_a(smem, buf);
        const float* tile_b_T_ptr = TilePolicy::get_tile_b_T(smem, buf);

        // k_step: number of k-indices consumed per mma() call.
        // ScalarFMAPolicy uses k_step=1 (one FMA per call); intrinsic-based
        // policies can use higher values (e.g. k_step=2 for dot-product ops).
        static_assert(TilePolicy::k_tile_size % ComputePolicy::k_step == 0,
                      "k_tile_size must be divisible by ComputePolicy::k_step");

        for(int ki = 0; ki < TilePolicy::k_tile_size; ki += ComputePolicy::k_step)
        {
            typename ComputePolicy::elem_a a_frag[ComputePolicy::thread_tile_m];
            typename ComputePolicy::elem_b b_frag[ComputePolicy::thread_tile_n];

            ComputePolicy::load_a(tile_a_ptr, thread_row, ki,
                                   TilePolicy::k_tile_size, lane_id, a_frag);
            ComputePolicy::load_b(tile_b_T_ptr, thread_col, ki,
                                   TilePolicy::k_tile_size, lane_id, b_frag);
            ComputePolicy::mma(acc, a_frag, b_frag);
        }

        // --- Release: synchronise after compute ---
        // Single-buffer:  __syncthreads() before next iteration overwrites tile
        // Double-buffer:  __syncthreads() covers both compute-done and prefetch-done
        TilePolicy::release(t);
    }
    // [Sphinx gemm kernel loop end]

    // [Sphinx gemm kernel store start]
    // Write the accumulator to global memory.
    ComputePolicy::store_c(acc, C,
                            block_row_start + thread_row,
                            block_col_start + thread_col,
                            m, n, lane_id);
    // [Sphinx gemm kernel store end]
}

Architecture dispatch

A compile-time dispatch block selects the appropriate ComputePolicy based on the target Instruction Set Architecture (ISA). The architecture-specific MFMA and WMMA policies are stubs to be filled by follow-up architecture-specific sections:

//
// Selecting a ComputePolicy at compile time based on device ISA:
//
//   gfx90a / gfx942   →  MFMAPolicy<float,float>   (CDNA2/3)
//   gfx1100 / gfx1101                   →  WMMAPolicy<__half,float>  (RDNA3)
//   gfx1200 / gfx1201                   →  RDNA4WMMAPolicy (relaxed lane rules)
//   fallback (all others)               →  ScalarFMAPolicy
//
// Usage pattern (not yet compiled; shown for documentation):
//
// #if defined(__gfx90a__) || defined(__gfx942__)
//   using DefaultComputePolicy = MFMAPolicy</* … */>;
// #elif defined(__gfx1100__) || defined(__gfx1101__) || defined(__gfx1102__)
//   using DefaultComputePolicy = WMMAPolicy</* … */>;
// #elif defined(__gfx1200__) || defined(__gfx1201__)
//   using DefaultComputePolicy = RDNA4WMMAPolicy</* … */>;
// #else
//   using DefaultComputePolicy = ScalarFMAPolicy<8, 8, 16>;
// #endif
//
// The kernel template (GemmKernel) is identical for all policies; only the
// type alias changes.

Policy aliases used in this example:

// Single-buffer variant (baseline – equivalent to register-tiled step 3)
using SingleBufPolicy = SingleBufferTilePolicy<BLOCK_TILE_M, BLOCK_TILE_N, K_TILE_SIZE>;

// Double-buffer variant (software ping-pong, same tile/compute parameters)
using DoubleBufPolicy = SoftwareDoubleBufferTilePolicy<BLOCK_TILE_M, BLOCK_TILE_N, K_TILE_SIZE>;

// Portable scalar outer-product compute policy
using ScalarPolicy = ScalarFMAPolicy<THREAD_TILE_M, THREAD_TILE_N, BLOCK_DIM_X>;

// Direct global→LDS tile policy (CDNA3/CDNA4, bypasses VGPRs during tile load)
#if defined(__gfx942__) || defined(__gfx950__)
using DirectLoadPolicy = DirectLoadTilePolicy<BLOCK_TILE_M, BLOCK_TILE_N, K_TILE_SIZE>;
#endif

Compile and run:

amdclang++ -O3 -std=c++17 matrix_multiply_generic.hip -o mm_generic
./mm_generic

The program launches three variants:

1. GemmKernel<SingleBufPolicy, ScalarPolicy> - single-buffer, scalar FP32

2. GemmKernel<DoubleBufPolicy, ScalarPolicy> - double-buffer, scalar FP32

3. GemmKernel<DirectLoadPolicy, ScalarPolicy> - direct global-to-LDS loads (CDNA3 and CDNA4 only)

The first two differ only in their TilePolicy; the third replaces the standard cooperative load with a hardware-specific builtin. All three share the same kernel template and ComputePolicy.

Profile wall-clock time with rocprofv3:

rocprofv3 --kernel-trace --output-format csv -- ./mm_generic

The --kernel-trace CSV shows VGPR_Count and LDS_Block_Size for each of the three dispatched kernels. Compare these columns across the SingleBufPolicy, DoubleBufPolicy, and DirectLoadPolicy variants to confirm the expected differences in register and LDS usage.

Builtin drop-in: DirectLoadTilePolicy

The DirectLoadTilePolicy demonstrates how an architecture-specific builtin can be used as a drop-in TilePolicy without touching the kernel template or the ComputePolicy.

On CDNA3 and CDNA4, the __builtin_amdgcn_global_load_lds intrinsic transfers data from global memory directly into LDS without staging in vector registers:

// Gather 64 floats from global memory into 64 contiguous LDS locations.
// Each lane provides its own global source address (per-lane VADDR).
// The hardware writes lane k's value to dst_chunk + k * sizeof(float).
__builtin_amdgcn_global_load_lds(
    src_lane,      // per-lane global address (VADDR)
    dst_chunk,     // wave-uniform LDS base (-> M0)
    4,             // size per lane in bytes (immediate)
    0,             // offset (immediate)
    0);            // cache policy (immediate)

There are three distinct benefits:

1. Instruction count reduction - each global_load_lds_dword replaces a global_load_dword into a VGPR followed by a ds_write_b32 from that VGPR into LDS. A single instruction does the work of two, halving the total instruction count for the tile-load phase. For the tile parameters in this example (BLOCK_TILE_M = 128, K_TILE_SIZE = 16, 4 wavefronts) each wavefront issues 32 global_load_lds_dword instructions instead of 32 global_load_dword + 32 ds_write_b32 = 64 instructions. This directly frees instruction-issue bandwidth for the FMA compute phase.

2. VGPR file bandwidth - the VGPR file is no longer used as a staging area for tile data during the load phase. Its read/write bandwidth is fully available to the outer-product FMA loop, reducing contention between the load and compute phases.

3. VGPR count - the loaded data never occupies vector registers. Fewer VGPRs allocated means more wavefronts can be resident per CU simultaneously (see Step 6), which improves the hardware’s ability to hide memory latency through wavefront switching.

On throughput-sensitive workloads the first effect dominates: if the scalar tile-load sequence is instruction-issue-bound, halving its instruction count can produce speedups larger than VGPR savings alone would suggest.

At the ISA level, global_load_lds_dword is a wavefront-wide gather:

• VADDR (a vector register) provides each lane’s global source address - the global pointer need not be wave-uniform.

• M0 (a scalar register) provides the LDS base address - wave-uniform.

• The LDS write destination is implicitly offset per lane: lane k writes to M0 + offset + k * 4 (for size <= 4), or M0 + offset + k * 16 (for size > 4).

A single global_load_lds_dword instruction therefore gathers 64 floats (256 bytes) from 64 potentially different global addresses into 64 contiguous LDS locations - all without touching VGPRs.

The DirectLoadTilePolicy treats the tile as a flat array of elements and processes it in chunks of 64 (one wavefront width per instruction). For BLOCK_TILE_M = 128, K_TILE_SIZE = 16: 2048 elements / 64 per instruction = 32 chunks. With 4 wavefronts: 8 instructions per wavefront to fill the entire tile.

template<int BlockTileM_, int BlockTileN_, int KTileSize_>
struct DirectLoadTilePolicy
{
    static constexpr int num_buffers  = 1;
    static constexpr int block_tile_m = BlockTileM_;
    static constexpr int block_tile_n = BlockTileN_;
    static constexpr int k_tile_size  = KTileSize_;

    static_assert((BlockTileM_ * KTileSize_) % 64 == 0,
                  "tile_a element count must be a multiple of 64 (wavefront width)");
    static_assert((BlockTileN_ * KTileSize_) % 64 == 0,
                  "tile_b_T element count must be a multiple of 64 (wavefront width)");

    struct SharedStorage
    {
        float tile_a[BlockTileM_][KTileSize_];
        float tile_b_T[BlockTileN_][KTileSize_];
    };

    // Helper: gather-load a flat tile from global memory to LDS.
    //
    // Elements are processed in chunks of 64 (one per wavefront-wide
    // instruction).  Each lane computes its own global source address;
    // the hardware writes lane k's value to M0 + k * sizeof(float).
    //
    // Parameters:
    //   lds_dest       – base of the flat LDS tile array
    //   global_src     – base of the global matrix
    //   num_rows/cols  – tile dimensions (LDS layout is [num_rows][num_cols])
    //   global_row/col – top-left corner of this tile in the global matrix
    //   global_stride  – row stride of the global matrix
    //   transpose      – if true, swap row/col in the global read pattern
    //                    (used for tile_b_T: LDS[col][ki], global B[ki][col])
    __device__ static void direct_load_flat(float*       lds_dest,
                                            const float* global_src,
                                            int          num_rows,
                                            int          num_cols,
                                            int          global_row_start,
                                            int          global_col_start,
                                            int          global_stride,
                                            bool         transpose,
                                            int          tid,
                                            int          num_threads)
    {
        constexpr int warp_size    = 64;
        const int     wave_id      = tid / warp_size;
        const int     lane         = tid % warp_size;
        const int     num_waves    = num_threads / warp_size;
        const int     total        = num_rows * num_cols;
        const int     total_chunks = total / warp_size;

        for(int chunk = wave_id; chunk < total_chunks; chunk += num_waves)
        {
            const int chunk_start = chunk * warp_size;

            // --- Per-lane global source address ---
            // Flat index of the element this lane handles.
            const int flat_idx = chunk_start + lane;
            const int lr       = flat_idx / num_cols; // row in LDS tile
            const int lc       = flat_idx % num_cols; // col in LDS tile

            // Map to global coordinates (swap row/col for transposed tiles).
            const int gr = global_row_start + (transpose ? lc : lr);
            const int gc = global_col_start + (transpose ? lr : lc);

            // Each lane computes its own global address (VADDR, per-lane).
            const float* src_lane = &global_src[gr * global_stride + gc];

            // --- Wave-uniform LDS destination ---
            // M0 points to the start of this 64-element chunk.
            // The hardware adds lane * sizeof(float) implicitly, so
            // lane k's value lands at lds_dest[chunk_start + k].
            float* dst_chunk = &lds_dest[chunk_start];

            __builtin_amdgcn_global_load_lds(
                const_cast<float*>(src_lane), // per-lane global address
                dst_chunk,                    // wave-uniform LDS base (→ M0)
                4,                            // size: 4 bytes per lane
                0,                            // offset: immediate, 0
                0);                           // aux: wave scope, temporal
        }
    }

    __device__ static void prologue(SharedStorage& /*smem*/,
                                    const float* /*A*/, const float* /*B*/,
                                    int /*block_row_start*/, int /*block_col_start*/,
                                    int /*m*/, int /*n*/, int /*k*/,
                                    int /*tid*/, int /*num_threads*/)
    {
    }

    __device__ static void prefetch(SharedStorage& smem,
                                    const float*   A,
                                    const float*   B,
                                    int            t_next,
                                    int block_row_start, int block_col_start,
                                    int m, int n, int k,
                                    int tid, int num_threads)
    {
        const int t = t_next; // single-buffer: t_next == t

        // tile_a: row-major gather [block_tile_m × k_tile_size]
        // LDS flat layout matches global row-major order, so transpose=false.
        direct_load_flat(&smem.tile_a[0][0], A,
                         block_tile_m, k_tile_size,
                         block_row_start, t * k_tile_size,
                         k,
                         /*transpose=*/false,
                         tid, num_threads);

        // tile_b_T: transposed gather [block_tile_n × k_tile_size]
        //   LDS layout: tile_b_T[col_within_tile][k_idx]
        //   Global read: B[k_idx_global][col_global]
        //   transpose=true swaps lr↔lc in the global coordinate mapping.
        direct_load_flat(&smem.tile_b_T[0][0], B,
                         block_tile_n, k_tile_size,
                         t * k_tile_size, block_col_start,
                         n,
                         /*transpose=*/true,
                         tid, num_threads);
    }

    __device__ static void acquire(int /*t*/) { __syncthreads(); }
    __device__ static void release(int /*t*/) { __syncthreads(); }
    __device__ static int  buf_idx(int /*t*/) { return 0; }

    __device__ static float* get_tile_a(SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_a[0][0];
    }
    __device__ static float* get_tile_b_T(SharedStorage& smem, int /*buf*/)
    {
        return &smem.tile_b_T[0][0];
    }
};

Note

Because the arithmetic is unchanged (full FP32 scalar outer-product), the DirectLoadPolicy variant produces results identical to the other two. The only difference is the data path during the tile load phase.

For the full builtin reference - signatures, parameter tables, address calculation formulas, and cache policy encoding - see Global-to-LDS builtins for AMD GPUs.

What to observe:

Counter	What to look for
Kernel duration (kernel-trace CSV)	Should be similar to or better than SingleBufPolicy (same tile parameters mean same DRAM traffic).
VGPR_Count (kernel-trace CSV)	From --kernel-trace CSV: should be lower for DirectLoadPolicy than for SingleBufPolicy (tile data bypasses VGPRs during load).
LDS_Block_Size (kernel-trace CSV)	From --kernel-trace CSV: same as SingleBufPolicy (both use one buffer pair).

Further reading

The following resources provide deeper coverage of the tools and hardware referenced in this tutorial.

• rocprofv3 documentation - detailed guide to timeline and counter profiling.

• AMD GPU architecture guides (ISA references) - VGPR budgets, LDS bank geometry, and wavefront scheduling details for each architecture family.