Add memory-bound GEMM kernel and update bandwidth model

docs/memory_bandwidth_model.md

@@ -77,7 +77,8 @@ keep simultaneously in flight depends on how quickly it issues new requests rela
 how long each request takes to complete:
 
 ```
-num_req_per_wave = min(hbm_latency / compute_latency, num_stages - 1)
+num_req_per_wave = min(hbm_latency / compute_latency,
+                       num_stages * hbm_latency / (hbm_latency + compute_latency))
 ```
 
 Note that `num_req_per_wave` is the *average* number of in-flight memory requests over
@@ -86,6 +87,27 @@ a period of `hbm_latency`. During that window, the wave issues
 remains outstanding for `hbm_latency` cycles. The average occupancy of the memory
 pipeline is therefore this ratio, capped by the available buffers.
 
+The first term (`hbm_latency / compute_latency`) is the issue-rate limit: the wave
+can issue one request every `compute_latency` cycles, and each request stays in flight
+for `hbm_latency` cycles.
+
+The second term is the buffer limit. Each of the `num_stages` buffers cycles through
+two phases: being filled by an in-flight memory request (`hbm_latency` cycles) and
+being consumed by compute (`compute_latency` cycles). The total cycle time per buffer
+is `hbm_latency + compute_latency`. The fraction of time each buffer has an in-flight
+request is `hbm_latency / (hbm_latency + compute_latency)`. With `num_stages` buffers,
+the average number that can hold in-flight requests at any moment is:
+
+```
+buffer_limit = num_stages * hbm_latency / (hbm_latency + compute_latency)
+```
+
+Note that when `compute_latency` is small relative to `hbm_latency`, the buffer limit
+approaches `num_stages` — almost all buffers are being filled at any given time because
+consumption is nearly instantaneous. Conversely, when `compute_latency` is comparable to
+`hbm_latency`, the buffer limit approaches `num_stages / 2` — each buffer spends roughly
+half its time being consumed and half being filled.
+
 The three parameters here are:
 
 - **`hbm_latency`** — how long a memory request remains outstanding (hardware property).
@@ -94,14 +116,11 @@ The three parameters here are:
   as the total loop iteration time, which may be longer if the loop is HBM-latency bound.
   `compute_latency` determines the fastest rate at which a wave *could* issue new memory
   requests.
-- **`num_stages`** — number of software pipeline buffers (e.g., LDS buffers). Of these,
-  one buffer is always occupied by the consumer (compute reads from it), so at most
-  `num_stages - 1` buffers are available for in-flight memory requests.
+- **`num_stages`** — number of software pipeline buffers (e.g., LDS buffers).
 
 When `compute_latency` is large relative to `hbm_latency`, the wave issues requests too
 infrequently to fill the memory pipeline. When `compute_latency` is small, the wave can
-issue at a high rate but becomes limited by the number of available buffers
-(`num_stages - 1`).
+issue at a high rate but becomes limited by the number of available buffers.
 
 ### 2.2 Request Size
 
@@ -134,7 +153,8 @@ effective_compute_latency = waves_per_simd × compute_latency
 This feeds back into the in-flight count:
 
 ```
-num_req_per_wave = min(hbm_latency / effective_compute_latency, num_stages - 1)
+num_req_per_wave = min(hbm_latency / effective_compute_latency,
+                       num_stages * hbm_latency / (hbm_latency + effective_compute_latency))
 ```
 
 This creates a tension: increasing occupancy adds more waves but slows each one down.
@@ -224,12 +244,13 @@ iter_latency = max(hbm_latency / effective_pipeline_depth, effective_compute_lat
 ```
 
 where `effective_pipeline_depth` accounts for the TCP size limit (Section 2.4).
-The kernel has `num_stages - 1` buffers available for in-flight requests per wave,
+The kernel has `num_stages * hbm_latency / (hbm_latency + effective_compute_latency)`
+buffers available on average for in-flight requests per wave (Section 2.1),
 but the total in-flight bytes across all waves on a CU cannot exceed the TCP capacity:
 
 ```
 effective_pipeline_depth = min(
-    num_stages - 1,
+    num_stages * hbm_latency / (hbm_latency + effective_compute_latency),
     floor(CU_inflight_limit / (num_active_waves_per_CU × data_per_request_per_wave))
 )
 ```
@@ -293,3 +314,140 @@ to the same steady-state formula with `num_req_per_wave = effective_pipeline_dep
 
 For smaller problems or kernels with deep pipelines, the fixed overhead of prologue
 and epilogue becomes significant and the E2E view reveals the actual performance gap.
+
+---
+
+## 4. TCP Utilization Efficiency and Redundant Data
+
+The analysis in Sections 2 and 3 treats all in-flight bytes equally, but not all
+bytes contribute equally to useful work. In many kernels — particularly GEMMs with
+skinny M — some of the data loaded per iteration is **redundant** across workgroups,
+while other data is **unique**. The TCP capacity is shared between both, so redundant
+data reduces the fraction of TCP available for useful work.
+
+### 4.1 Redundant vs. Unique Data in GEMM
+
+In a GEMM C = A × B with shape [M, N] = [M, K] × [K, N], each workgroup computes
+a tile of C of size [BLOCK_M, BLOCK_N]. Per K-loop iteration, it loads:
+
+- **A tile** [BLOCK_M, BLOCK_K]: shared across all workgroups along the N dimension.
+  Every workgroup loads the same A data, making it redundant.
+- **B tile** [BLOCK_K, BLOCK_N]: unique per workgroup. Each workgroup loads a
+  distinct slice of B.
+
+For memory-bound GEMMs with small M (skinny matrices), we typically set BLOCK_M = M
+so that a single tile covers the entire M dimension. This means all workgroups
+partition work only along the N dimension. Every workgroup loads the same A rows
+from HBM but loads a distinct slice of B. As a result, A loads are fully redundant
+across all workgroups while B loads are fully unique.
+
+The total HBM traffic is:
+
+```
+actual_hbm_bytes = num_workgroups × (A_tile_bytes + B_tile_bytes) × num_iters
+                 + num_workgroups × C_tile_bytes
+```
+
+But the unique problem bytes are:
+
+```
+unique_bytes = M × K × elem_bytes + K × N × elem_bytes + M × N × elem_bytes
+```
+
+The ratio between these defines the **data reuse efficiency**:
+
+```
+reuse_efficiency = unique_bytes / actual_hbm_bytes
+```
+
+When reporting achieved bandwidth as `unique_bytes / time`, the kernel can never
+exceed `peak_hbm_bandwidth × reuse_efficiency`.
+
+### 4.2 TCP Utilization Efficiency
+
+Since the TCP (32 KB on GFX9) limits the total in-flight bytes per CU, and both
+redundant and unique data compete for this capacity, only a fraction of TCP
+is doing useful work:
+
+```
+tcp_efficiency = B_tile_bytes / (A_tile_bytes + B_tile_bytes)
+```
+
+This creates a direct relationship between tile shape and achievable bandwidth.
+Since A tile and B tile share the same BLOCK_K, the TCP efficiency depends on
+the ratio of BLOCK_M (= M for skinny GEMMs) to BLOCK_N:
+
+```
+tcp_efficiency = BLOCK_N / (M + BLOCK_N)
+```
+
+Increasing BLOCK_N improves TCP efficiency:
+
+| M   | BLOCK_N | BLOCK_K | A tile | B tile | TCP efficiency |
+|-----|---------|---------|--------|--------|----------------|
+| 32  | 32      | 128     | 8 KB   | 8 KB   | 50%            |
+| 32  | 64      | 128     | 8 KB   | 16 KB  | 67%            |
+| 32  | 128     | 64      | 4 KB   | 16 KB  | 80%            |
+| 32  | 256     | 64      | 4 KB   | 32 KB  | 89%            |
+
+However, increasing BLOCK_N has a cost: it reduces the number of workgroups
+for a given problem size (`num_workgroups = N / BLOCK_N`), which may leave CUs
+idle. This creates a **three-way tradeoff**:
+
+1. **TCP efficiency** — wants large BLOCK_N (relative to M)
+2. **CU utilization** — wants many workgroups (small BLOCK_N or large N)
+3. **TCP capacity** — the total tile data in the pipeline, roughly
+   `(A_tile + B_tile) × num_stages`, must fit within the TCP size (32 KB on
+   GFX9). After M and BLOCK_N are determined, BLOCK_K is the remaining knob
+   to control the per-iteration data size within this limit
+
+### 4.3 Strategies for Improving TCP Efficiency Without Reducing CU Utilization
+
+The tradeoff in Section 4.2 arises because each workgroup processes exactly one
+output tile. Several strategies can break this coupling:
+
+**A-reuse via inner N-tile loop.** Within each K iteration, load the A tile once
+into registers, then process multiple B tiles sequentially:
+
+```
+for k in range(K / BLOCK_K):
+    a = load_A(k)                        # load once
+    for t in range(N_TILES_PER_WG):
+        b = load_B(k, base_n + t)        # load each B tile
+        acc[t] = mfma(a, b, acc[t])      # reuse a registers
+```
+
+This amortizes the A tile load over `N_TILES_PER_WG` B tile loads, giving an
+effective TCP efficiency of:
+
+```
+tcp_efficiency = (N_TILES × B_tile) / (A_tile + N_TILES × B_tile)
+```
+
+which approaches 100% as `N_TILES` grows. The cost is additional accumulator
+registers — each N-tile requires `BLOCK_M × BLOCK_N × 4` bytes of fp32 VGPRs.
+The number of workgroups also decreases by `N_TILES_PER_WG`, requiring
+proportionally larger N to fill all CUs.
+
+**Split-K.** Partition the K dimension across multiple workgroups. Each workgroup
+computes a partial sum over a subset of K, and a final reduction combines them.
+This increases the workgroup count without changing BLOCK_N, allowing both high
+TCP efficiency (large BLOCK_N) and high CU utilization. The cost is the reduction
+overhead and additional global memory traffic for partial sums.
+
+### 4.4 Integrating TCP Efficiency into the Bandwidth Model
+
+The steady-state bandwidth from Section 2.4 gives the raw HBM bandwidth the kernel
+can sustain. To predict the *useful* (unique-bytes) bandwidth, multiply by the
+data reuse efficiency:
+
+```
+BW_useful = BW_raw × reuse_efficiency
+```
+
+where `reuse_efficiency` accounts for redundant A loads across workgroups.
+Equivalently, in the E2E model, `total_bytes` should be the unique problem bytes
+(A + B + C counted once), while `total_cycles` reflects the actual execution time
+which includes cycles spent loading redundant A data. The gap between raw HBM
+bandwidth and useful bandwidth is entirely determined by the A/B tile ratio and
+the number of workgroups sharing the same A data.

kernels/memory_bound/gemm/README.md

@@ -0,0 +1,149 @@
+# Memory-Bound GEMM Kernel for AMD GFX950 (MI350/CDNA4)
+
+This directory contains a Triton Gluon kernel optimized for **memory-bound GEMM**
+on AMD GFX950 (MI350). The target workload is skinny matrix multiplication with
+small M (e.g., M=32), where the kernel is bottlenecked by HBM bandwidth rather
+than compute throughput.
+
+## Problem Setup
+
+```
+C[M, N] = A[M, K] x B[K, N]    (fp16, fp32 accumulation)
+```
+
+With M=32, the arithmetic intensity is low and the kernel becomes memory-bound.
+The goal is to maximize HBM bandwidth utilization on a GPU with 8.1 TB/s peak
+HBM bandwidth and 256 CUs.
+
+## Kernel Configuration
+
+| Parameter | Value |
+|-----------|-------|
+| BLOCK_M   | 32    |
+| BLOCK_N   | 256   |
+| BLOCK_K   | 64    |
+| Warps     | 4 (2x2) |
+| Pipeline stages | 2 |
+| A tile size | 32 x 64 = 4 KB |
+| B tile size | 64 x 256 = 32 KB |
+| Total per iteration | 36 KB |
+
+## Key Optimizations
+
+### 1. Async Copy Pipeline
+
+The kernel uses a 2-stage software pipeline with CDNA4 async copy instructions:
+
+```
+buffer_load_dwordx4 → LDS → ds_read_b128 → MFMA
+```
+
+The prologue fills both pipeline buffers before entering the main loop. Each
+iteration overlaps the current MFMA compute with the next iteration's HBM loads.
+
+Ablation testing (removing all `ds_read_b128` and `s_waitcnt lgkmcnt` from the
+assembly) confirmed that the kernel is **purely HBM-bandwidth bound** — LDS
+reads are completely hidden behind HBM load latency.
+
+### 2. `.cg` Cache Modifier on B Loads
+
+The B tile is unique per workgroup (each workgroup loads a different slice of B
+along the N dimension). The `.cg` (cache global) modifier bypasses L2 caching for
+B loads, avoiding cache pollution since B data is not reused across workgroups.
+
+The A tile is shared across all workgroups along the N dimension (redundant loads),
+so it benefits from caching and uses the default cache policy.
+
+### 3. Large BLOCK_N for TCP Utilization Efficiency
+
+On GFX9, the TCP (Texture Cache per Pipe, 32 KB per CU) limits the total in-flight
+bytes per CU. Both A and B tile data compete for TCP capacity, but only B tile
+data contributes unique (non-redundant) work. The TCP utilization efficiency is:
+
+```
+tcp_efficiency = B_tile_bytes / (A_tile_bytes + B_tile_bytes)
+```
+
+With BLOCK_N=256, BLOCK_K=64:
+- A tile = 4 KB, B tile = 32 KB
+- TCP efficiency = 32 / (4 + 32) = **89%**
+
+This is a significant improvement over smaller BLOCK_N values:
+
+| BLOCK_K | BLOCK_N | A tile | B tile | TCP efficiency |
+|---------|---------|--------|--------|----------------|
+| 128     | 32      | 8 KB   | 8 KB   | 50%            |
+| 128     | 64      | 8 KB   | 16 KB  | 67%            |
+| 64      | 128     | 4 KB   | 16 KB  | 80%            |
+| 64      | 256     | 4 KB   | 32 KB  | 89%            |
+
+### 4. Per-Thread Load Width Constraint
+
+Each thread must issue at least one full-width `buffer_load_dwordx4` (8 fp16
+elements = 16 bytes). This constrains the minimum tile size:
+
+```
+BLOCK_M × BLOCK_K / num_threads ≥ 8 elements
+```
+
+With BLOCK_M=32, BLOCK_K=64, and 256 threads (4 warps × 64 lanes):
+- A tile: 32 × 64 / 256 = 8 elements → 1 buffer_load per thread
+- B tile: 64 × 256 / 256 = 64 elements → 4 buffer_loads per thread
+- Total: 5 buffer_loads per iteration per thread
+
+## Performance Results
+
+Benchmark configuration: M=32, K=8192, fp16, 256 CUs on GFX950.
+
+### Bandwidth vs. Workgroup Count (BLOCK_N=256)
+
+N is chosen so that `N / BLOCK_N` gives the target number of workgroups:
+
+| Workgroups | N      | Bandwidth (GB/s) |
+|------------|--------|-------------------|
+| 32         | 8192   | 1529              |
+| 64         | 16384  | 2864              |
+| 96         | 24576  | 4067              |
+| 128        | 32768  | 4580              |
+| 160        | 40960  | 4828              |
+| 192        | 49152  | 5628              |
+| 224        | 57344  | 6198              |
+| 256        | 65536  | **6442**          |
+
+Bandwidth scales nearly linearly with workgroup count, reaching **6442 GB/s**
+(79% of 8.1 TB/s peak) at 256 workgroups (one per CU).
+
+### Pipeline Stages Comparison (256 WGs, K=8192)
+
+| Stages | Bandwidth (GB/s) |
+|--------|-------------------|
+| 2      | 6430              |
+| 3      | 6442              |
+
+The two configurations perform identically because the TCP (32 KB) is the
+binding constraint — not the number of software pipeline buffers.
+
+## Running the Benchmark
+
+```bash
+# Run with default sizes (M=32, N=65536, K=1024..16384)
+python bench.py
+
+# Run a specific K value
+python bench.py --K 8192
+
+# Run with rocprofv3 for external timing (uses rotating buffers to defeat caching)
+rocprofv3 --kernel-trace -- python bench.py --K 8192 --rocprof
+```
+
+## Files
+
+- `matmul_kernel.py` — Triton Gluon kernel and Python wrapper
+- `bench.py` — Benchmark harness with correctness checking and rocprof support
+
+## Further Reading
+
+See [`docs/memory_bandwidth_model.md`](../../../docs/memory_bandwidth_model.md)
+for the bandwidth model used to analyze this kernel, including the TCP utilization
+efficiency framework and strategies for improvement (A-reuse via inner N-tile
+loop, Split-K).