reasoning
Let me analyze the problem carefully.
The reference model is doing a matrix multiplication: A (M, N) @ B (N, M) where:
- M = 32768 (16384 * 2)
- N = 32 (16 * 2)
Wait, looking at get_inputs():
```python
A = torch.rand(M, N) # 32768 x 32
B = torch.rand(N, M) # 32 x 32768
```
So the matmul is A @ B = (32768 x 32) @ (32 x 32768) = (32768 x 32768)
Let me think about this more carefully. The shapes are:
- A: (32768, 32)
- B: (32, 32768)
- Output: (32768, 32768)
This is a tall-skinny matrix multiplication. But actually the output is a square matrix of size 32768 x 32768. That's quite large!
For this kind of GEMM, cuBLAS should already be very optimized. Let me think about whether we can beat it with a custom kernel.
Looking at the compute:
- M = 32768, K = 32, N = 32768
- Flops = 2 * M * K * N = 2 * 32768 * 32 * 32768 = 2 * 2^15 * 2^5 * 2^15 = 2 * 2^35 = 2^36 = 68.7 GFLOPs
To read A: 32768 * 32 * 4B = 4 MB
To read B: 32 * 32768 * 4B = 4 MB
To write C: 32768 * 32768 * 4B = 4 GB
Wait, writing 4GB is the dominant memory transfer! The arithmetic intensity is very low.
Actually let me re-read the problem description:
"A (torch.Tensor): Input matrix of shape (M, K) or (K, M) where M >> N or N >> M."
Looking at get_inputs:
A = torch.rand(M, N)
B = torch.rand(N, M)
So A is (32768, 32) and B is (32, 32768). The matmul A @ B gives (32768, 32768).
But wait - the dtype defaults to float32, so 4GB output. This is memory bound.
Hmm, can a custom kernel beat cuBLAS for this? The output is huge and the compute is relatively small. Let me think...
Actually, I should just try the reference first. But I don't have a way to run the reference. Let me think about what kernel would be optimal.
For (32768 x 32) @ (32 x 32768), each output element is a dot product of length 32. With 4GB of output, this is definitely memory-bound on the write side.
The arithmetic intensity is:
68.7 GFLOPs / (8MB read + 4GB write) ≈ 68.7 / 4.008 ≈ 17 FLOPs/byte
Wait, H100 memory bandwidth is 3.35 TB/s = 3350 GB/s
Compute is 67 TFLOPS = 67000 GFLOPS
Roofline turning point = 67000 / 3350 ≈ 20 FLOPs/byte
So we're slightly below the roofline turning point, maybe slightly compute-bound or memory-bound depending on caching.
Actually, since the output is 4GB, and we need to write all of it, if we're memory bound:
Time = 4GB / 3.35 TB/s = 4 / 3350 s = 1.19 ms
GFLOPS = 68.7 / 0.00119 = 57.7 TFLOPS
Wait that's actually very close to peak FP32. Let me recalculate:
4GB / 3.35TB/s = 4 / 3350 seconds = 0.001194 seconds = 1.194 ms
68.7 GFLOPs / 1.194 ms = 57.5 GFLOPS... wait, I said TFLOPS but it's GFLOPs. Let me convert:
68.7 GFLOPs / 0.001194 s = 57.5 GFLOPS, not TFLOPS.
57.5 GFLOPS vs 67 TFLOPS peak. So this is definitely memory bound.
Wait, I miscalculated. 2 * 32768 * 32 * 32768:
32768 = 2^15
2^15 * 2^5 = 2^20 = 1,048,576
2^20 * 2^15 = 2^35 = 34,359,738,368
Times 2 = 68,719,476,736 operations = 68.7 GOPS (not GFLOPS, it's just operations, but each is a multiply-add so it's FLOPs).
68.7 GFLOPs / 1.194 ms = 57.5 GFLOPS.
Peak FP32 is 67 TFLOPS = 67,000 GFLOPS.
So utilization is 57.5 / 67000 = 0.08%. Very very memory bound.
For memory-bound kernels, the key is to maximize memory throughput. A custom kernel that writes coalesced and reads from shared memory might help, but cuBLAS should already be very good at this.
Wait, but let me reconsider the shapes. Actually, the code says A is (M, N) and B is (N, M):
A = torch.rand(M, N) # (32768, 32)
B = torch.rand(N, M) # (32, 32768)
torch.matmul(A, B) = (32768, 32768)
But is there any chance they meant something else? Let me re-read the model docstring:
"Input matrix of shape (M, K) or (K, M) where M >> N or N >> M"
Hmm, the inputs use M and N differently. A is (M, N) and B is (N, M) in the get_inputs. So K = N for A and K = N for B as well? That would mean:
A : (M, K) = (32768, 32)
B : (K, N) = (32, 32768)
C : (M, N) = (32768, 32768)
Yes.
For this shape, each element of C is computed from a dot product of two vectors of length 32. With 1 billion output elements and only 32 multiplications each, the computation is tiny compared to the memory traffic.
Actually, for a kernel: each thread could compute one or more output elements. We want to:
1. Read A and B in a coalesced way
2. Compute the dot product
3. Write C coalesced
For this case, since K=32 is very small, one interesting approach is to have each thread read 32 elements from A (one row) and compute a stride of outputs. But B needs to be read multiple times.
Actually, since it's memory bound, the most important thing is writing C efficiently. Each element of C requires reading one row of A and one column of B. Since columns of B are very long (32768), reading them is strided and inefficient.
Wait, let me reconsider. B is (32, 32768). If we read a column of B, that's reading with a stride of 32768 elements, which is very bad.
But cuBLAS will handle this by transposing B first, or by doing some blocked algorithm.
Actually, in the get_inputs function, B = torch.rand(N, M). But the matmul docstring says B is of shape (K, N) or (N, K). So B could be (32, 32768), meaning K=32, N=32768.
For C = A @ B with A:(32768,32) and B:(32,32768):
- A is row-major, each row is contiguous
- B is row-major, each column has stride 32768
A column-major read from B is extremely inefficient. One obvious optimization is to transpose B first so that each "column" in the transposed B is contiguous. But that requires an extra memory operation.
Actually, B is just (32, 32768), so it's only 4MB. We could load all of B into shared memory or registers!
Let me think about a custom kernel:
- Load all of B (32 x 32768 = 1M floats = 4MB) into shared memory or cache? That's too big for shared memory (227KB max).
- But we can tile B.
Actually, 4MB could fit in L2 cache easily (50MB on H100). But for a cooperative kernel, we need to think about how to use it efficiently.
Let me think about a different interpretation. What if A were transposed or B? The get_inputs creates:
A = torch.rand(M, N) # (32768, 32)
B = torch.rand(N, M) # (32, 32768)
And C = A @ B = (32768, 32768)
For a thread computing an output tile:
- Read a row of A: 32 floats, contiguous
- For each column j in the output tile, read B[k, j] for k in 0..31
- B[k, j] with k varying is contiguous (since B is row-major, consecutive k means... wait, B is (32, 32768). B[k, j] for fixed j and varying k means we read B[0,j], B[1,j], B[2,j], etc. The stride between B[k,j] and B[k+1,j] is 32768. So columns are strided by 32768.
This is terrible for memory access! Reading a column of B means 32 memory accesses with 32768-element strides.
So actually a huge opportunity is to transpose B first! B is just 4MB. If we transpose B to get B.T of shape (32768, 32), then:
- We need to compute C[i,j] = sum_k A[i,k] * B.T[j,k]
- Now for a fixed j, the values B.T[j, :] are contiguous (32 floats)
- For reading A[i, :], they are also contiguous
So both A[i,:] and B.T[j,:] are contiguous 32-float vectors. We can compute each output element with two contiguous reads!
But the naive approach would transpose B first. However, transposing B is another 4MB memory operation. For a memory-bound overall operation, maybe we can fuse the transpose with the computation?
Or we can work with tiles:
- A thread block computes a tile of C, say 128x128
- It loads 128 rows of A (128 x 32 floats = 16KB)
- It loads 128 columns of B... but these are strided!
- Instead, load rows of B? Wait, B is (32, 32768). Rows of B are contiguous! So for computing C[i,j] = sum_k A[i,k] * B[k,j].
For a tile of C of size TileM x TileN:
- For each k in 0..31:
- Load A[i_offset : i_offset+TileM, k] - TileM floats, strided by 32... actually for A (32768, 32), A[:, k] has stride 32. Not great.
- Load B[k, j_offset : j_offset+TileN] - TileN floats, contiguous!
Hmm wait, A[i,k] for fixed k and varying i: A is (32768, 32), so A[:,k] has stride 32. That's actually decent for coalescing since threads in a warp will read consecutive i's. If each thread handles one i, then consecutive threads read consecutive memory locations (offset by 32 * 4 = 128 bytes).
For B[k, j_offset:j_offset+TileN], consecutive j's are contiguous.
So for each k, we can:
- Load a slice of A (for all i in the tile)
- Load a slice of B (for all j in the tile)
- And perform outer product to accumulate into C tile!
This is actually a very efficient approach! For each k from 0 to 31:
- Thread block loads A_tile[i] = A[i_tile_start + i, k] for i = 0..TileM-1
- Thread block loads B_tile[j] = B[k, j_tile_start + j] for j = 0..TileN-1
- Each thread accumulates A_tile[local_i] * B_tile[local_j] into its register
This is an outer-product formulation of GEMM!
Since K=32 is so small, we only need 32 iterations. Each iteration loads O(TileM + TileN) elements and does O(TileM * TileN) work.
For memory bound with large output, this approach is very good because:
1. We read each element of A exactly once (32 * 32768 = 1M floats = 4MB)
2. We read each element of B exactly once (32 * 32768 = 1M floats = 4MB)
3. We write each element of C exactly once (1B floats = 4GB)
The reads are coalesced! For A[:,k], consecutive threads read consecutive elements. For B[k,:], we load contiguous segments.
Wait, is this better than what cuBLAS does? cuBLAS likely also does something similar. But the outer-product formulation with K=32 being very small means we can:
- Keep the entire accumulator tile in registers
- Only need to load one "panel" from A and one from B per k iteration
Since K is only 32 and we need to write 4GB, the bottleneck is truly the write bandwidth. But we want to maximize throughput.
Let me think about optimal tile sizes. If TileM=128 and TileN=128, then each thread block computes 128*128 = 16384 outputs. The accumulator is 16384 floats = 64KB. Each thread would hold 16384/256 = 64 floats (if 256 threads). That's fine.
Actually, we can use much larger tiles if we want. But we also need to load A and B into shared memory.
For each k, we load:
- A panel: 128 floats = 512 bytes (for TileM=128)
- B panel: 128 floats = 512 bytes (for TileN=128)
Total shared memory per k: 1KB. We need double buffering? Not really necessary since K=32 is small and we can just load both into shared memory each iteration.
Wait, we can also load A into shared memory once and reuse it across all j tiles! That's a standard optimization in GEMM. If we have a grid where each block handles a TileM x TileN tile of C, blocks in the same row of the grid share the same A rows. On H100 with L2 cache, this might happen automatically.
But for K=32, A is only 4MB total. All of A easily fits in L2. B is also 4MB. So both inputs fit in cache!
Actually, A is (32768, 32) = 1,048,576 floats = 4MB
B is (32, 32768) = 1,048,576 floats = 4MB
Total input = 8MB
L2 cache on H100 is 50MB. The entire working set of inputs fits in L2!
So the key is to write C efficiently and maximize memory parallelism.
Let me try writing a kernel that computes an outer product for each k. Let's use:
- Each block computes a TileM x TileN tile
- Thread block size = TileM (one thread per row in the tile) or some 2D configuration
- For each k in 0..K-1 (K=32):
- Load A[i_block + tid, k] into shared memory (or registers)
- Load B[k, j_block + ...] into shared memory
- Compute outer product and accumulate
Actually, with a 2D thread block, say 16x16 = 256 threads:
- For each k:
- Thread (tx, ty) loads A[i0 + tx, k] and B[k, j0 + ty]... but we need many more than 16 elements.
Let me think more carefully. If TileM=128, TileN=128:
- We have 128 rows and 128 columns to compute.
- We can use 256 threads organized as 16x16.
- Each thread computes a sub-tile of 8x8 = 64 outputs.
- For each k in 0..31:
- The block collaboratively loads 128 elements of A column (A[i0:i0+128, k]) into shared memory
- The block collaboratively loads 128 elements of B row (B[k, j0:j0+128]) into shared memory
- Each thread reads one A value and one B value from shared memory, multiplies, and accumulates into one output. But each thread needs 64 outputs...
Actually, if each thread computes 64 outputs (8x8 sub-tile), then for each k:
- Thread reads A_val for its 8 rows? No, for each output C[local_i, local_j], it needs A[local_i, k] * B[k, local_j].
- For an 8x8 sub-tile with fixed local_i and varying local_j, all outputs in the same row share the same A value! And all outputs in the same column share the same B value.
- So for each k, thread needs 8 different A values (one per row of its sub-tile) and 8 different B values (one per column of its sub-tile).
- Total: 8A + 8B reads from shared memory per k, and 64 multiply-adds.
- For K=32: 256*(8+8) shared memory reads and 256*64*32 = 524,288 FMAs per block.
Actually the FMA count is TileM * TileN * K = 128 * 128 * 32 = 524,288 per block.
Shared memory reads: 128 * 32 A-reads + 128 * 32 B-reads = 4096 + 4096 = 8192 reads of float. With 256 threads, that's 32 reads per thread on average. Actually if each thread reads 8 A's and 8 B's per k, that's 16 * 32 = 512 reads per thread over all k.
Wait, let me think about the parallel loading:
- 256 threads, need to load 128 A elements and 128 B elements.
- 128 threads load A elements (each loads 1 float)
- 128 threads load B elements (each loads 1 float)
- All threads have their data in shared memory.
Then each thread computes 8x8 outputs. For each k, it fetches 8 A-values and 8 B-values from shared memory.
But we could also let each thread hold its A and B values in registers. Thread with local tile starting at (bi, bj):
- For each k: A_vals[0..7] = A[i0+bi+0..7, k], B_vals[0..7] = B[k, j0+bj+0..7]
- Then do nested loops: for a in 0..7, for b in 0..7: accum[a][b] += A_vals[a] * B_vals[b]
This is good because it amortizes the shared memory reads (or we can load directly from global memory if cached).
Wait, actually we can load A and B directly into registers without shared memory for the tiles!
- A[:, k] has good locality for threads in a warp. If threads in a warp handle different i's but same bi sub-tile offset, they read strided by 32 floats = 128 bytes. For a warp of 32 threads, consecutive threads read 128 bytes apart. Is that coalesced?
Actually, for A of shape (32768, 32), if thread i reads A[i0 + i, k], the addresses are:
&A[0, k] + (i0 + i) * 32 * 4 bytes
=&A[0, k] + i0 * 128 + i * 128 bytes
Consecutive threads in a warp access addresses 128 bytes apart! For H100, coalescing requires accesses within a 128-byte cache line. So 32 threads * 128 bytes = 4KB span per warp. That's way larger than a cache line. This is not fully coalesced!
Actually wait, for a warp reading A[i0 + warp_tid, k] where warp_tid = 0..31:
Address = base + warp_tid * 128 bytes
H100 has 128-byte L2 cache lines. So each thread is reading a different cache line! That's terrible. We'd get 32 cache line reads for 32 floats = 128 bytes of useful data, but 32 * 128 = 4096 bytes transferred! 32x waste.
So we should transpose how threads load A. We should have threads in a warp read contiguous elements.
Actually, for A with shape (32768, 32), a column A[:, k] consists of elements at indices [0*32+k, 1*32+k, 2*32+k, ...]. Consecutive elements are 32 floats = 128 bytes apart.
To read a column coalesced, we'd want each warp to read a segment where the 32 consecutive elements are read by 32 threads. But they're 128 bytes apart...
Wait, on H100, global memory loads are done at L2 cache line granularity (128 bytes). When threads in a warp access memory, the hardware groups accesses to the same 128-byte segment. If each thread accesses a different 128-byte segment, we get 32 transactions.
So a column of A cannot be read fully coalesced with one float per thread. We need vectorized loads!
If each thread loads 32 consecutive floats? No, consecutive floats in the column are 128 bytes apart.
Hmm, wait. Let me reconsider A's layout. A = torch.rand(M, N) = torch.rand(32768, 32). In PyTorch, this is row-major.
A[i, k] is at address A + (i * N + k) * 4 = A + i * 128 + k * 4.
For a fixed k and varying i:
Address = base + i * 128
So consecutive rows are 128 bytes apart. To read a column segment of 32 rows, we need 32 cache lines.
However, we can read 32 rows with vectorized loads differently. What if we read 4 rows at a time per thread with float4? The addresses would still be far apart.
Actually, for this specific case, we might want to transpose A first. A is 4MB. Transposing A makes A.T of shape (32, 32768), so rows of A.T are columns of A and are contiguous in memory.
But transposing is an extra kernel launch and memory operation. For 4MB it's cheap, but still adds overhead.
Alternatively, we can have the thread block load a chunk of A column using vectorized loads in a smarter way. But I don't see how to fully coalesce it since the stride is exactly 128 bytes = cache line size.
Wait! If the stride is exactly 128 bytes, then loading 32 consecutive elements hits exactly 32 consecutive cache lines. If adjacent warps load adjacent 32-element chunks, those are also consecutive cache lines. The memory controller can handle this efficiently because the cache lines are contiguous in address space!
Actually, reading 32 consecutive cache lines is as efficient as reading a single 4KB block. The memory system on H100 handles bursts of sequential cache line requests very well. So even though each thread reads a different cache line, the overall pattern is sequential at the warp/block level. This should actually achieve close to peak bandwidth!
Let me verify: A[:, k] for k=0 is elements at addresses: 0, 128, 256, ..., 4096. These are cache lines 0, 1, 2, ..., 32. A warp reading 32 elements triggers 32 cache line loads. But these are to cache lines 0-31, which are perfectly sequential. The L2 and memory controllers will merge these into efficient burst reads.
In fact, a memory request for 32 sequential cache lines is probably the best pattern possible! The memory controller can issue one big read request.
Wait, but each thread in a warp is requesting a different cache line. Is the memory controller smart enough to see this? Modern GPUs do have coalescing logic that groups accesses. If accesses fall in sequential cache lines, they should be merged.
Actually, H100 has improved memory subsystems. I believe sequential cache line patterns are handled well. But let me think about whether we can do even better.
What about reading A in a transposed pattern? If we read A as float4 from A[i, k:k+4], the addresses are:
A + i*128 + k*4, A + i*128 + (k+1)*4, A + i*128 + (k+2)*4, A + i*128 + (k+3)*4.
For a warp where each thread has a different i but same k:
Addresses: base_i + i*128 for i = 0..31. Again 32 different cache lines.
Hmm. What if we have a warp read a 32x32 tile of A? A[i0..i0+31, k0..k0+31]. That would require reading 32x32 = 1024 floats. Each thread could read 32 floats along the k dimension for a fixed i. Those 32 floats would be contiguous in memory (128 bytes, 2 cache lines). So each thread reads 2 cache lines, and 32 threads read 64 cache lines. But the cache lines for consecutive i are spaced...
Actually, A[i, k0..k0+31] for thread i: addresses from i*128 + k0*4 to i*128 + k0*4 + 124. That's 128 bytes, exactly one cache line if aligned, or 2 if straddling. Wait, 32 floats = 128 bytes. Yes! Exactly one cache line.
So if thread i loads A[i, k0:k0+32] as a float4×8 or something, it reads 128 bytes from address i*128 + k0*4. For consecutive i, these are consecutive cache lines. A warp of 32 threads reads 32 consecutive cache lines. Perfect!
But in our outer product algorithm, for each k we only need A[:, k], not A[:, k:k+32].
Hmm, for K=32, we could process 32 k's at once! Instead of iterating k=0..31 one at a time, iterate over the entire K dimension in one go:
- Load A[i0..i0+TileM-1, :] (32 floats per row) into shared memory or registers
- Load B[:, j0..j0+TileN-1] (32 floats per column, but wait, B is (32, 32768), so B[k, j0..j0+TileN-1] is contiguous for each k)
Actually yes! For a tile:
- Load A_panel[TileM][32] - each row is 32 contiguous floats = 128 bytes.
- Load B_panel[32][TileN] - each row is TileN contiguous floats.
If TileM=128, A_panel is 128 rows × 32 floats = 4096 floats = 16KB.
If TileN=1024, B_panel per k is 1024 floats = 4KB. For all 32 k's, that's 32 × 4KB = 128KB. Total = 144KB. That's a lot for shared memory (max ~227KB per block, and we need multiple blocks per SM).
But we can process K in chunks. Or better yet, we can load A rows into registers since each thread only needs its own row's 32 floats.
Let me think about a thread-centric design. If each thread computes a TileM_sub × TileN_sub chunk of outputs:
- Say thread computes 1 × TileN_sub outputs (one row, many columns)
- Thread loads A[i, 0..31] into registers (32 floats)
- For each j in its sub-tile, thread loads B[0..31, j] - but B[:, j] is strided by 32768!
- That's terrible.
Alternative: each thread computes TileM_sub × 1 outputs (many rows, one column):
- Thread loads B[0..31, j] into registers (32 floats, strided terribly)
- For each i in its sub-tile, loads A[i, 0..31] - contiguous 32 floats per row.
- Multiply-add.
- But B[:, j] load is terrible.
Alternative: each thread block collaborates to load B.T tiles into shared memory.
- B is (32, 32768). B.T is (32768, 32).
- For a tile of B.T of shape TileN × 32, each row is 32 contiguous floats = 128 bytes.
- This is very efficient to load!
So if we pre-transpose B in a separate kernel, the main GEMM becomes:
C = A @ B = A @ (B.T).T
But actually if B.T is (32768, 32), then matmul with A (32768, 32):
C[i,j] = sum_k A[i,k] * B[k,j] = sum_k A[i,k] * B.T[j,k]
This is a batched dot product between rows of A and rows of B.T!
- A row i: 32 contiguous floats
- B.T row j: 32 contiguous floats
- Dot product gives C[i,j]
So if we have pre-transposed B.T, then each output element is just a dot product of two contiguous 32-float vectors. This is very efficient!
But pre-transposing costs 4MB of memory traffic. For a 4GB output, this is 0.1% overhead. That's totally worth it!
Wait, but the challenge is to write a single kernel that replaces the PyTorch forward. We can launch a custom CUDA kernel that does the transpose internally and then the matmul. Or do a two-kernel approach if the overhead is acceptable.
Actually, we might be able to fuse it. Or we could use a single kernel that computes C using the outer product method directly, without explicit transpose.
Let me think about what cuBLAS actually does. For A (32768, 32) × B (32, 32768):
- If B is marked as needing transpose internally, cuBLAS might do:
1. Transpose B to B.T (32768, 32)
2. Then compute C[i,j] = dot(A[i,:], B.T[j,:])
cuBLAS has heuristics for when to transpose. It might automatically transpose if beneficial.
For our custom kernel, if we do the outer product approach:
- Don't transpose B explicitly
- For each k in 0..31:
- Load A[:, k] (this is strided by 128 bytes, but sequentially cache-line-friendly)
- Load B[k, :] (this is contiguous! 32768 floats)
- Do outer product accumulation
Wait, B[k, :] has length 32768. If a thread block computes a 128×128 tile of C:
- For each k, B[k, :] needs to be loaded.
- Each block only needs B[k, j0:j0+128], which is just 128 contiguous floats = 512 bytes.
- A[:, k] needs A[i0:i0+128, k], which is 128 floats strided by 128 bytes.
So for each k, we load ~512 bytes for the tile. Over K=32, that's ~16KB per block. Very manageable.
And each block writes 128×128 floats = 64KB.
Total memory per block = 16KB read + 64KB write = 80KB.
Total blocks = (32768/128) × (32768/128) = 256 × 256 = 65536 blocks.
Total memory = 65536 × 80KB = 5GB.
But inputs are read multiple times. Each A[:, k] is read by 256 blocks (one for each column tile). So A is read 256 times: 4MB × 256 = 1GB. Each B[k, :] is read by 256 blocks (one for each row tile). B is read 256 times: 4MB × 256 = 1GB. Output is 4GB. Total = 6GB.
Hmm, that's actually quite a bit of extra reads. But with L2 cache, subsequent reads of A[:, k] should hit cache since 256 blocks reading the same A[:, k] will be scheduled at similar times. And A[:, k] is only 128KB (32768 floats = 128KB). So it easily stays in L2!
B[k, j0:j0+128] segments are 512 bytes. There are 256 such segments per k, total 128KB per k. All segments for a given k also fit in L2. And consecutive blocks will likely have these in cache.
So the outer product approach with proper tiling could be very efficient, and the redundant reads might actually be well-cached.
But I wonder if we can do better with the dot product approach after explicit transpose.
With transposed B.T (32768, 32), each output C[i,j] is a dot product of A[i,:] (32 floats) and B.T[j,:] (32 floats).
- We can organize threads to compute C in tiles.
- A block computing 128×128 tile:
- Loads A[i0:i0+128, :] = 128 × 32 = 4KB (each row is 128 bytes, contiguous)
- Loads B.T[j0:j0+128, :] = 128 × 32 = 4KB (each row is 128 bytes, contiguous)
- Then each thread does 32 multiply-adds for each of its outputs.
Wait, if each thread computes one output, it needs 32 FMAs. A 128×128 tile = 16384 outputs. With 256 threads, each does 64 outputs = 2048 FMAs.
For the 128×128 tile:
- Load 128 rows of A, each 32 floats. All contiguous.
- Load 128 rows of B.T, each 32 floats. All contiguous.
- Keep them in shared memory.
- Each thread reads the appropriate A row and B.T row from shared memory.
- But this means shared memory needs 4KB + 4KB = 8KB.
- And we only need 32 FMAs per output...
Actually, this seems inefficient in terms of shared memory bandwidth. Each thread reads 32 floats from A and 32 from B.T just to do 32 FMAs. That's 2:1 load-to-compute ratio from shared memory.
Can we improve this? What if we vectorize the dot product? Each thread computes multiple outputs and reuses the A row data.
If thread computes 1 output: reads 32 A + 32 B.T = 64 loads for 32 FMAs.
If thread computes 8 outputs in same row: reads 32 A + 8×32 B.T = 288 loads for 256 FMAs. Worse!
If thread computes 8 outputs in same column: reads 8×32 A + 32 B.T = 288 loads for 256 FMAs.
Wait, I think I'm confusing myself. Let me think again.
For a tile approach where the whole tile fits in shared memory, we can do:
- Load A_tile[128][32] and B_tile[128][32]
- Thread computes C[i, j] by reading A_tile[i][0..31] and B_tile[j][0..31]
- But these are 32 elements each. Could we store them in registers instead?
What if K=32 is small enough that each thread just holds the A row in registers?
- Thread computing row i loads A[i, 0..31] into 32 registers.
- Then for each j in its output range, thread loads B.T[j, 0..31]... but B.T[j, :] is contiguous, and for different j, B.T[j, :] is strided by 32768 floats = 128KB.
Hmm, that's also a bad access pattern for B.T.
But in the outer product approach, B[k, :] is contiguous! We process one k at a time and each block only loads a slice of B[k, :].
I think the outer product approach is more natural for this shape. Let me think about how to implement it efficiently.
Outer Product Kernel Design:
- Grid: (NumTileM, NumTileN) where NumTileM = M / TileM, NumTileN = N / TileN
- Each block computes TileM × TileN outputs
- Thread block: e.g., 256 threads = (16, 16) or 32×8 or whatever
- For k = 0 to K-1 (K=32):
- Load A[TileM_start : TileM_start+TileM, k] into shared memory
- This is 128 floats, each 128 bytes apart.
- We can load them with vectorized or coalesced patterns.
- Actually, we can have 128 threads each load one float. The addresses are base + tid * 128. This spans 128 * 128 = 16384 bytes = 128 cache lines. Sequential cache lines! The memory system should handle this well.
- Load B[k, TileN_start : TileN_start+TileN] into shared memory
- This is TileN contiguous floats. Very coalesced.
- __syncthreads()
- Each thread computes its outputs using the loaded A and B values.
- If thread computes an 8×8 sub-tile, it needs 8 A-values and 8 B-values from shared memory per k.
- Wait, for each k, the contribution is A[i,k] * B[k,j] to C[i,j].
- For a thread with local sub-tile at (bi, bj), for each local_i in bi..bi+7 and local_j in bj..bj+7: accum[local_i-bi][local_j-bj] += A_smem[local_i] * B_smem[local_j]
- So it needs 8 distinct A-values and 8 distinct B-values.
- 16 shared memory reads per k, 64 FMAs.
- __syncthreads()
Total shared memory: TileM floats + TileN floats. If TileM=128, TileN=128: 512 + 512 = 1KB. Very small! We can use much larger tiles.
What if TileM=256, TileN=256?
- Shared memory: 1KB + 1KB = 2KB
- Threads: 256
- Each thread computes 4x4 = 16 outputs
- Per k: 4 A-reads, 4 B-reads, 16 FMAs. 8 reads for 16 FMAs.
- For K=32: 256 reads per thread, 512 FMAs.
- Total blocks: (32768/256)^2 = 128^2 = 16384 blocks.
What about register pressure? Each thread has 16 accumulators = 16 floats. Plus 4 A temps and 4 B temps. Total ~24 registers. Very comfortable.
Can we go bigger? TileM=512, TileN=256? Then shared memory = 2KB + 1KB = 3KB.
Threads = 512 × 256 / outputs_per_thread? Let's say 256 threads, each computes 2x1 outputs? That's very few FMAs.
Actually, we want each thread to have enough work. Let's target 512-1024 FMAs per thread.
With TileM=128, TileN=256, 256 threads:
- Each thread: 128*256/256 = 128 outputs.
- If organized as 16x16 threads, each thread handles an 8x16 sub-tile.
- Per k: 8 A-reads + 16 B-reads, 128 FMAs.
- Total: 24 reads per k × 32 = 768 reads, 4096 FMAs. Good balance.
Shared memory = 128 + 256 = 384 floats = 1536 bytes. Very small.
Actually, we don't even need shared memory if we can load directly from global memory with caching! On H100, L1 cache is large. If each warp loads its own A and B values into L1 cache, we might not need explicit shared memory.
But without shared memory, we'd have bank conflicts? No, global memory to registers.
Let me think about a version without shared memory:
- For each k in 0..31:
- Warp collaboratively loads A segment into registers? Not straightforward.
- Each thread could load its own A-values from global memory. A[i0+local_i, k]. For 8 values, addresses are base + local_i * 128. Since each thread loads 8 values spaced 128 bytes apart, that's 8 cache lines per thread.
- But these are the same cache lines for all threads in the same row group...
Actually, shared memory is beneficial here because it allows us to share the loaded A and B values among all threads in the block. Without it, each thread would redundantly load the same A and B values multiple times.
With shared memory:
- 128 A-values loaded once, shared by all threads in the block.
- 256 B-values loaded once, shared by all threads.
Total loads from global memory per k = 128 + 256 = 384 floats.
Without shared memory, if each thread loads its own: 256 threads × (8 A + 16 B) = 6144 loads per k.
So shared memory reduces global memory loads by 16x! That's huge.
Given A[:, k] might be in L2 cache (only 128KB), repeated loading isn't terrible, but shared memory is still better.
Let me design the kernel with shared memory.
Actually wait - for A[:, k], if 256 threads redundantly load the same 128 values, they'd all hit L2 cache. But 256 × 128 cache line requests might cause some contention. Better to use shared memory.
Kernel design:
```cuda
template<int TileM, int TileN, int BlockDimX, int BlockDimY>
__global__ void gemm_outer_product(const float* A, const float* B, float* C, int M, int K, int N) {
// A: M x K, B: K x N, C: M x N
int block_i = blockIdx.x * TileM;
int block_j = blockIdx.y * TileN;
int tid = threadIdx.y * BlockDimX + threadIdx.x;
// Each thread computes a sub-tile
constexpr int ItemsPerThreadM = TileM / BlockDimY;
constexpr int ItemsPerThreadN = TileN / BlockDimX;
float accum[ItemsPerThreadM][ItemsPerThreadN] = {0};
__shared__ float A_smem[TileM];
__shared__ float B_smem[TileN];
for (int k = 0; k < K; ++k) {
// Load A column segment
// We need TileM elements, each from A[block_i + idx, k]
// Address: A + (block_i + idx) * K + k
for (int idx = tid; idx < TileM; idx += BlockDimX * BlockDimY) {
A_smem[idx] = A[(block_i + idx) * K + k];
}
// Load B row segment
// B[k, block_j + idx], address: B + k * N + block_j + idx
for (int idx = tid; idx < TileN; idx += BlockDimX * BlockDimY) {
B_smem[idx] = B[k * N + block_j + idx];
}
__syncthreads();
// Compute
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = threadIdx.y * ItemsPerThreadM + ii;
float a = A_smem[i];
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = threadIdx.x * ItemsPerThreadN + jj;
accum[ii][jj] += a * B_smem[j];
}
}
__syncthreads();
}
// Write out
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = block_i + threadIdx.y * ItemsPerThreadM + ii;
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = block_j + threadIdx.x * ItemsPerThreadN + jj;
C[i * N + j] = accum[ii][jj];
}
}
}
```
Wait, C[i * N + j]. For C of shape (32768, 32768), writing 128 rows × 256 columns = 32768 floats = 128KB per block. But each warp writes consecutive j's for a fixed i, so contiguous floats. Then next row. This is coalesced.
Let me think about whether this kernel is better than cuBLAS.
Actually, there's another approach. Since K=32 is very small, we might be able to use tensor cores! Even though the compute is small, tensor cores might still help if we can structure the data well.
But for K=32, we'd need to accumulate across multiple K groups to fill a tensor core MMA instruction. On H100, the smallest MMA might be m16n8k8 or similar. We'd need to tile K into chunks of 8.
Actually, using tensor cores for this problem might be overkill and have setup overhead. Also, since it's memory bound, tensor cores won't help much.
Let me consider another approach: what if we transpose the problem?
C = A @ B where A is (32768, 32), B is (32, 32768).
What if we reshape A as (32, 1024, 32) and do batched operations? Not helpful.
Let's just try implementing the outer product kernel and see.
Tile sizes:
- TileM = 128, TileN = 128, threads = 16×16 = 256. Each thread does 8×8 = 64 outputs.
- Or TileM = 64, TileN = 256, threads = 8×32 = 256. Each thread does 8×8 = 64 outputs.
- Or TileM = 256, TileN = 128, threads = 32×8 = 256. Each thread does 8×16 = 128 outputs.
For the outer product, B[k, block_j:block_j+TileN] needs to be loaded. Larger TileN means more contiguous B data per load, which is good. But each thread does more work.
Actually, let me think about warps. With 32 threads in a warp, if we organize as (8, 4), then each warp covers 64 × TileN/4 or something.
I want B_smem[j] accesses to not have bank conflicts. Since threads in a warp read different j values, and shared memory is banked by 4 bytes, if consecutive threads read consecutive j's, that's a broadcast? No, B_smem[j] for different j values - if j0 and j1 are such that threads read j = base + lane_id, then each thread reads a different bank. No conflicts as long as the stride doesn't cause multiple threads to hit the same bank.
For B_smem with TileN=128, if threadIdx.x goes from 0..15 and each handles 8 j's: thread tx reads B_smem[tx*8 + jj]. For jj=0, threads read B_smem[0], B_smem[8], B_smem[16], etc. With 32 banks, B_smem[0] and B_smem[32] are in the same bank. For 16 threads reading even multiples of 8: banks = 0, 8, 16, 24, 0, 8, ... for tx=0..15, banks are 0, 8, 16, 24, 0, 8, 16, 24, ...
tx=0 reads bank 0, tx=4 reads bank 0, tx=8 reads bank 0, tx=12 reads bank 0. That's 4 threads in same bank! Conflict.
To avoid this, we can pad B_smem. Or better, have threads in a warp read consecutive elements. If thread tx handles a contiguous range of j's, then for a fixed jj in the inner loop, all threads read different elements spaced by ItemsPerThreadN. If ItemsPerThreadN is odd or not a multiple of 32, there might be conflicts.
Actually, a simpler solution: use a 1D thread block and let each thread compute one row of outputs. Then for each k, thread i loads A_smem[i] and all threads in a warp load the same B_smem[j] values but with some pattern...
Let me think about 1D threads computing rows:
- Block size = 256 threads.
- TileM = 256, TileN = something.
- Each thread computes one entire row of the tile (TileN outputs).
- Thread i loads A[i, k] from shared memory (or global).
- Then for j = 0..TileN-1: C[i, j] += A[i, k] * B[k, j]
- Here B[k, j] is read from shared memory by all threads. This is a broadcast! No bank conflicts.
- But B[k, :] is 32768 floats. If TileN = 32768 / num_blocks, say 256, then B_smem is 256 floats = 1KB. 256 threads access 256 different elements. For consecutive j's, consecutive banks. If each thread accesses a unique j, no conflicts!
Wait, if thread i computes C[i, j] for a specific j, and we have 256 threads computing 256 different outputs in the same row... no, then we'd need 256 threads per row and only 1 row per block. That would mean TileM = 1, which is terrible for A loading.
What if each thread computes a full row but for a subset of columns?
- 256 threads, 16 rows (TileM=16), each thread computes 16 columns (TileN=16*16=256).
- For each k: load A_smem[16], load B_smem[256].
- Thread in row r, column group c: computes C[r, c*16 + offset].
- In the compute loop: accum[offset] += A_smem[r] * B_smem[c*16 + offset].
- Threads in same warp: same r? If warp has 32 threads and we have 16 columns per row group, then 2 rows fit in a warp. Or if we organize as 16×16, then a warp of 32 threads covers 2 rows × 16 columns.
- For B_smem access: threads read B_smem[base + lane_id] or similar. Consecutive lanes read consecutive elements = different banks. No conflict.
Hmm, let me try a specific configuration.
TileM=128, TileN=128, Block=(16, 16).
For a warp (32 threads), if arranged as 8 rows × 4 cols of threads:
- Covers 8 × 8 outputs per thread = 64×32 output area? No.
- Thread (ty, tx) with ty in 0..15, tx in 0..15.
- Warp covers thread rows 0..7 and thread cols 0..3? Actually warps are formed in thread ID order. If threads are (tx, ty) with tx varying fastest, then a warp has threads (0..15, ty0) and (0..15, ty0+1) for ty0=0,2,4,...
- So a warp covers 2 thread-rows and 16 thread-cols, which is 2×16 outputs-per-thread in height? No, each thread covers 8×8 outputs. So a warp covers 16 thread-rows? Wait.
With 16×16 threads, thread IDs are 0..255. A warp is 32 consecutive threads. If x varies fastest: threads 0..15 are (0..15, 0), threads 16..31 are (0..15, 1). A warp covers rows y=0 and y=1, all x=0..15.
For compute: thread (x, y) handles output block at (y*8, x*8).
In the compute loop, for each k:
- Thread reads A_smem[y*8 + ii] for ii=0..7.
- Thread reads B_smem[x*8 + jj] for jj=0..7.
For a warp (y=0..1, x=0..15), threads read:
A_smem: for y=0, A_smem[0..7]; for y=1, A_smem[8..15]. These are all different elements, no shared memory bank conflict.
B_smem: for x=0, B_smem[0..7]; x=1, B_smem[8..15]; etc. Also all different. And since the addresses are consecutive, they map to different banks (bank = address % 32). No conflicts!
This is actually quite clean. Let me go with this.
TileM=128, TileN=128, Block=16×16. Each thread does 8×8 outputs = 64 FMAs per k. K=32, so 2048 FMAs per thread. Total blocks = 256×256 = 65536.
Wait, output is 32768×32768. With TileM=TileN=128:
Num blocks in M = 32768 / 128 = 256
Num blocks in N = 32768 / 128 = 256
Total blocks = 65536.
H100 has 132 SMs. Each SM can run multiple blocks. If each block has 256 threads and moderate register usage, occupancy might be decent.
Can we increase tile size to reduce block count and maybe launch overhead? But larger tiles mean more registers. If TileM=256, TileN=256, threads=16×16, each thread does 16×16=256 outputs. That's 256 accumulators, each float. 256 registers just for accum! Plus loop variables, temp a, temp b. Probably ~300 registers. That's over the 255 limit per thread!
So we need to be careful. 8×8=64 outputs is only 64 registers for accum. Very safe.
What about TileM=128, TileN=256, threads=16×16? Each thread does 8×16 = 128 outputs. 128 accum registers. Probably ~160-180 total registers. Still under 255.
Total blocks = 256 × 128 = 32768. Better.
Or TileM=256, TileN=128, threads=32×8 = 256? Each thread does 8×16=128 outputs. Wait, with threads=32×8: ItemsPerThreadM = 256/8 = 32, ItemsPerThreadN = 128/32 = 4. So 32×4 = 128 outputs. But 32 A-reads per k! That's a lot of registers for A temps if we unroll. We can keep A in registers though: 32 floats for A row indices, 128 accumulators = 160 regs. OK.
But loading A_smem: 256 elements for 256 threads, each thread loads 1. Then each thread needs 32 different A-values. That's fine.
Actually, I think we can also try a warp-level approach without shared memory. Since K=32 is small and A columns may be cached...
Let me try another idea. What if we process 32 K values simultaneously using vector operations?
- For an output tile, load A_tile[TileM][32] into shared memory.
- Load B_tile[32][TileN]...but B is (32, 32768), so B_tile[k][j] = B[k*N + j]. For a tile of B, we'd load B[k, j0..j0+TileN-1] for k=0..31. That's 32 rows of length TileN. Total = 32*TileN floats.
- If TileN=128, that's 4096 floats = 16KB of B data.
- A data = TileM*32 floats. If TileM=64, that's 8KB.
- Total shared memory = 24KB. Fine.
- Then each output C[i,j] = dot(A_tile[i][0..31], B_tile[0..31][j]).
Wait, this is the batched dot product approach! And it uses shared memory more heavily but global memory less heavily. Each element of A and B is loaded exactly once from global memory per tile.
If TileM=64, TileN=128:
- C tile size = 64×128 = 8192 outputs.
- With 256 threads, each thread computes 32 outputs.
- Each thread can compute 32 outputs as, say, a 4×8 sub-tile.
- For each output, need a dot product of two 32-float vectors.
- We can load one vector from A_smem (32 floats) and one from B_smem (32 floats) and compute dot product = 32 FMAs.
- But the dot products for different outputs in the sub-tile share the same A_smem row or B_smem column...
Actually, for a thread computing a 1×32 strip: it needs A_smem[i][0..31] (32 floats) and B_smem[0..31][j] for each j in its strip. That's 32 vectors of 32 floats = 1024 floats from B_smem per row.
Shared memory bandwidth could be the bottleneck. For K=32 and 32 outputs, 32 FMAs per output means 1024 FMAs. But loads from shared memory: 32 (for A) + 32×32 (for B) = 1056 float loads. That's barely 1 FMA per load. Not great.
What if each thread computes a 32×1 strip (32 rows, 1 column)?
- Needs B_smem[0..31][j] (32 floats) once.
- Needs A_smem[i][0..31] for each of 32 rows = 1024 floats.
- Same issue.
What if we compute the dot product in a blocked manner? For kk=0..31 in chunks of 4:
- Load A_val[kk] and B_val[kk]
- For all outputs: accum += A_val * B_val
Wait, this is exactly what the outer product approach does! Process one k at a time, outer product. So the outer product is fundamentally the right approach for small K.
For the outer product with TileM=128, TileN=128:
- Global memory reads per block per k: 128 A + 128 B = 256 floats.
- Over K=32: 8192 float reads per block.
- FMAs per block: 128×128×32 = 524,288.
- Read-to-compute ratio: 8192 / 524288 = 1/64.
- Output writes: 128×128 = 16384 floats.
- Total memory traffic per block: 8192 + 16384 = 24576 floats = 96KB.
- Arithmetic intensity: 524288 / (24576 * 4 bytes)? Wait, 524288 FMAs / (24576 × 4 bytes) = 524288 / 98304 = 5.3 FLOPs/byte.
- This is memory bound (turning point at ~20 FLOPs/byte).
So the bottleneck is writing C. We need to maximize C write bandwidth. Writing 128×128 floats per block, if the writes are coalesced, that's good.
Can we reduce the number of blocks to reduce C write overhead? Well, larger tiles would help, but register limits apply.
What about using vectorized stores for C? Writing float4 or float2? That might help slightly with instruction efficiency.
Actually, for output writes, each thread writes 64 floats (8×8). Can we write 8 float4s per thread? Each float4 is 16 bytes. 8 float4s = 128 bytes. Thread writes at addresses: C + i*N + j, C + i*N + j+1, etc. These are contiguous for consecutive j. So we can vectorize to float4!
Let's think about vectorization:
- Thread computes (8, 8) = 64 floats.
- We can store as 16 float4s. For each ii in 0..7: store 2 float4s at C[(i0+ii)*N + j0 + jj*4].
- This requires N to be divisible by 4. For N=32768, yes!
Vectorized loads for A and B:
- Loading A: A[(block_i + idx) * K + k]. For 128 elements strided by K=128 bytes. Float4 won't help since the stride is 128 bytes.
- Loading B: B[k * N + block_j + idx]. 128 contiguous floats. Can load as 32 float4s!
So B loading can be vectorized by 4x.
What about A loading? If K=32, A[(block_i + idx) * 32 + k]. The stride between A[i,k] and A[i+1,k] is 32*4=128 bytes. So we can load 4 consecutive rows as a float4? Let's check:
A[(block_i + idx) * 32 + k] is a float.
If thread idx loads float4 starting at A[(block_i + idx*4) * 32 + k]:
Addresses: base, base+128, base+256, base+384. These are 128 bytes apart. A float4 requires 16 contiguous bytes. So this doesn't work.
But we could load with a different pattern. Have 32 threads each load 4 floats from A using a strided pattern? No, float4 requires contiguous addresses.
We could just load floats individually. For 128 elements, 128 loads. It's fine.
Actually, let me reconsider if we can pack A differently in memory. A has shape (32768, 32). If we transpose A first (4MB), then rows become contiguous in what was the K dimension. But then for the outer product approach we'd need A[:, k] which would be columns of the original... no, same problem.
Wait, I have an idea. Instead of loading A[:, k] for each k, what if we load A rows and keep them in registers/shared memory, then iterate k in a way that accesses B row by row? That's the batched dot product approach, and it requires loading A once and B row by row.
Actually let's go back to the dot product with transposed B. If we transpose both A and B so that K becomes the leading dimension:
A has shape (32768, 32). If we view A as (32, 32768) with appropriate strides, then the dot product C[i,j] = sum_k A[k, i] * B[k, j]. This is a standard matrix multiplication of A.T (32, 32768) @ B (32, 32768). No wait, A.T is (32, 32768), B is (32, 32768). For matmul, we'd need A.T.T @ B = A @ B. That's (32768, 32) @ (32, 32768).
If we store A as (32, 32768) column-major? Not how PyTorch works.
Actually, what about grouping the K dimension into warp-level operations? Since K=32 and a warp has 32 threads, we could have each warp compute an output tile where each thread handles one k value!
Consider: a warp of 32 threads collaborates to compute a TileM × TileN output tile.
- For each i in TileM and j in TileN:
- Thread k loads A[i, k] and B[k, j]
- All threads do a warp-level reduction (sum) to get C[i,j]
But this is terrible because we'd need 32 loads for each output element. Very inefficient.
Alternatively, for a fixed i and j:
- Thread k computes A[i,k] * B[k,j] for its assigned k's.
- Then warp sum to get C[i,j].
If each thread handles 1 k (K=32, 32 threads), then each thread does 1 multiply per output. To compute one output, 32 threads load 32 A values and 32 B values. That's 64 loads for 1 output. Terrible.
So the outer product approach is much better: 1 load of A[i,k] is reused for TileN outputs, and 1 load of B[k,j] is reused for TileM outputs.
I think the outer product kernel is the way to go. Let me refine it.
To maximize performance:
1. Vectorized B loads (float4)
2. Coalesced C writes (float4)
3. Good occupancy with small shared memory
4. Possibly use large TileN to amortize A loads
Let me also consider warp-level mma instructions. H100 has tensor cores. The smallest MMA on Hopper is probably m16n8k8 for fp32 accumulation with fp16 input, but for fp32 input, it's m16n8k8 with TF32 or something. But for fp32 accumulation of fp32, it might be m16n8k8 -> 2 accumulations? I need to check.
Actually, for K=32, even with tensor cores, we'd be limited by memory bandwidth. Tensor cores do 8 multiplies per thread per instruction, but if we're memory bound, it won't help. The bottleneck is writing 4GB.
What about using fewer threads per block but larger tiles? Say 128 threads, TileM=128, TileN=256. Then each thread does 1×2 or something.
Actually, here's another idea. What if we have each thread compute a full 32-dot product using a tree reduction in shared memory? Or better, since K=32, we can have each thread compute 4 dot products with K/8=4 using float4 vector operations, then accumulate.
Actually, I should check what the actual PyTorch performance is, but I can't run the reference. I need to make my best guess.
Let me also think about an alternative: using CUTLASS or a batched GEMV. Since K=32 is so small, this is essentially 32768 independent dot products of length 32 (for each row of A with rows of B.T). Or equivalently, 32768 independent outer products? No.
Actually, C = A @ B can be viewed as:
For each row i of A: C[i,:] = A[i,:] @ B
Where A[i,:] is 1×32 and B is 32×32768.
This is 32768 independent GEMV operations! Each GEMV is a vector (32) times a matrix (32×32768) to produce a row vector (32768).
But all these GEMVs share the same B matrix. If we process them in parallel, each thread block handles a set of rows.
For GEMV: C_row[j] = sum_k A_row[k] * B[k,j]
This is exactly the outer product approach! For each k, we broadcast A_row[k] to all j and multiply with B[k,:].
So our kernel is doing multiple GEMVs in parallel.
Can we make each GEMV more efficient? With K=32, a dot product of two 32-element vectors. If we vectorize to float4, each dot product needs 8 multiply-adds of float4s. But B[k,:] is strided by N=32768 for consecutive k. So we can't use float4 across k.
However, we can rearrange B in memory by explicit transpose. If B.T is (32768, 32), then B.T[j,:] is contiguous. Then each output is dot(A[i,:], B.T[j,:]). Both are 32 contiguous floats!
Now we can vectorize the dot product with float4:
- Load A[i, 0..3], A[i, 4..7], etc. as float4
- Load B.T[j, 0..3], B.T[j, 4..7], etc. as float4
- 8 multiply-adds of float4 = 32 scalar FMAs
- Sum the 8 float4 partials into a scalar result
This is very efficient! The core computation per output is just 8 vector FMAs.
But to do this, we need B transposed. We can either:
1. Pre-transpose B in a separate kernel (4MB, very fast)
2. Have the main kernel transpose B on the fly into shared memory
Option 1 is simple. Launch a transpose kernel, then do the dot-product GEMM.
Transpose: B (32, 32768) -> B_T (32768, 32).
Each output element B_T[j,k] = B[k,j].
This can be done by blocks that load a chunk of B rows into shared memory and write transposed. But actually for this shape, a simple kernel where each thread handles one column of B (which becomes one row of B_T) is trivial. Or we can just call torch.transpose or use PyTorch's transpose (which is just a view with stride change, not actually moving data!).
Wait! PyTorch transpose is just a view. B.T would have shape (32768, 32) with strides (1, 32768). When we access B.T[j, k], the memory access pattern is still strided!
We need a contiguous transpose. So we'd need to actually copy the data into a new tensor with different layout. That's `B.contiguous()` after transpose, or `permute + contiguous`.
So a custom transpose kernel that writes B to B_T contiguously is useful.
But then the total work is:
1. Transpose B: 4MB read + 4MB write = 8MB
2. GEMM: 4MB read A + 4MB read B_T + 4GB write C = ~4GB
Total: still ~4GB, plus 8MB overhead.
This might be faster than the outer product approach because:
- Outer product reads A 256 times (1GB) and B 256 times (1GB) for the tiled approach, plus 4GB output = 6GB total.
- Dot product reads A once (4MB) and B_T once (4MB), plus 4GB output = ~4GB total.
Wait, for the dot product approach, if each block loads A rows and B_T rows, how many times are they read?
- A has 32768 rows, each 32 floats = 128 bytes. Total 4MB. Each row is read by all blocks that compute a column tile. If we tile N dimension into 256 tiles of 128 columns each, each A row is read 256 times. So A is read 1GB.
- Similarly for B_T rows.
Oh! In both approaches, the inputs are read multiple times because each block only computes a subset of outputs and needs all the rows that contribute to its outputs.
To read A only once, we'd need each block to compute a full row of C (all N outputs). Then each block needs only one A row. With 32768 rows, that's 32768 blocks. Each block computes 32768 outputs (one row of C).
- Block loads A[i, 0..31] = 128 bytes
- Block loads all of B_T (or all needed rows): but B_T has 32768 rows, each 32 floats. For a row of C, we need ALL rows of B_T! That's 4MB per block.
- Block computes 32768 dot products, each of 32 FMAs = 1M FMAs.
- Block writes 32768 floats = 128KB.
This is actually very interesting! A block computing a full row of C:
- Loads 128 bytes of A (once, stays in registers)
- Loads 4MB of B_T... but does it need all of it at once?
- B_T is (32768, 32). To compute C[i,j] = dot(A[i,:], B_T[j,:]), block needs all B_T rows.
- Block can iterate over chunks of B_T rows. Say chunk size = 1024 rows.
- Load 1024 × 32 floats = 128KB of B_T into shared memory
- Compute 1024 dot products using A[i,:] in registers
- Write result
- Total iterations = 32768 / 1024 = 32.
- Each iteration: load 128KB, do 1024 × 32 = 32768 FMAs, write 1024 floats = 4KB.
- Arithmetic intensity: 65536 FLOPs / (128KB + 4KB) = 65536 / 132KB = 0.5 FLOPs/byte. Very memory bound on reads.
Hmm, B_T reading is 4MB per row of C. With 32768 rows, total B_T reads = 128GB! That's way too much.
But wait - if blocks for different rows of C run concurrently, they all read the same B_T data. With L2 cache of 50MB, B_T is only 4MB. So once a chunk of B_T is loaded by one block, subsequent blocks can hit L2 cache!
So if we schedule all row-blocks to execute and they share B_T through L2 cache, total B_T traffic from DRAM is only 4MB (if perfectly cached). But L2 cache bandwidth is high too.
Actually, the key insight: if we can get all blocks for different rows to share the same B_T data in L2 cache, the total read bandwidth is minimized.
But scheduling 32768 blocks with all of them competing for L2 might not be efficient either.
Alternatively, we can process multiple rows of A in one block. Say TileM=32 blocks process 32 rows at a time.
- Block loads A[32 rows × 32 floats] = 4KB into shared memory or registers.
- Block loads a chunk of B_T [1024 rows × 32 floats] = 128KB.
- Computes 32 × 1024 = 32768 outputs.
- Thread assignment: 256 threads, each computes 128 outputs.
- For B_T chunk, iterate j in 0..1023, thread does dot products for its assigned (i,j) pairs.
This reduces block count to 32768/32 × 32768/1024 = 1024 × 32 = 32768 blocks. Still a lot.
Let me go back to the outer product approach with small tiles. It has the advantage of very regular access patterns and small shared memory. And the redundant reads of A and B are well-cached because:
- A[:, k] is 128KB. 256 blocks in a row of the grid read the same A[:, k]. The 128KB easily fits in L2 and can even fit in L1 (50KB? No, L1 is smaller). But L2 cache (50MB) can hold all A[:, k] for all k simultaneously (4MB).
- B[k, :] segment per block is only 512 bytes (for TileN=128). Not an issue.
So the outer product approach might actually have close to optimal memory traffic when accounting for caching.
Actually let me think about this more carefully. With the outer product approach and TileM=128, TileN=128:
- The grid is 256×256 = 65536 blocks.
- All blocks in row `bi` of the block grid (i.e., same blockIdx.x, varying blockIdx.y) read the same A data: for each k, they read A[bi*128 : (bi+1)*128, k].
- There are 256 such blocks per row. All 256 read the same 128 A floats for each k. The 128 floats = 512 bytes. With 256 blocks, that's 256 × 512 = 128KB of A read traffic per k per row of blocks. But with L2 caching, the actual DRAM traffic is 512 bytes (the first block reads it, subsequent blocks hit L2).
- For the 256 rows of blocks (blockIdx.x = 0..255), total A reads from DRAM = 256 × 32 × 512 bytes = 4MB. Exactly the size of A! Perfect cache utilization.
- Similarly for B: all blocks in column `bj` read the same B segment for each k. 256 blocks per column, each reads B[k, bj*128:(bj+1)*128] = 512 bytes. Total B DRAM traffic = 256 × 32 × 512 = 4MB. Exactly the size of B!
So with perfect L2 caching, total input DRAM traffic is just 8MB! Output is 4GB. Total = ~4GB.
This is actually optimal.
Now, can a custom kernel achieve this and beat cuBLAS? cuBLAS likely also achieves close to optimal caching. The difference might be in:
1. Kernel launch overhead
2. Thread block scheduling
3. Register usage and occupancy
4. Vectorization of stores
5. Avoiding unnecessary transpose workspace
For cuBLAS, it might create a workspace for transposing B internally, or use less optimal algorithms. A custom kernel tailored to K=32 might have an edge.
Let me also consider: what if we use the NVIDIA CUTLASS or similar? But we need a custom kernel.
Another idea: since K=32 and A[:, k] reads are strided by exactly 128 bytes (cache line size), we can optimize this by having each thread load a full cache line using float4, but from different rows! Actually, float4 is 16 bytes. A cache line is 128 bytes = 8 float4s.
If thread tid loads a float4 from A[(block_i + tid*4) * 32 + k], that's 4 floats from every 4th row. Not helpful for the outer product.
What if we process 4 k values at once? For k=0,1,2,3:
- A[:, k+kk] for kk=0,1,2,3: four columns of A. For a row i, these are stored contiguously: A[i*32 + k], A[i*32 + k+1], A[i*32 + k+2], A[i*32 + k+3]. So if we load 4 columns at a time, each row has 4 contiguous floats!
- Thread tid can load float4 from A[(block_i + tid) * 32 + k]. This is a single float4 (4 contiguous floats) from one row.
- This is perfectly coalesced! Because thread tid reads address base + tid * 128 + k * 4. For k fixed and tid=0..31, consecutive threads read consecutive cache lines, and within each cache line, they read a single float4. This is a very good pattern.
So if we process K in chunks of 4, for each chunk kk in 0..7:
- Load A_sub[TileM][4]: each thread loads float4 A[(block_i + tid) * 32 + kk*4 + offset]
- Load B_sub[4][TileN]: for each of the 4 k values, load B[kk*4 + offset, :] segment. This is 4 × TileN contiguous floats.
Actually, B[k, :] for k=kk*4, kk*4+1, kk*4+2, kk*4+3. These are 4 separate rows, each of length TileN. Total 4×TileN floats. We can load them as 4 separate segments.
Then for each i, j in the tile:
C[i,j] += A[i, kk*4+0] * B[kk*4+0, j]
+ A[i, kk*4+1] * B[kk*4+1, j]
+ A[i, kk*4+2] * B[kk*4+2, j]
+ A[i, kk*4+3] * B[kk*4+3, j]
We can compute this using 4 independent outer products in one iteration! This increases the compute-to-synchronization ratio (fewer __syncthreads()) and allows vectorized loads of A.
With K=32, we'd have 8 iterations instead of 32. Each iteration does 4× more FMAs per thread per syncthreads. This might help by reducing synchronization overhead and improving instruction-level parallelism.
Let me design this:
```cuda
for (int kk = 0; kk < K; kk += 4) {
// Load A[TileM][4]
// Each thread loads 4 floats as float4 from its assigned row(s)
for (int idx = tid; idx < TileM; idx += BlockSize) {
A_smem[idx][0..3] = A[(block_i + idx) * K + kk] as float4;
}
// Load B[4][TileN]
for (int idx = tid; idx < 4 * TileN; idx += BlockSize) {
int k_off = idx / TileN;
int j_off = idx % TileN;
B_smem[k_off][j_off] = B[(kk + k_off) * N + block_j + j_off];
}
__syncthreads();
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = threadIdx.y * ItemsPerThreadM + ii;
float4 a = A_smem[i]; // 4 floats
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = threadIdx.x * ItemsPerThreadN + jj;
// Compute 4 FMAs
accum[ii][jj] += a.x * B_smem[0][j] + a.y * B_smem[1][j] + a.z * B_smem[2][j] + a.w * B_smem[3][j];
}
}
__syncthreads();
}
```
Wait, this would require A_smem to be an array of float4. Or we can store as `float A_smem[TileM][4]`. But shared memory indexing would be `A_smem[i * 4 + k]`.
For B_smem, we can do `float B_smem[4][TileN]`. Bank conflicts: thread (x, y) for fixed x, y varies. But in a warp, if threads access B_smem[k][j] for different j, no conflict as long as j varies. But here all threads access B_smem[0][j], B_smem[1][j], etc. in the inner loop with fixed j per thread per kk iteration. Since each thread accesses a unique j (or small set of j's with unroll), there should be no bank conflicts.
Let me check: In a warp, threads have different threadIdx.x. If each handles 8 consecutive j's via unroll, for a specific B_smem[k][j_base + offset], all threads access different j. Bank = j % 32. Since j varies across threads, banks vary. No conflict!
This grouped-k approach is very promising. It halves the number of syncthreads and allows vectorized A loads.
Let's see how much shared memory we need:
- A_smem: TileM × 4 floats = TileM × 16 bytes. For TileM=128: 2KB.
- B_smem: 4 × TileN floats = TileN × 16 bytes. For TileN=128: 2KB.
Total: 4KB. Still very small.
For C writes, we can vectorize to float4. Each thread writes ItemsPerThreadM × ItemsPerThreadN / 4 float4s.
One concern: the inner loop now has a long dependency chain for the 4 FMAs. But modern GPUs can pipeline this well. We can unroll or interleave.
Actually, let me compute FMAs per thread for K=32:
With TileM=128, TileN=128, threads=16×16=256, ungrouped: 32 iterations × 64 FMAs = 2048 FMAs.
Grouped by 4: 8 iterations × (64 × 4) = 2048 FMAs. Same total compute, fewer syncs.
This should definitely be better. Let me code this up.
Another optimization: process K in chunks of 8 using float2 or even all 32? We can't load 32 floats as a single vector easily, but K=32 means A row is 128 bytes = exactly one cache line! So A[i, 0..31] is one cache line. We could load it with `ld.global.ca.L2::128B` or similar. But float4 loads of 8 consecutive floats (32 bytes) might be good enough.
Actually, A row is 32 contiguous floats = 128 bytes. We can load it as two float4 loads. But for the outer product, we can load the full row of A into registers and use it for the entire tile computation!
Wait, if each thread computes a sub-tile and loads its A row data into registers, then for each k iteration it already has A[i,k] in a register. But for multiple k's, it can load the full row vector and use it over the entire tile computation without reloading.
Let me think: for a thread block computing TileM × TileN outputs, with K=32:
- Load A[block_i : block_i+TileM, 0:32] into shared memory. That's TileM × 32 floats = TileM × 128 bytes.
- For TileM=128: 16KB.
- For TileM=64: 8KB.
- Load B[k, block_j:block_j+TileN] into shared memory in chunks? But B has 32 rows, each TileN floats. Total B data needed = 32 × TileN floats.
- For TileN=128: 16KB.
- Total shared memory = 32KB for TileM=64, TileN=128. Or 48KB for TileM=128, TileN=128. That's getting close to limits.
If we load all B rows for the tile into shared memory at once, we don't need to iterate over K in the outer product sense. Instead:
- A_smem holds TileM × 32 floats
- B_smem holds 32 × TileN floats
- Then compute: for each i, j: accum = dot(A_smem[i][0:32], B_smem[0:32][j])
But computing dot products directly means we need to access 32 elements per output. With accumulators for TileM_sub × TileN_sub outputs, and accessing A_smem[i][k] and B_smem[k][j] for each k, we have shared memory bandwidth concerns.
For the dot product computation:
- Thread computes 8×8 outputs.
- For k=0..31:
- Read A_smem[i][k] for each of 8 rows = 8 reads
- Read B_smem[k][j] for each of 8 cols = 8 reads
- Do 8×8 = 64 FMAs (with broadcasting/batching)
Wait, for a fixed k, A_smem[i0+ii][k] is read and B_smem[k][j0+jj] is read. If each thread handles 8 rows and 8 columns, it needs 8 A values and 8 B values. 16 reads for 64 FMAs. Same as the outer product approach!
But the difference is that all A and B data is preloaded in shared memory. So there are no global memory loads in the compute loop, only shared memory loads. The outer product approach also had only shared memory loads in the compute loop after loading global -> shared each k.
With the preloaded approach:
- One global load phase: load A_smem (TileM × 32 floats) and B_smem (32 × TileN floats)
- Compute phase: 32 iterations over k
- Write phase: write C
Outer product grouped approach:
- 8 global load phases (for K=32, group=4)
- Compute phase within each group
- Write phase
The preloaded approach might be slightly better because it has fewer global load transactions (no redundant syncthreads overhead) and the global loads are more consolidated. But shared memory is larger.
For TileM=64, TileN=256:
- A_smem: 64 × 32 = 2048 floats = 8KB
- B_smem: 32 × 256 = 8192 floats = 32KB
- Total: 40KB. Under the 48KB limit (or 228KB limit on H100 for certain configs). This works!
- Threads: 256. Each thread computes 64×256/256 = 64 outputs. Or say 8×8 = 64.
- For each k in 0..31: 8 A-reads + 8 B-reads, 64 FMAs.
But the A_smem load pattern: A[block_i + idx, 0..31] for idx=0..63. Each row is 32 floats = 128 bytes. 64 rows = 8KB. Can be loaded coalesced with float4.
B_smem load pattern: B[k, block_j : block_j+256] for k=0..31. Each row is 256 contiguous floats = 1KB. 32 rows = 32KB. Very coalesced.
This preloaded approach with TileM=64, TileN=256 seems very good!
Actually, we can even do TileM=64, TileN=512:
- A_smem: 8KB
- B_smem: 32 × 512 = 16K floats = 64KB
- Total: 72KB. Still under 228KB!
- Threads: say 512? Or 256 threads each handling more outputs. With 256 threads, each does 64×512/256 = 128 outputs. Say 16×8.
- But 512 threads gives more parallelism per block.
However, larger tiles mean fewer blocks. For M=32768, N=32768:
- TileM=64, TileN=256: blocks = 512 × 128 = 65536
- TileM=64, TileN=512: blocks = 512 × 64 = 32768
- TileM=128, TileN=256: A_smem=16KB, B_smem=32KB, total=48KB. Blocks = 256 × 128 = 32768.
TileM=128, TileN=256 is interesting:
- 128 rows × 32 cols of A = 4096 floats = 16KB
- 32 rows × 256 cols of B = 8192 floats = 32KB
- Total shared memory = 48KB.
- Threads: 256 or 512. With 512 threads, each handles 128×256/512 = 64 outputs. 8×8.
- Block count: 256 × 128 = 32768. Still plenty.
Can we do TileM=128, TileN=512?
- A_smem: 16KB
- B_smem: 32 × 512 = 16K floats = 64KB
- Total: 80KB. Still OK!
- Threads: 512, each handles 128 outputs. 16×8 or similar.
- Blocks: 256 × 64 = 16384.
TileM=256, TileN=256:
- A_smem: 32KB
- B_smem: 32KB
- Total: 64KB.
- Threads: 1024? But max threads per block is 1024. With 1024 threads, each handles 64 outputs.
- Blocks: 128 × 128 = 16384.
What about TileM=256, TileN=512?
- A_smem: 64KB
- B_smem: 64KB
- Total: 128KB. Possible on H100!
- Threads: 1024, each handles 128 outputs.
- Blocks: 128 × 64 = 8192.
Actually, we need to be careful about register limits. More outputs per thread = more accumulators. And we also have loop variables. If a thread handles 16×16 = 256 outputs, that's 256 accum registers. Too many.
For 128 outputs (16×8), that's 128 accum registers. Plus other registers. Might be over 255 or close to it.
Let's target 64 outputs per thread = 64 accum registers. Very safe.
With TileM=128, TileN=256, threads=512:
- 128×256/512 = 64 outputs per thread.
- A_smem: 16KB, B_smem: 32KB, total: 48KB.
- 32768 blocks.
- Each thread: 64 accum + some temps + loop counters. Probably ~90 registers. Good occupancy.
With TileM=128, TileN=512, threads=1024:
- 128×512/1024 = 64 outputs per thread.
- A_smem: 16KB, B_smem: 64KB, total: 80KB.
- 16384 blocks.
- 1024 threads per block. On H100, max threads per SM is 1536 or 2048. So at most 1-2 blocks per SM. Might limit occupancy, but 1024 threads is plenty of warps (32 warps). Good occupancy actually.
Hmm, 1024 threads means at most 1 block per SM (since 1024 < 1536 but > 1024? Actually H100 can run 32 blocks per SM concurrently, but limited by thread count. 2048 threads max per SM. So 1024 threads means max 2 blocks per SM. And shared memory of 80KB per block, with 228KB per SM, means max 2 blocks per SM. So occupancy would be 2 blocks × 1024 threads = 2048 threads per SM. That's actually full occupancy!
But register file is 64K per SM. If each thread uses ~100 registers, that's 102400 per block. 2 blocks = 204800 < 256K (64K × 4 bytes? Wait, 64K 32-bit registers per SM, but each thread can use up to 255). 1024 threads × 100 regs = 102400. 2 blocks = 204800. Under 64K? No, 64K = 65536 registers per SM. 204800 > 65536!
Wait, I misread. The spec says: "Register File Size: 64K 32-bit registers per SM". That's 65536 registers total per SM!
If each thread uses 100 registers, that's 1024 × 100 = 102400 per block. Already exceeds 65536!
Oh wow, I completely misread that. 64K registers per SM is tiny! This means occupancy is severely limited by registers.
For 65536 registers per SM and max 2048 threads per SM, average registers per thread = 32. But maximum is 255. If each thread uses 96 registers, then per block with 1024 threads: 98304 registers. Exceeds SM limit! So only 1 block per SM, and it would need to use ≤64 registers per thread if we want full 1024 threads? No, if 1 block has 1024 threads, it can use up to 65536/1024 = 64 registers per thread on average.
Actually, that's not how it works. If a block requests 1024 threads and uses 96 registers each, the hardware will try to allocate 98304 registers. But there's only 65536. So either:
1. The kernel won't launch (compile-time error if known), or
2. The hardware reduces occupancy by not running as many threads.
Wait, actually on NVIDIA GPUs, the register file is partitioned. The maximum registers per thread is 255. The compiler will allocate what's needed. At runtime, if a block needs more registers than available, it can't run. But actually the compiler spills certain decisions. Let's think...
No, if a kernel is compiled with, say, 64 registers per thread, a 1024-thread block needs 65536 registers, which exactly fills the SM. So only 1 block can run. With 96 registers, 1024 threads would need 98304, which is impossible. The hardware would run fewer threads? No, it runs fewer blocks per SM. A 1024-thread block with 96 registers cannot run at all because it needs more than an SM's register file. Wait, that can't be right because many kernels use 64+ registers with 1024 threads.
Let me reconsider: "64K 32-bit registers per SM". Is that 64K = 65536, or does "K" mean something else? Looking at the spec more carefully: "64K 32-bit registers per SM". I think this might be a typo or I might be misreading. Actually, Ampere SMs have 64K 32-bit registers per SM? No, Ampere has 65,536 registers per SM. And Hopper also has 64K? That seems very small. Let me check: Ampere GA100 has 64K registers per SM? That doesn't sound right. I thought GA102 had 16,384 registers per SM or something.
Actually, I think "64K" might mean 64 × 1024 = 65536. And yes, Ampere and Hopper SMs indeed have 64K 32-bit registers. With 128 CUDA cores per SM, that's 512 registers per core. It seems small but that's actually the spec.
But then how do 1024-thread blocks work with >64 registers? They might not! Let's check: if each thread uses 64 registers, 1024 threads need 65536. Exactly fills the SM. So a 1024-thread block with 64 registers per thread uses the entire register file. Only 1 concurrent block.
If a thread uses 128 registers, then only 512 threads can run concurrently per SM. But a block might still have 1024 threads - the SM just can't run the full block at once? Wait, actually on GPUs, a block is not split across SMs. The entire block runs on one SM. So if a 1024-thread block needs 128 registers per thread, that's 131072 registers. But the SM only has 65536. The kernel would fail to launch!
Hmm, but that contradicts common practice. Many kernels use 128 registers with 1024 threads. What's going on?
Oh! I think I misread the architecture. Let me look again. The spec says "Register File Size: 64K 32-bit registers per SM". But maybe this is per SM partition? On Hopper, each SM might have multiple partitions. Or maybe "64K" stands for something else.
Actually, looking up Ampere: each SM has 64K 32-bit registers. A GA102 SM can indeed run 1536 threads concurrently. If each thread uses 128 registers, max threads would be 65536/128 = 512. But the hardware reports occupancy limited by registers.
For a block of 1024 threads, max registers per thread is 65536/1024 = 64. If the kernel is compiled with more than 64 registers, the kernel would not be able to run a 1024-thread block.
Wait, but what about 32-register limit? No, max is 255. But practical max for a large block is 64.
OK so we need to be very careful. Let's target ≤64 registers per thread for 1024-thread blocks, or use smaller blocks (512 threads → max 128 registers per thread, or 256 threads → max 256 registers per thread).
For our kernel, let's use 256 threads per block. Then we can comfortably use up to 255 registers. With 64 outputs per thread (64 accum registers) and some loads, total might be ~80-100. That's fine with 256 threads (max ~256 registers/thread).
Wait, with 256 threads per block, we can have up to 255 registers. And 256 threads × 255 = 65280, just under 65536. So a 256-thread block can fully utilize the register file if needed.
OK so smaller blocks give us more register breathing room. And we can have multiple blocks per SM. With 256 threads × 80 regs = 20480 regs. 65536 / 20480 ≈ 3 blocks per SM. With 32 max blocks per SM, occupancy limited by registers to 3 blocks = 768 threads out of 2048 max. That's 37.5% occupancy. Not terrible, but not great.
If we can keep registers to ~64 per thread: 256 × 64 = 16384. 65536 / 16384 = 4 blocks per SM. 1024 threads = 50% occupancy.
I should write the kernel with compiler hints to manage register usage, though nvcc mostly handles this.
Let me go with TileM=128, TileN=128, BlockSize=256 (16×16). Each thread handles 8×8 = 64 outputs.
- A_smem: 128 × 32 floats = 16KB.
- B_smem: 32 × 128 floats = 16KB.
- Total shared: 32KB. Under the 48KB limit.
- 256 threads per block.
- Block count: 256 × 256 = 65536.
In the compute loop, each thread does for each k:
a = A_smem[i0 + ii][k] -- wait, A_smem is 2D: A_smem[i][k]
b = B_smem[k][j0 + jj]
If A_smem is stored as `float A_smem[128][32]`, then `A_smem[i][k]` is row-major. Row i is 128 bytes. With 32 banks, accesses in a warp might conflict if not careful.
Let me think about the memory layout for A_smem. We store A block as `A[i][k]` where i=0..127 and k=0..31. In shared memory, we want consecutive threads accessing consecutive elements. For the compute loop, for a fixed k and varying i, we want A_smem[i][k] for consecutive i. In row-major layout `A_smem[i][k]`, the address of A_smem[i][k] is i*32 + k. For consecutive i, addresses differ by 32. Bank = (i*32 + k) % 32 = k % 32. So ALL threads in a warp accessing the same k but different i will hit the SAME bank! Massive bank conflicts!
To avoid this, we need to transpose A in shared memory. Store A_smem as `A_smem[k][i]`. Then A_smem[k][i] for fixed k and varying i has address k*128 + i. Banks = i % 32. Consecutive i means consecutive banks. No conflict!
Similarly for B_smem, we store as `B_smem[k][j]`. For fixed k and varying j, address is k*TileN + j. Banks = j % 32. If consecutive threads access consecutive j, no conflict.
Wait, but in our compute loop, thread (tx, ty) accesses A_smem[k][i0 + ty*8 + ii] and B_smem[k][j0 + tx*8 + jj]. For a warp of threads with varying tx and ty, the i values are not consecutive in the warp...
If a warp consists of threads (tx=0..15, ty=0) and (tx=0..15, ty=1), then the i values accessed for A_smem are from two groups: i = ty*8 + ii. For ty=0: i=0..7; ty=1: i=8..15. So 16 consecutive i values! Banks = i % 32, all different for i=0..15. No conflict.
For B_smem, j values are tx*8 + jj. For tx=0..15, j spans 0..127 (spread across multiple values). For a specific jj, j = tx*8 + jj. Banks = j % 32. For jj=0: j = 0, 8, 16, ..., 120. Banks = 0, 8, 16, 24, 0, 8, ...
tx=0: bank 0, tx=4: bank 0, tx=8: bank 0, tx=12: bank 0. 4 threads in bank 0! Conflict!
To avoid B_smem bank conflicts, we can transpose B_smem as well! Store as B_smem[j][k]. Then for fixed j, B_smem[j][k] for varying k... but in our compute loop, k is fixed (outer loop over k) and j varies.
Hmm, let me think about the compute pattern differently. In the outer product approach (preloaded), the outer loop is k=0..31:
for each thread:
for each ii, jj:
C[ii][jj] += A_smem[k][i0+ii] * B_smem[k][j0+jj]
For a warp, different threads have different (i0+ii, j0+jj) pairs. For fixed k, they all read A_smem[k][i] and B_smem[k][j].
If we store A_smem as [k][i] and B_smem as [k][j], and we have 32 threads accessing A_smem[k][i] with i values from a set like {0,1,2,3,4,5,6,7, 8,9,10,11,12,13,14,15, ...}, then banks are i % 32. If i values are consecutive or spread with gaps less than 32, there's no conflict. With our 16×16 thread layout and ItemsPerThreadM=8, a warp covers i=0..15 (first 16 threads) and i=64..79 (next 16 threads, if ty=8 for the second half). Wait, warps are formed by thread ID order. If threadIdx = ty*16 + tx, then:
- tid 0..15: ty=0, tx=0..15, i=0..7 (for a specific ii)
- tid 16..31: ty=1, tx=0..15, i=8..15
So in a warp, i values are 0..15. Banks 0..15. No conflict. Great!
For B_smem[k][j], j values in warp: tx=0..15, j = tx*8 + jj for a specific jj.
If jj=0: j = 0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120.
Banks = j % 32 = 0, 8, 16, 24, 0, 8, 16, 24, 0, 8, 16, 24, 0, 8, 16, 24.
Only 4 unique banks, 4 threads per bank. 4-way conflict!
To fix this, we can:
1. Pad B_smem rows: `float B_smem[32][TileN + 4]`. Then j=tx*8 + jj, address = k*(TileN+4) + j. Banks = (k*(TileN+4) + j) % 32. If TileN+4=132, then k*132 % 32... this might still cause conflicts.
2. Transpose B_smem to [j][k]. Then for fixed k and varying j, we access B_smem[j][k]. Address = j*32 + k. Banks = (j*32 + k) % 32 = k % 32. ALL threads hit the same bank! That's worse.
3. Change thread layout so that consecutive threads access consecutive j values. If each thread handles 8 consecutive j's and we iterate jj from 0 to 7, then for jj=0 all threads access j_start + lane_id. If j_start is a multiple of TileN/BlockSize...
Actually, let's reorganize the thread assignment. Instead of thread (tx, ty) handling j = tx*8 + jj, what if we assign threads to outputs in a 1D fashion in the j dimension?
For TileN=128 and 256 threads, we could have 256 threads each handling 1 output j (but then they'd all need all 128 i values, which means each thread computes 128 outputs, one per i... wait, that means each thread handles an entire column! No, each thread handles a full column of 128 outputs. Then for each k, thread reads A_smem[k][i] for all i and B_smem[k][j] once. 128 A-reads + 1 B-read per k. 129 reads for 128 FMAs. Not great for shared mem bandwidth.
What if threads in a warp handle consecutive j's? With 32 threads in a warp, and each thread handling 4 j's (e.g., ItemsPerThreadN=4), then for jj=0, j = base_j + lane_id. Banks = (base_j + lane_id) % 32 = lane_id % 32. All different banks! No conflict.
But the outer loop over k means each thread only accesses one B_smem[k][j] per inner iteration. If we have a 2D thread block and handle the j assignment differently...
Let me try: Block = 32×8 = 256 threads. tx in 0..31, ty in 0..7.
TileM=128, TileN=128. Each thread handles 4×16 or 16×4 outputs.
If ItemsPerThreadM=16, ItemsPerThreadN=4:
- Thread (tx, ty) handles i = ty*16 + ii, j = tx*4 + jj.
- Warp: 32 consecutive tids. Since tx varies fastest, a warp covers (tx=0..31, ty=0) and maybe (tx=0..31, ty=1) depending on how CUDA forms warps. Actually, with 2D threads, warps are formed in x-dimension first. So a warp has threads with the same ty but tx=0..31 if blockDim.x >= 32. Here blockDim.x=32. So warp covers 32 threads in one row of the block.
- Warp has ty fixed, tx=0..31.
- For jj=0: j = tx*4 + 0 = 0, 4, 8, ..., 124. Banks = j % 32 = 0, 4, 8, ..., 28. All different! No conflict!
What about A_smem? In the warp, ty is fixed, so all threads have the same ty. i = ty*16 + ii. For different ii in the inner loop, all threads access the same set of i values. For a specific ii, all threads read A_smem[k][ty*16 + ii]. They all read the same address! This is a broadcast, which is handled specially by shared memory and is fully parallel (no bank conflicts). Excellent!
So with Block=(32, 8), TileM=128, TileN=128, ItemsPerThreadM=16, ItemsPerThreadN=4:
- A_smem accesses: all threads in a warp read the same element (broadcast). No conflict.
- B_smem accesses: consecutive lanes read j = 0, 4, 8, ... Banks 0, 4, 8, ..., 28. All different if the bank calculation is j % 32. Wait, 0, 4, 8, 12, 16, 20, 24, 28, 0, 4, 8, ... for tx=0..31. tx=0: bank 0, tx=8: bank 0. Conflict!
Hmm wait: tx goes from 0 to 31. j = tx * 4. For tx=0: j=0, bank=0. tx=8: j=32, bank=0. tx=16: j=64, bank=0. tx=24: j=96, bank=0. 4 threads hit bank 0. 4-way conflict.
To fix: make ItemsPerThreadN not a multiple of 32's factors with 4? Actually, conflict occurs when 4 * tx differs by a multiple of 32. This happens when tx differs by 8. So 4-way conflict.
Can we pad B_smem? If B_smem[k][j] is padded, address = k*(TileN + pad) + j. Bank = (k*(TileN+pad) + j) % 32. For k=0, bank = j % 32. Same issue. For k>0, if (TileN+pad) is odd, then banks shift by k. For a fixed warp with varying j but same k, if k*(TileN+pad) ≡ 0 mod 32? We'd need this to shift banks.
Let's try TileN + pad = 129. Then for k=1, address offset = 129. 129 % 32 = 1. So banks for k=1 are (j + 1) % 32. This shifts by 1. But within a warp, the relative j spacing is still 4, so conflicts still happen (just shifted).
What if we do TileN=128 and ThreadsPerBlockX=32 but each warp handles a sub-range? Actually, warps always have 32 threads.
The 4-way conflict for B_smem might not be terrible. Shared memory bandwidth on H100 is very high. A 4-way conflict means 4 cycles instead of 1 for that access. But since the compute is memory-bound on output writes, shared memory bandwidth might not be the bottleneck.
Alternatively, we can use 4 warps within a block to load different B values and compute using a different approach. Or we can use `__shfl_sync` to share B values.
Actually, another approach: instead of each thread loading B_smem individually, we can load B values into registers via warp shuffles or using cooperative loading. But that's more complex.
Let me check if a 4-way shared memory bank conflict is actually a problem. H100 shared memory bandwidth is ~10-20 TB/s. A 4-way conflict on a read gives effective bandwidth of ~2.5-5 TB/s. But each warp reads 32 floats = 128 bytes. At 1 cycle per 128 bytes, 4 cycles = 4 clock cycles. For H100 at ~2 GHz, that's 2 ns. The compute loop has many FMAs, so the shared memory latency is hidden by the FMAs and the unroll.
Actually, I think the 4-way conflict is minor. Let's not over-optimize it unless profiling shows it's a problem. Let's also make sure we use padding or a different layout if needed.
Wait, another idea: what if we reorder B_smem as `__shared__ float B_smem[TileN][32];`? Then B_smem[j][k] for fixed k and varying j has address j*32 + k. Bank = (j*32 + k) % 32 = k % 32. Same bank for all j! Conflict.
What about `__shared__ float B_smem[TileN/4][32][4];` or using float4? Store B_smem as array of float4: `float4 B_smem[32][TileN/4]`. Then each float4 is 16 bytes = 4 consecutive banks? Actually, shared memory banks are 4 bytes each. A float4 span is 16 bytes = 4 banks. When a thread reads a float4, it uses 4 consecutive banks. If warp threads read float4s with offset 1 (i.e., j values 0,1,2,3,... per thread), then thread 0 reads banks 0,1,2,3, thread 1 reads banks 4,5,6,7, ..., thread 7 reads banks 28,29,30,31, thread 8 reads banks 0,1,2,3. Conflict! 4-way conflict for float4.
But our j pattern is j = tx * 4, so thread tx reads B_smem[k][tx*4], B_smem[k][tx*4+1], B_smem[k][tx*4+2], B_smem[k][tx*4+3] in an unrolled loop. For tx=0: banks 0,1,2,3. tx=1: banks 4,5,6,7. tx=7: banks 28,29,30,31. tx=8: banks 0,1,2,3. Again 4-way conflict.
Hmm. What if we process B in float4 and ensure tx accesses are to non-overlapping bank sets? If we have only 8 threads accessing 8 float4s, covering 32 banks, that's fine. But a warp has 32 threads.
What if we process B row in chunks where 8 threads collaborate on one set of 32 banks? For example, 8 threads read 8 float4s covering 32 banks. The other 24 threads in the warp do something else? We can't easily do that.
Actually, let me reconsider the preloaded approach. Shared memory bank conflicts might not be a significant issue compared to the overall memory throughput. And cuBLAS might also have similar issues or use different strategies.
Let me go back to the outer product approach with iteration over k. In that approach, A_smem and B_smem are 1D arrays: A_smem[TileM] and B_smem[TileN].
A warp accesses:
- A_smem[i] where i = ty*8 + ii. For fixed ii, all threads with same ty access same value? No, wait. In the outer product, thread (tx, ty) computes outputs for i = ty*ItemsPerThreadM + sub_i and j = tx*ItemsPerThreadN + sub_j. For a specific sub_i, thread reads A_smem[i]. In a warp (tx=0..15, ty=0,1), for sub_i=0: i values are 0 (for ty=0) and 8 (for ty=1). Wait, if ItemsPerThreadM=8, ty=0 handles i=0..7, ty=1 handles i=8..15. For a specific sub_i, say sub_i=0: ty=0 reads A_smem[0], ty=1 reads A_smem[8]. These are different elements. Banks 0 and 8. No conflict... unless other threads also hit these. In a warp of 32 threads: ty=0, tx=0..15 all read A_smem[0]. Ty=1, tx=0..15 all read A_smem[8].
Oh! For a fixed sub_i, all threads with the same ty read the same A_smem element! This is a broadcast. On modern GPUs, broadcasts are handled efficiently and don't cause bank conflicts (all threads read the same bank, but hardware can broadcast it). So A_smem is fine.
For B_smem[j] where j = tx*ItemsPerThreadN + sub_j. For fixed sub_j: tx=0 reads B_smem[0], tx=1 reads B_smem[8] (if ItemsPerThreadN=8). Banks 0 and 8. For tx=0..15: banks 0, 8, 16, 24, 0, 8, 16, 24, ... 4-way conflict!
So in the outer product approach, B_smem has a 4-way conflict with ItemsPerThreadN=8. But if we process 4 K values at once, the inner loop has 4 reads from B_smem per K-group iteration? No, in the grouped-k approach, for a K-group of 4, we still access B_smem[0][j], B_smem[1][j], B_smem[2][j], B_smem[3][j]. Each of these is a separate 4-way conflict.
But again, is this a big deal? Probably not for a memory-bound kernel.
To reduce B_smem conflicts, we can make ItemsPerThreadN odd or not a multiple of 32's factors. For example, ItemsPerThreadN=5. Then j = tx*5 + sub_j. Banks = (tx*5 + sub_j) % 32. For tx=0..31 and sub_j=0: banks 0, 5, 10, 15, 20, 25, 30, 3, 8, 13, 18, 23, 28, 1, 6, 11, 16, 21, 26, 31, 4, 9, 14, 19, 24, 29, 2, 7, 12, 17, 22, 27. All 32 banks, no conflicts!
But ItemsPerThreadN=5 means threads compute 5 outputs in the N dimension. For TileN=128, we'd need 128/5 = 25.6 threads. That's not integer. We could use TileN=128 with 32 threads and ItemsPerThreadN=4. Or TileN=160 with 32 threads and ItemsPerThreadN=5? But N=32768, and TileN needs to divide N evenly. 32768 / 160 = 204.8, not integer.
What about TileN=128, threads in N dimension = 32, ItemsPerThreadN=4? Then j=tx*4, banks 0,4,8,... 4-way conflict.
TileN=256, threads in N dim = 32, ItemsPerThreadN=8. j=tx*8, banks 0,8,16,24. 4-way conflict.
What if we use 64 threads in N dimension and ItemsPerThreadN=2? j=tx*2 + sub_j. For sub_j=0: banks 0,2,4,...,30,0,2,... 2-way conflict. For sub_j=1: banks 1,3,5,...,31,1,3,... 2-way conflict.
So no matter what, if ItemsPerThreadN is even and threads are powers of 2, we'll have conflicts. The minimum conflict is 2-way for ItemsPerThreadN=2, 4-way for ItemsPerThreadN=4, etc.
Actually, what if we add padding to B_smem? `float B_smem[TileN + 1]`. Then address = j + something. Bank = (j + k_offset) % 32. Doesn't help within a single k iteration.
What if we transpose the threads-to-j mapping? Instead of j = tx * ItemsPerThreadN + sub_j, we use j = sub_j * BlockDimX + tx. Then for sub_j=0, j=tx. Banks = tx % 32. All different, no conflict!
But then the memory layout in registers/accumulator is different. We can iterate sub_j in the outer loop and tx in the inner loop:
```cuda
for (int sub_j = 0; sub_j < ItemsPerThreadN; ++sub_j) {
int j = sub_j * BlockDimX + threadIdx.x;
float b = B_smem[j];
for (int sub_i = 0; sub_i < ItemsPerThreadM; ++sub_i) {
int i = sub_i * BlockDimY + threadIdx.y;
float a = A_smem[i];
accum[sub_i][sub_j] += a * b;
}
}
```
For a warp (ty fixed, tx=0..31):
- For each sub_j, all threads read B_smem[sub_j * 32 + tx]. Banks = (sub_j * 32 + tx) % 32 = tx % 32. All different! No conflict.
- For each sub_i, all threads read A_smem[sub_i * BlockDimY + ty]. All threads in the warp have the same ty, so they all read the same A_smem element. Broadcast! No conflict.
This is perfect! And the compute order changes slightly but is equivalent.
So the compute loop should be:
```cuda
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = jj * BlockDimX + threadIdx.x;
float b = B_smem[j];
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = ii * BlockDimY + threadIdx.y;
float a = A_smem[i];
accum[ii][jj] += a * b;
}
}
```
With Block=(32, 8), ItemsPerThreadM=16, ItemsPerThreadN=4:
- TileM = 8 * 16 = 128
- TileN = 32 * 4 = 128
- accum per thread: 16 * 4 = 64 floats.
For this to work, A_smem must have no conflicts for broadcast reads, which it doesn't.
B_smem[jj*32 + tx] has banks (jj*32 + tx) % 32 = tx % 32. Since tx ranges from 0 to 31 across the warp, each thread hits a different bank. Perfect!
Now, can we vectorize the inner loop with this reordering? The inner loop over ii does 16 iterations, each with the same b but different a. We can unroll this. The compiler will handle it.
Now let's write the complete kernel with all these optimizations.
Parameters: I'll try multiple configurations or settle on one that seems solid. Let's use:
- TileM = 128
- TileN = 128
- BlockDimX = 32
- BlockDimY = 8
- ItemsPerThreadM = 16
- ItemsPerThreadN = 4
- K = 32 (from input)
Wait, K isn't fixed at 32 for the kernel - it depends on the input shapes. But in the problem, get_inputs() creates A (32768, 32) and B (32, 32768), so K=32. However, the Model should be general. Let me check if K is always 32 or if it varies. The reference says "where one of the matrices is tall and skinny". The get_inputs() specifically uses N=32 for one dimension and M=32768 for the other. So in this test, K=32. But the Model should handle general K for correctness in other tests.
To make the kernel general, we need to handle arbitrary K. We can iterate k from 0 to K-1.
For vectorized A loads: A[i, k:k+4] is contiguous if we process 4 k's at once. But in the outer product, we load A[:, k] for one k at a time. To vectorize A loads, we need to change the iteration pattern to process 4 k values together.
With grouped K processing in chunks of 4:
```cuda
for (int k = 0; k < K; k += 4) {
// Load A_val[TileM][4] and B_val[4][TileN]
// A_smem[i][kk] for kk=0..3
// B_smem[kk][j] for kk=0..3
// In compute:
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
for (int kk = 0; kk < 4; ++kk) {
float b = B_smem[kk][jj * BlockDimX + threadIdx.x];
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
float a = A_smem[ii * BlockDimY + threadIdx.y][kk];
accum[ii][jj] += a * b;
}
}
}
}
```
For B_smem[kk][j], stored as B_smem[kk][j] = B[(k+kk)*N + block_j + j]. We want B_smem to avoid bank conflicts. B_smem is accessed as B_smem[kk][jj*32 + tx]. Address = kk * TileN + jj*32 + tx. Bank = (kk * TileN + jj*32 + tx) % 32. For TileN=128, kk*128 % 32 = 0. So bank = tx % 32. No conflict! But for kk=1, address = 128 + j, bank = j % 32. Same pattern. Good.
For A_smem[ii*8 + ty][kk]: address = (ii*8 + ty) * 4 + kk in a [TileM][4] array. Banks = ((ii*8 + ty)*4 + kk) % 32. For fixed ii and kk, varying ty: banks = (ty*4 + const) % 32. ty ranges 0..7. Banks 0,4,8,12,16,20,24,28 or shifted. These are all different! No conflict.
Great! So with this layout:
- `float A_smem[TileM][4]` - actually we need `float A_smem[TileM][4]` where the inner dimension is the kk group.
- `float B_smem[4][TileN]`
For A_smem, threads in a warp have the same ii and kk but different ty (since ty varies in y-dim and a warp spans tx=0..31 for fixed ty... wait, with Block=(32,8), threadIdx = ty*32 + tx. A warp has tid=0..31, which is ty=0, tx=0..31. Next warp is ty=1, tx=0..31? No, 256 threads = 8 warps. With tx varying fastest, warps are:
- Warp 0: ty=0, tx=0..31
- Warp 1: ty=1, tx=0..31
- ...
- Warp 7: ty=7, tx=0..31
So a warp has fixed ty! And tx varies 0..31.
In the compute loop, for a specific ii and kk, thread reads A_smem[ii*8 + ty][kk]. Since ty is FIXED in a warp, ALL threads in a warp read the exact same A_smem element! Broadcast. No conflict.
For B_smem[kk][jj*32 + tx], tx varies 0..31. Banks = (kk*128 + jj*32 + tx) % 32.
- For kk*128 % 32: 128 = 4*32, so kk*128 % 32 = 0.
- jj*32 % 32 = 0.
- Bank = tx % 32 = tx.
All different. No conflict!
This is actually a perfect configuration!
Now, for loading A_smem from global memory:
A_smem[idx][kk] = A[(block_i + idx) * K + k + kk]
For idx in 0..127 and kk in 0..3.
Total elements: 512. With 256 threads, each thread loads 2 elements.
Thread tid can load:
idx = tid, kk = 0 if tid < 128? Or we can iterate.
for (int idx = tid; idx < TileM * 4; idx += BlockSize) {
int i = idx / 4;
int kk = idx % 4;
A_smem[i][kk] = A[(block_i + i) * K + k + kk];
}
But wait, A[i][k+kk] for i varying by 1 and kk varying: addresses are (base + i) * K + k + kk. For consecutive idx values, we might have i changing and kk wrapping. The memory accesses might not be fully coalesced.
Actually, if we structure loads differently:
for (int i = tid; i < TileM; i += BlockSize) {
// Load 4 floats from A[i, k:k+4] as float4
*reinterpret_cast<float4*>(&A_smem[i][0]) =
*reinterpret_cast<const float4*>(&A[(block_i + i) * K + k]);
}
Since each row's k:k+4 floats are contiguous, float4 loads are perfectly coalesced!
Consecutive threads (tid=0,1,2,...) load consecutive i's, so consecutive float4s from memory. This is excellent!
For B_smem:
B_smem[kk][j] = B[(k+kk)*N + block_j + j]
For each kk in 0..3, the segment B[(k+kk)*N + block_j : (k+kk)*N + block_j + TileN] is TileN contiguous floats.
We can load each of the 4 segments with float4s.
for (int kk = 0; kk < 4; ++kk) {
for (int j = tid; j < TileN; j += BlockSize) {
B_smem[kk][j] = B[(k + kk) * N + block_j + j];
}
}
With 256 threads and TileN=128, each thread loads 0.5 elements on average for each kk. Or we can use vectorized loads:
for (int kk = 0; kk < 4; ++kk) {
float4* B_smem4 = reinterpret_cast<float4*>(&B_smem[kk][0]);
const float4* B_global4 = reinterpret_cast<const float4*>(&B[(k + kk) * N + block_j]);
for (int j = tid; j < TileN / 4; j += BlockSize) {
B_smem4[j] = B_global4[j];
}
}
TileN=128, TileN/4=32. With 256 threads, each thread loads... wait, only 32 float4s total for each kk. With 256 threads, most threads won't participate. We need to be careful.
Actually, we can load all 4 kk segments together:
int total_float4s = 4 * (TileN / 4) = 128 float4s = 512 floats.
With 256 threads, each thread loads 2 float4s = 8 floats. Very good coalescing.
for (int idx = tid; idx < 4 * (TileN / 4); idx += BlockSize) {
int kk = idx / (TileN / 4);
int j4 = idx % (TileN / 4);
reinterpret_cast<float4*>(B_smem)[idx] =
reinterpret_cast<const float4*>(B)[(k + kk) * (N / 4) + block_j / 4 + j4];
}
Wait, B is (K, N). B[(k+kk)*N + block_j + j4*4 ... j4*4+3]. In terms of float4, if N is divisible by 4 and block_j is divisible by 4, then B[(k+kk)*N + block_j + j*4] is a valid float4 load.
For N=32768, yes, divisible by 4. For block_j = blockIdx.y * TileN, and TileN=128, block_j is divisible by 128, hence by 4. Great.
So vectorized loads of both A and B are possible!
For A: each thread loads float4 from `A + (block_i + i) * K + k`.
For B: each thread loads float4 from `B + (k + kk) * N + block_j + j*4`.
With 256 threads, we can load the entire A and B segments very efficiently.
For the compute loop:
```cuda
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
for (int kk = 0; kk < 4; ++kk) {
float b = B_smem[kk][jj * BlockDimX + threadIdx.x];
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
float a = A_smem[ii * BlockDimY + threadIdx.y][kk];
accum[ii][jj] += a * b;
}
}
}
```
This has 4 * ItemsPerThreadN * ItemsPerThreadM FMAs per outer k-iteration.
For ItemsPerThreadM=16, ItemsPerThreadN=4: 4 * 4 * 16 = 256 FMAs per k-group.
With K=32, there are 8 k-groups. Total 2048 FMAs per thread.
Wait, earlier I confirmed 2048 FMAs = 64 outputs × 32 k values = 2048. But with kk loop of 4, it's ItemsPerThreadM × ItemsPerThreadN × 4 per group, and 8 groups: 16 × 4 × 4 × 8 = 2048. Correct.
The inner loop over ii (16 iterations) accesses A_smem with fixed kk and varying ii. But since all threads in a warp have the same ty, and the ii loop is the same for all threads in the warp, A_smem accesses are broadcast. No bank conflicts.
B_smem is accessed as B_smem[kk][jj*32 + tx]. For fixed kk and jj, tx varies. Banks = tx % 32. No conflict.
FMAs per shared memory read: total reads per k-group = ItemsPerThreadN * 4 B-reads + ItemsPerThreadM * 4 * ItemsPerThreadN A-reads? Wait, B_smem is read ItemsPerThreadN * 4 times (one per jj and kk). A_smem is read ItemsPerThreadM * 4 * ItemsPerThreadN times (one per ii, jj, kk). Total reads = 16 + 16*4*4 = 16 + 256 = 272 reads. FMAs = 256. Not great.
Wait, there's reuse that the compiler should optimize. For fixed jj and kk, b is loaded once and reused across the ii loop. So B reads = 4 * 4 = 16. For fixed ii, a is loaded once per kk,jj combination? Actually in the current loop order, for each kk, b is loaded and then the ii loop runs with a different a each time. A reads = 4 (kk) * 4 (jj) * 16 (ii) = 256. Total reads = 272 for 256 FMAs. Almost 1:1.
Can we improve the A reuse? What if we reorder:
```cuda
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
for (int kk = 0; kk < 4; ++kk) {
float a = A_smem[ii * BlockDimY + threadIdx.y][kk];
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
float b = B_smem[kk][jj * BlockDimX + threadIdx.x];
accum[ii][jj] += a * b;
}
}
}
```
A reads = 16 * 4 = 64. B reads = 4 * 4 = 16. Total reads = 80 for 256 FMAs. Much better! 3.2:1 ratio.
But does this have bank conflicts?
- For A_smem[ii*8 + ty][kk]: fixed ii, fixed ty, varying kk across iterations. In a warp, all threads read the same A_smem[ii*8 + ty] element but different kk. Wait, A_smem[ii*8 + ty] is a float4 (or 4 floats). If stored as `float A_smem[128][4]`, then A_smem[i][0], A_smem[i][1], A_smem[i][2], A_smem[i][3] are 4 consecutive floats in memory. Banks = (i*4 + kk) % 32. For a warp (fixed i = ii*8 + ty), all threads read the same i. So they read addresses i*4 + kk for kk=0,1,2,3 across iterations. These are consecutive banks. But within a single kk iteration, all threads read the same address? No, they read A_smem[i][kk] where i is the same for all threads in the warp (since ty is fixed). So it's a broadcast. No conflict.
But wait, with the reordered loop, for fixed ii and kk, all threads in the warp read A_smem[ii*8 + ty][kk]. This is the same address for all threads. Broadcast. Good.
For B_smem[kk][jj*32 + tx]: fixed kk, varying jj and tx. For a specific jj, tx varies 0..31. Banks = (kk*128 + jj*32 + tx) % 32 = tx % 32. All different. No conflict.
This is perfect! The reordered loop is better.
However, I need to check: in the reordered loop, for fixed ii and kk, we iterate jj from 0..3. Each iteration loads a new B_smem value and accumulates into a different output. The compiler can unroll this and interleave loads and FMAs.
But there's a concern: with 16 ii values and 4 kk values, the inner jj loop runs 4 times. So 16*4*4 = 256 iterations, each doing 1 FMA. Plus 64 A reads and 16 B reads. Actually, 64 A reads and 64 B reads (wait: for each ii and kk, jj=0..3, so B is read 16*4*4 = 64 times).
Wait, the jj loop is inside:
for ii (16):
for kk (4):
a = A_smem[...][kk]
for jj (4):
b = B_smem[kk][jj*32 + tx]
accum[ii][jj] += a * b;
B reads: 16 * 4 * 4 = 64.
A reads: 16 * 4 = 64.
Total reads: 128 for 256 FMAs. 2:1 ratio. Good.
And the accumulators are 16*4 = 64 floats. Each one is updated 4 times (for kk=0..3) per ii? Actually, for each (ii,jj), accum is updated for kk=0,1,2,3. So 4 FMAs per (ii,jj) pair, one per kk. The reordered loop has better A reuse because a is loaded once and used for 4 jj values.
This should be very efficient.
Now for writing C. Each thread has 64 outputs. We want to write them efficiently.
C[i, j] where i = block_i + ii*8 + ty, j = block_j + jj*32 + tx.
For fixed ii and jj, threads in a warp have consecutive j values (jj*32 + tx for tx=0..31). These are contiguous floats! We can write them as float4 or float2 or just coalesced float stores.
Actually for jj=0, j = tx. For jj=1, j = 32 + tx. These are separated by 32 floats = 128 bytes. So each jj group writes a contiguous 32-float segment.
If we write jj=0 as a float4×8 or just 32 float stores, they're perfectly coalesced.
Then jj=1 writes the next 32 floats. Also coalesced.
But can we vectorize to float4? Each thread has 4 outputs across jj=0..3. For a fixed ii, thread tx writes:
C[i][tx], C[i][32+tx], C[i][64+tx], C[i][96+tx].
These are NOT contiguous! They're spaced by 32 floats.
So within a thread, vectorization is not straightforward. But across the warp, the writes are coalesced. Each warp writes 32 contiguous floats for each jj. So just regular float stores are fine and fully coalesced.
However, for better instruction efficiency, we could have threads write float4 by rearranging the output. But that's complex and probably not worth it.
Actually, what if we change ItemsPerThreadN to 4 but ensure that for each thread, its 4 j values are contiguous? That would mean j = threadIdx.x * ItemsPerThreadN + jj. But then we'd have bank conflicts on B_smem.
Trade-offs... I think avoiding bank conflicts in the compute loop is more important than vectorizing stores. The compute loop runs many more times than the store loop (8 iterations of k-group vs 1 store). And stores are already coalesced.
Let me write up the kernel. I'll define a template or just hardcode some values. Since the kernel needs to handle general M, K, N (as long as K is small), I'll pass them as arguments.
Actually, looking at the problem again: the model is just `torch.matmul(A, B)` where A is 2D and B is 2D. So M, K, N are determined at runtime from the input shapes. We need a kernel that handles arbitrary M, K, N.
For tiling, we need M divisible by TileM and N divisible by TileN. We can add bounds checking or require alignment. For the test case, M=32768 and N=32768, both divisible by 128. For generality, I'll add bounds checking.
Actually, for simplicity, I'll write the kernel with if-guards for the boundaries. Given the test M and N are powers of 2, they should align.
Let me draft the CUDA code:
```cuda
#include <torch/extension.h>
#include <cuda_runtime.h>
template <int TileM, int TileN, int BlockDimX, int BlockDimY, int VecK>
__global__ void gemm_kernel(const float* __restrict__ A, const float* __restrict__ B, float* __restrict__ C,
int M, int K, int N) {
static_assert(TileM % BlockDimY == 0, "TileM must be divisible by BlockDimY");
static_assert(TileN % BlockDimX == 0, "TileN must be divisible by BlockDimX");
static_assert(VecK == 4, "Only VecK=4 supported for float4 loads");
constexpr int ItemsPerThreadM = TileM / BlockDimY;
constexpr int ItemsPerThreadN = TileN / BlockDimX;
int block_i = blockIdx.x * TileM;
int block_j = blockIdx.y * TileN;
int tid = threadIdx.y * BlockDimX + threadIdx.x;
__shared__ alignas(16) float A_smem[TileM][VecK];
__shared__ alignas(16) float B_smem[VecK][TileN];
float accum[ItemsPerThreadM][ItemsPerThreadN] = {0.0f};
for (int k_base = 0; k_base < K; k_base += VecK) {
// Load A_smem
// Each row has VecK floats = VecK * 4 bytes. Use float4.
int num_a_vec4 = TileM; // TileM rows, each loads 1 float4
for (int idx = tid; idx < num_a_vec4; idx += BlockDimX * BlockDimY) {
int i = idx;
if (block_i + i < M && k_base < K) {
const float4* A_ptr = reinterpret_cast<const float4*>(A + (block_i + i) * K + k_base);
float4 val = *A_ptr;
A_smem[i][0] = val.x;
A_smem[i][1] = val.y;
A_smem[i][2] = val.z;
A_smem[i][3] = val.w;
}
}
// Load B_smem
// 4 rows of B, each with TileN floats. Load as float4.
int num_b_vec4 = VecK * (TileN / 4);
for (int idx = tid; idx < num_b_vec4; idx += BlockDimX * BlockDimY) {
int kk = idx / (TileN / 4);
int j4 = idx % (TileN / 4);
if (k_base + kk < K && block_j + j4 * 4 < N) {
const float4* B_ptr = reinterpret_cast<const float4*>(B + (k_base + kk) * N + block_j + j4 * 4);
float4 val = *B_ptr;
B_smem[kk][j4 * 4 + 0] = val.x;
B_smem[kk][j4 * 4 + 1] = val.y;
B_smem[kk][j4 * 4 + 2] = val.z;
B_smem[kk][j4 * 4 + 3] = val.w;
}
}
__syncthreads();
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
#pragma unroll
for (int kk = 0; kk < VecK; ++kk) {
int i = ii * BlockDimY + threadIdx.y;
float a = A_smem[i][kk];
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = jj * BlockDimX + threadIdx.x;
float b = B_smem[kk][j];
accum[ii][jj] += a * b;
}
}
}
__syncthreads();
}
// Write output
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = block_i + ii * BlockDimY + threadIdx.y;
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = block_j + jj * BlockDimX + threadIdx.x;
if (i < M && j < N) {
C[i * N + j] = accum[ii][jj];
}
}
}
}
```
Wait, I need to make sure float4 alignment is correct. A + (block_i + i) * K + k_base. For this to be aligned to 16 bytes, (block_i + i) * K + k_base must be divisible by 4 (since float4 requires 16-byte alignment). K is a runtime value. In the test, K=32, so (block_i + i) * 32 + k_base is divisible by 4 if k_base is divisible by 4. For k_base starting at 0 and incrementing by 4, k_base is divisible by 4. And 32 is divisible by 4, so (block_i + i) * 32 is always divisible by 4. Great!
For B: B + (k_base + kk) * N + block_j + j4 * 4. N=32768, which is divisible by 4. k_base is divisible by 4. block_j = blockIdx.y * TileN = blockIdx.y * 128, divisible by 4. j4 * 4 is divisible by 4. So total offset is divisible by 4. Alignment satisfied!
For A_smem[i][VecK], the inner dimension VecK=4 means row i starts at i*4. If we want to load float4 directly into A_smem[i], we need A_smem[i] to be aligned to 16 bytes. Since A_smem is `alignas(16) float[TileM][4]`, the base is aligned. And each row is 16 bytes, so row i is at i*16, which is 16-byte aligned. Great!
But wait, `float A_smem[TileM][4]` - when I do `reinterpret_cast<float4*>(&A_smem[i][0])`, is that valid? In C++, yes, because A_smem[i] is an array of 4 floats, which is contiguous. And with the alignment, it should be fine.
Actually, for the float4 store into A_smem, I could do:
```cuda
float4* A_smem_ptr = reinterpret_cast<float4*>(&A_smem[i][0]);
*A_smem_ptr = val;
```
instead of scalar assignments. This is more efficient.
Similarly for B_smem:
```cuda
float4* B_smem_ptr = reinterpret_cast<float4*>(&B_smem[kk][j4 * 4]);
*B_smem_ptr = val;
```
Now, one issue: the kernel assumes K is divisible by VecK. For general K, we need a tail loop or bounds checking. Since K might not be a multiple of 4, we should check `k_base + kk < K` in the compute loop. Let me add a mask or handle it.
For simplicity, I'll handle arbitrary K by processing in chunks of 4 with masking. In the compute loop, if k_base + kk >= K, we should skip or treat as 0.
Actually, if K is not a multiple of 4, we can still load the full 4 floats from A, but some might be out of bounds (reading past K). This is bad. We need to guard the loads.
Let me add a mask in the compute loop:
```cuda
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
#pragma unroll
for (int kk = 0; kk < VecK; ++kk) {
int i = ii * BlockDimY + threadIdx.y;
float a = A_smem[i][kk];
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = jj * BlockDimX + threadIdx.x;
float b = B_smem[kk][j];
if (k_base + kk < K) {
accum[ii][jj] += a * b;
}
}
}
}
```
Wait, but the loads might have read out of bounds. For A, if K is not a multiple of 4 and k_base + 3 >= K, then A + (block_i + i) * K + k_base + 3 is out of bounds. We should only load valid elements.
Better approach: for A loads, check bounds. The float4 load is tricky because we can't easily do a masked float4 load. We can check if k_base + 3 < K, and if so, do float4. Otherwise, do scalar loads for valid elements.
For simplicity and since K is likely to be a multiple of 4 in the test (K=32), I'll handle the general case but with a scalar fallback for the final chunk.
Actually, let me look at the test inputs again. `A = torch.rand(M, N)` where M=32768, N=32. And `B = torch.rand(N, M)` where N=32, M=32768. So K for A is N=32, and the inner dimension for B is also N=32. So K=32, which is divisible by 4. Great.
For robustness, I'll add the general case handling just in case.
Actually, I realize the float4 load of A might be problematic if K isn't divisible by 4. I'll add a scalar fallback code path.
Or better yet, I'll just use a simpler scalar loading approach since K is only 32 and the amount of data is small. Vectorization is nice but might add complexity. However, scalar loads for A might be fine since each thread only loads 1 float (A_smem has 512 floats, 256 threads, 2 loads per thread). Actually with the scalar approach:
```cuda
for (int idx = tid; idx < TileM * VecK; idx += BlockDimX * BlockDimY) {
int i = idx / VecK;
int kk = idx % VecK;
if (block_i + i < M && k_base + kk < K) {
A_smem[i][kk] = A[(block_i + i) * K + k_base + kk];
}
}
```
For TileM=128, VecK=4: 512 elements. 256 threads each load 2 elements. Consecutive threads load consecutive elements? idx=0: i=0, kk=0. idx=1: i=0, kk=1. These are consecutive in memory for A[0*K + k_base + 0] and A[0*K + k_base + 1]. Yes! For K=32, these are consecutive floats. So consecutive threads load consecutive floats. This is coalesced!
For B_smem:
```cuda
for (int idx = tid; idx < VecK * TileN; idx += BlockDimX * BlockDimY) {
int kk = idx / TileN;
int j = idx % TileN;
if (k_base + kk < K && block_j + j < N) {
B_smem[kk][j] = B[(k_base + kk) * N + block_j + j];
}
}
```
For TileN=128, VecK=4: 512 elements. 256 threads each load 2 elements.
idx=0: kk=0, j=0. Address: B[k_base * N + block_j + 0].
idx=1: kk=0, j=1. Address: B[k_base * N + block_j + 1].
Consecutive threads load consecutive floats. Coalesced!
So scalar loads are actually perfectly coalesced for this configuration. And they're simpler and safer for bounds checking. No need for float4.
The only advantage of float4 is fewer instructions. But with 2 loads per thread, it's negligible. Let's use scalar loads for simplicity.
But actually, wait. For A_smem loads with idx = tid and 256 threads:
idx=0: A[0*32 + k_base + 0]
idx=1: A[0*32 + k_base + 1]
idx=2: A[0*32 + k_base + 2]
idx=3: A[0*32 + k_base + 3]
idx=4: A[1*32 + k_base + 0]
...
These are not all contiguous! After the first 4 elements, the next is 32 floats away. For a warp of 32 threads:
Thread 0: A[0*32 + k]
Thread 1: A[0*32 + k+1]
...
Thread 3: A[0*32 + k+3]
Thread 4: A[1*32 + k+0]
...
Thread 7: A[1*32 + k+3]
Thread 8: A[2*32 + k+0]
...
Addresses: base, base+4, base+8, base+12, base+128, base+132, base+136, base+140, base+256...
These are strided by 4 bytes for the first 4 threads, then jump by 112 bytes, then stride by 4, etc.
In a warp, threads 0..3 access addr 0,4,8,12 (all in cache line 0).
Threads 4..7 access addr 128,132,136,140 (all in cache line 1).
Threads 8..11 access addr 256,260,264,268 (cache line 2).
...
Threads 28..31 access addr 896,900,904,908 (cache line 7).
So a warp accesses 8 cache lines. The addresses are perfectly sequential! The memory controller should merge these into one big burst read. This is actually very efficient.
For B_smem with scalar loads:
Thread 0: B[k_base * N + block_j + 0]
Thread 1: B[k_base * N + block_j + 1]
...
Thread 255: B[(k_base + 7) * N + block_j + 127] (if VecK=4, TileN=128, then idx goes 0..511).
Wait, with 256 threads and 512 elements, threads 0..127 handle kk=0,j=0..127 and threads 128..255 handle kk=1,j=0..127.
Thread 0: B[k_base*N + block_j + 0]
Thread 1: B[k_base*N + block_j + 1]
...
Thread 127: B[k_base*N + block_j + 127]
Thread 128: B[(k_base+1)*N + block_j + 0]
These are two contiguous regions of 128 floats each. Thread 128 starts a new region. If thread 128 is in the same warp as thread 96 (which is warp 3: threads 96..127), then warp 3 has threads 96..127 and also... wait, warps are formed by consecutive thread IDs. So warp 0: threads 0..31, warp 1: 32..63, warp 2: 64..95, warp 3: 96..127, warp 4: 128..159, warp 5: 160..191, warp 6: 192..223, warp 7: 224..255.
Warp 3 (threads 96..127) all access B[k_base*N + block_j + 96..127]. Contiguous! One cache line (actually 128 bytes = 32 floats, so 4 cache lines for 128 floats). Wait, 32 floats = 128 bytes = 1 cache line. 96 floats = 384 bytes = 3 cache lines. Threads 96..127 access 32 floats = 1 cache line. Actually 32 threads × 1 float each = 128 bytes. Yes, one cache line if aligned.
Warp 4 (threads 128..159) access B[(k_base+1)*N + block_j + 0..31]. This is far away (offset by N floats from warp 3's address). For N=32768, this is 131072 bytes away. So warp 4 accesses a different region.
Each warp's accesses are fully coalesced (one cache line per warp). This is very good.
So scalar loads are fine and simpler.
Now let me think about whether there's an even better tile configuration.
Actually, one more optimization: C should use `__restrict__` and we should compile with `-O3` and `--use_fast_math` if appropriate. But for fp32 accumulation, we don't want fast math to change precision. Actually, `--use_fast_math` wouldn't affect simple FMAs. It's fine.
We should also set `__launch_bounds__` to help the compiler with register allocation.
For our kernel with 256 threads, if we set `__launch_bounds__(256)`, the compiler knows there are 256 threads and can optimize register usage accordingly. Actually, we might want to hint maximum registers to improve occupancy.
Another optimization: since K is small, we might not need the outer k loop at all if K is known at compile time. But K is dynamic. If we want, we can template on K or handle specific K values with if-else. But that's micro-optimization.
Let me write the full Python module now. I'll create a kernel with TileM=128, TileN=128, Block=(32,8), VecK=4. But wait, ItemsPerThreadN = TileN / BlockDimX = 128 / 32 = 4. ItemsPerThreadM = TileM / BlockDimY = 128 / 8 = 16. Accum = 64 registers.
But I want to try something: what if TileM=64, TileN=256, Block=(32,8)?
- ItemsPerThreadM = 64 / 8 = 8
- ItemsPerThreadN = 256 / 32 = 8
- Accum = 64 registers. Same.
- A_smem: 64 × 4 = 256 floats = 1KB (per k-group of 4, but actually we load all K? No, per k-base. Wait, with VecK=4, A_smem is 64×4=256 floats=1KB).
- B_smem: 4 × 256 = 1024 floats = 4KB.
- Total shared: 5KB. Very small.
- Blocks: (32768/64) × (32768/256) = 512 × 128 = 65536.
- Each block writes 64×256 = 16384 floats = 64KB.
Wait, that's interesting. Same amount of shared memory per block but different tile shape. More blocks might help with load balancing.
What about TileM=64, TileN=512, Block=(32,8)?
- ItemsPerThreadM = 8, ItemsPerThreadN = 16. Accum = 128 registers. Might be too many.
- Let's not risk it.
What about TileM=256, TileN=128, Block=(32,8)?
- ItemsPerThreadM = 32, ItemsPerThreadN = 4. Accum = 128 registers. Might be borderline.
- A_smem: 256 × 4 = 1024 floats = 4KB.
- B_smem: 4 × 128 = 512 floats = 2KB.
- Total: 6KB.
- Blocks: 128 × 256 = 32768.
I'll test TileM=128, TileN=128 first. If it's not fast enough, I'll try others.
Actually, I realize that one potential issue is that with K=32 and VecK=4, we have 8 iterations. Each iteration does `__syncthreads()` twice (after load and after compute). That's 16 barriers per block. A barrier has overhead. Can we remove the second syncthreads? Actually, we only need syncthreads before overwriting shared memory in the next iteration. Since all threads read from shared memory in the compute phase, and the next iteration overwrites it, we need a syncthreads after compute to ensure all reads are done. But wait, threads in the same block don't write to shared memory in the compute phase. So do we need syncthreads after compute?
Yes! Because some threads might finish the compute phase and move to the load phase of the next k_base iteration, while other threads are still reading from shared memory in the compute phase. We need the syncthreads to prevent this race condition.
Wait, actually: in the compute phase, all threads only read from shared memory (A_smem and B_smem). In the next iteration's load phase, threads write to shared memory. If a fast warp finishes the compute phase and starts writing to A_smem, a slow warp still reading from A_smem would get corrupted data. So syncthreads IS needed between compute and next load.
So we have load -> syncthreads -> compute -> syncthreads, repeated 8 times. 16 syncthreads total. This is a lot of overhead.
Is there a way to reduce this? One approach: double-buffer the shared memory. Use two buffers and swap them. Then:
- Load into buffer 0
- syncthreads
- Compute from buffer 0 while next group loads into buffer 1? No, we can't compute and load simultaneously within the same kernel without more complex pipelining.
Actually, since we have barrier overhead, maybe a different algorithm that doesn't use shared memory at all is better? Or one that loads directly from global memory with L2 caching?
For the outer product without shared memory, each thread would need to redundantly load A and B values. But with L2 caching, the first warp loads A[:, k] into L2, and subsequent warps hit L2. However, without shared memory, each thread computes independently and doesn't need syncthreads.
But for a thread block of 256 threads, we can assign each thread a subset of outputs and have each thread load the A and B values it needs directly from global memory.
- Thread handles 8×8 outputs.
- For each k in 0..31:
- Thread loads A[block_i + i_local, k] for i_local in its 8 rows. These are 8 floats, strided by 128B each.
- Thread loads B[k, block_j + j_local] for j_local in its 8 cols. These are 8 contiguous floats.
- Compute accumulations.
No syncthreads needed! But each thread loads 8+8=16 floats per k from global memory. With 256 threads, total global loads per k = 4096 floats. But each float is loaded redundantly (only cached in L2, not shared).
For A: 256 threads × 8 floats = 2048 A-floats per k. But there are only 128 unique A-floats per k in the block. So 16x redundancy!
For B: 256 threads × 8 floats = 2048 B-floats per k. Only 128 unique. 16x redundancy!
With L2 cache, the redundancy means L2 bandwidth is 16x higher. L2 bandwidth is ~10 TB/s, so 16x overhead on 512 bytes per k = 8KB per k. Not terrible.
But global memory reads from L2 are still 16x. 8KB × 32 = 256KB per block.
For 65536 blocks, total L2 traffic from A and B = 16GB. L2 bandwidth is ~10 TB/s, so time = 16GB / 10TB/s = 1.6ms. Plus 4GB DRAM write = additional 1.2ms. Total ~2.8ms.
With shared memory, global memory reads are optimal (only ~8MB total from DRAM, cached perfectly). But we have 16 syncthreads per block. If syncthreads takes ~100 cycles each (very rough estimate), that's 1600 cycles per block. At 2GHz, that's 0.8 microseconds per block. For 65536 blocks running on 132 SMs, each SM runs ~500 blocks. 500 × 0.8us = 0.4ms overhead. Actually, this is a rough estimate.
The shared memory approach should be faster because it minimizes memory traffic. Let's go with it.
Can we reduce syncthreads? Wait, actually I miscounted. For VecK=4 and K=32, we have 8 iterations. Each iteration:
1. Load A and B into shared memory
2. __syncthreads()
3. Compute
4. __syncthreads()
That's 8 iterations × 2 syncthreads = 16 syncthreads. But actually, for the first iteration, we don't need a syncthreads before loading (no one is reading yet). And for the last iteration, we don't need a syncthreads after computing (we're writing to C, not shared memory). But the loop structure makes this tricky.
Actually, with proper loop unrolling or peeling, we can eliminate 2 syncthreads:
- Iteration 0: load, syncthreads, compute
- Iterations 1..6: syncthreads, load, syncthreads, compute
- Iteration 7: syncthreads, load, syncthreads, compute
Wait, that's still 1 + 6*2 + 2 = 15. Or:
- Load A, B for k=0
- Loop k=0..K-1 step VecK:
- __syncthreads() // ensure previous compute done (skip for k=0)
- Compute from shared memory
- __syncthreads() // ensure compute done before next load
- Load A, B for next k (skip after last)
This requires pipelining the loads. For k=0, load before loop. Then in the loop, compute and load next. This is more complex.
Actually, for simplicity and since 16 syncthreads for a block that's doing a lot of compute is probably fine, I'll just keep both syncthreads.
Another idea: can we process the entire K dimension in one go without shared memory iteration? A_smem[128][32] = 4096 floats = 16KB. B_smem[4][128] per iteration, or B_smem[32][128] = 4096 floats = 16KB. Total shared = 32KB. This fits!
If we load the entire A and B tiles into shared memory at once:
- Load A_smem[TileM][32] = 16KB
- Load B_smem[32][TileN] = 16KB
- Then compute with no intermediate syncthreads!
But B_smem[32][128] has 32 rows, and in the compute loop we iterate k=0..31 accessing B_smem[k][j]. For consecutive k and fixed j, B_smem[k][j] addresses are k*128 + j. Banks = (k*128 + j) % 32 = j % 32 (since 128 % 32 = 0). For fixed j and varying k across iterations, each iteration hits the same bank! But within a single k iteration, different threads access different j, which might be spaced.
Wait, for our compute pattern with B_smem[k][jj*32 + tx], banks = (k*128 + jj*32 + tx) % 32 = tx % 32. No conflict within a k iteration. Across k iterations, the bank is the same, but that's fine since they're sequential.
So preloading all of B_smem and A_smem into shared memory at once, then computing with one syncthreads after load, sounds very appealing.
Let's recalculate: A_smem[128][32] = 4096 floats. B_smem[32][128] = 4096 floats. Total = 32KB. Under the 48KB limit for most configurations. On H100, we can use up to 227KB, so 32KB is fine.
With preloading:
- Load all A and B data (32KB total) using 256 threads.
- Scalar loads as before:
- A: 128 rows × 32 cols = 4096 elements. 256 threads, 16 loads each. Wait, 4096/256 = 16 per thread. That's quite a few loads.
- B: 32 rows × 128 cols = 4096 elements. 16 loads each.
- Total loads per thread: 32 scalar loads. At 2 loads per clock or so, this takes ~16 cycles. Very fast.
- __syncthreads()
- Compute: 32 iterations over k. Each iteration does 16 A-reads + 4 B-reads (for ItemsPerThreadM=16, ItemsPerThreadN=4) = 20 shared mem reads, 64 FMAs.
- Total compute: 32 * 64 = 2048 FMAs.
- Total shared reads: 32 * 20 = 640 reads. 3.2 FMAs per shared read. Good.
- Then write C.
This has only 1 syncthreads! Much better than 16. And all the data is in shared memory, no global memory traffic during compute.
Wait, actually, if we preload everything, we need to make sure shared memory is large enough. 32KB is fine.
Can we make the tiles bigger? TileM=128, TileN=256:
- A_smem: 128 × 32 = 4096 floats = 16KB
- B_smem: 32 × 256 = 8192 floats = 32KB
- Total: 48KB. Exactly the 48KB limit (which is often the default max for some block configs). On H100 we can go higher, but 48KB should be safe.
- Threads: 256.
- ItemsPerThreadM = 128 / BlockDimY. If BlockDimY=8, ItemsPerThreadM=16.
- ItemsPerThreadN = 256 / BlockDimX. If BlockDimX=32, ItemsPerThreadN=8.
- Accumulators = 16 × 8 = 128. More registers, but still probably OK (128 accum + some temps + loop ctrs = ~160? Maybe too many for 255 limit, but 256 threads would allow it?).
Actually, let's compute registers more carefully. With 128 float accumulators in registers: 128 registers. Plus loop variables (ii, jj, k), plus temporaries for a and b. The compiler might be able to unroll everything and keep it under 200. 256 threads × 200 = 51200. Under 65536. It might be OK.
But to be safe, let me try TileM=128, TileN=128 first with preloading.
Kernel with preloading:
```cuda
__global__ void gemm_kernel(const float* __restrict__ A, const float* __restrict__ B, float* __restrict__ C,
int M, int K, int N) {
constexpr int TileM = 128;
constexpr int TileN = 128;
constexpr int BlockDimX = 32;
constexpr int BlockDimY = 8;
constexpr int ItemsPerThreadM = TileM / BlockDimY; // 16
constexpr int ItemsPerThreadN = TileN / BlockDimX; // 4
int block_i = blockIdx.x * TileM;
int block_j = blockIdx.y * TileN;
int tid = threadIdx.y * BlockDimX + threadIdx.x;
__shared__ float A_smem[TileM][32]; // K dimension assumed <= 32? No, K is dynamic.
__shared__ float B_smem[32][TileN];
// But K could be larger than 32! We can't statically allocate [TileM][K].
```
Oh right! I forgot K is dynamic. We can only statically allocate with a compile-time constant maximum. If K is always small (like 32 or 64), we can set a max. But for a general kernel, this is problematic.
Actually, we can allocate 1D shared memory and index manually:
```cuda
__shared__ float smem[TileM * MaxK + MaxK * TileN];
```
Then `A_smem = smem` and `B_smem = smem + TileM * MaxK`. But we need MaxK at compile time.
If K is typically small (≤ 64 or 128), we can set MaxK=128. But then shared memory = 128*128 + 128*128 = 32768 floats = 128KB. That's getting large but still under 227KB.
But if K can be larger, say 1024, then we'd need MaxK=1024, shared mem = 1024*128 + 1024*128 = 262144 floats = 1MB. Way too much.
So dynamic shared memory is better. We can compute the shared memory size at kernel launch:
```cuda
extern __shared__ float smem[];
float* A_smem = smem;
float* B_smem = smem + TileM * K;
__shared__ ... // no, extern __shared__ must be the only shared array
```
Actually, with dynamic shared memory, we'd specify the size in the kernel launch. But then we need to know K at launch time.
However, dynamic shared memory with `extern __shared__` is a bit tricky with templates and load_inline. But we can pass the size to the kernel launch as the third argument to `<<<grid, block, shared_mem_size>>>`.
If K can be arbitrary, and the typical test case has K=32, we can still use dynamic shared memory. For K=32, shared mem = 128*32 + 32*128 = 8192 floats = 32KB.
Wait, but we need to be careful. For K larger than say 128, shared memory might not fit. In that case, we'd need to chunk K, which brings us back to the iterative approach.
Let me think about a hybrid approach:
- For small K (≤ 64), preload everything. Use dynamic shared memory.
- For large K, chunk it.
Actually, the iterative outer product approach with VecK=4 only needs shared memory of TileM*VecK + VecK*TileN = 128*4 + 4*128 = 1024 floats = 4KB. This works for any K. And it's simpler to implement with static shared memory.
Given the tradeoff, I think the iterative approach with 16 syncthreads is still good. But to reduce syncthreads, we can use double buffering with static shared memory:
```cuda
__shared__ float A_smem[2][TileM][VecK];
__shared__ float B_smem[2][VecK][TileN];
```
Then:
- Load into buffer 0
- syncthreads()
- for k = 0; k < K; k += VecK:
- if k + VecK < K: load next into buffer (k/VecK + 1) % 2
- compute from buffer k/VecK % 2
- syncthreads()
Wait, this requires overlapping loads with compute. But loads are in separate code from compute within the same kernel. Without async copies, we can't truly overlap. However, we can at least remove one syncthreads by doing:
```cuda
// Load buffer 0
load_buffer(0, 0);
for (int k = 0; k < K; k += VecK) {
__syncthreads();
// Next load (prefetch)
if (k + VecK < K) {
load_buffer((k/VecK + 1) % 2, k + VecK);
}
compute_from_buffer(k/VecK % 2);
}
__syncthreads(); // only if we wrote to shared memory after compute? No.
```
Wait, this isn't right. We need to ensure buffer is loaded before compute. Let's think again:
```cuda
// Initially load first chunk
load_to_smem(0, 0);
__syncthreads();
for (int k_base = 0; k_base < K; k_base += VecK) {
// Compute from current buffer
compute(smem_current);
// Load next chunk into other buffer (after ensuring previous compute done)
if (k_base + VecK < K) {
__syncthreads(); // Ensure compute done before writing to other buffer
load_to_smem(smem_other, k_base + VecK);
__syncthreads(); // Ensure load done before next compute
}
}
```
This still has 2 syncthreads per iteration (except last). No improvement.
To truly double buffer without extra syncthreads, we need async copies or have loads and computes in different warps. That's a software pipeline, which is complex.
Given the complexity, I'll stick with the simple iterative approach. 16 syncthreads for 2048 FMAs is probably fine. The occupancy should be good enough to hide the barrier latency.
Actually, I just realized: `__syncthreads()` is a block-level barrier. On modern GPUs, it's very fast when all threads reach it quickly. With a well-balanced workload and no divergence, it takes maybe 10-20 cycles. 16 × 20 cycles = 320 cycles. At 2 GHz, that's 0.16 microseconds. For a block doing 2048 FMAs per thread × 256 threads = 524k FMAs, the compute time is much larger. So syncthreads overhead is negligible.
Wait, but there are also loads. 512 loads per block per iteration × 8 iterations = 4096 loads. At one load per 4 cycles (roughly), that's 16k cycles. Plus 524k FMAs at 2 per cycle = 262k cycles. Total maybe ~280k cycles per block. Syncthreads is ~1% overhead.
OK, the iterative approach is fine. Let me finalize.
I want to try one more thing: use `constexpr int TileM = 64, TileN = 256`. This gives fewer blocks (65536 vs 16384? No, (32768/64)*(32768/256) = 512×128 = 65536. Same as 128×128 actually: 256×256=65536. Wait, 32768/128=256, so 256×256=65536. Yes same block count).
What about TileM=64, TileN=128, Block=(32,8)?
- ItemsPerThreadM=8, ItemsPerThreadN=4. Accum=32.
- Very low register pressure. High occupancy.
- Blocks = 512 × 256 = 131072. More blocks, smaller tiles.
- Might be better for load balancing.
I'll go with a moderate config: TileM=128, TileN=128 or maybe try both. But for the first version, let me start with 128×128.
Actually, let me parameterize the main kernel or just hardcode values and try different ones via compile-time constants. Since this is a single Python file with inline CUDA, I'll hardcode a config.
One more thought: since M and N are both 32768, the grid is square. But A is (32768, 32) and B is (32, 32768). In the C output, M dimension corresponds to A rows. We have 32768 rows of A, each 32 floats. For the outer product, processing M in chunks of 128 and N in chunks of 128 makes sense.
Wait, is there a possibility that the test could pass A as (K, M) and B as (K, N)? The forward docstring says:
"A (torch.Tensor): Input matrix of shape (M, K) or (K, M) where M >> N or N >> M.
B (torch.Tensor): Input matrix of shape (K, N) or (N, K) where M >> N or N >> M."
In the get_inputs(), A is (M, N) = (32768, 32) and B is (N, M) = (32, 32768).
So A is (M, K) and B is (K, N). The matmul is straightforward.
What if they pass transposed inputs? The kernel assumes A is (M, K) and B is (K, N). For the general case, we'd need to handle transposes, but the reference just calls torch.matmul which handles all cases.
Since the test in get_inputs() specifically uses M×N and N×M, I'll assume A is (M,K) and B is (K,N) with K being the small dimension. But to be safe, we might want to check and transpose if needed.
Actually, torch.matmul(A, B) for A(M,K) and B(K,N) is straightforward. The dtype is float32.
Let me write the code now. I'll create the load_inline CUDA extension.
One more thing: I should use `-O3` in the extra CUDA flags. And maybe `--ptxas-options=-v` to see register usage, but that would just add output without helping performance.
Actually, I might want to add `-maxrregcount=128` or similar to constrain registers and increase occupancy. But if I use too few, the compiler might spill. With 64 accumulators, 128 is probably safe. Let me not add it initially and see.
Here's the plan for the Python code:
```python
import torch
import torch.nn as nn
from torch.utils.cpp_extension import load_inline
cuda_source = r'''
#include <torch/extension.h>
#include <cuda_runtime.h>
// ... kernel code ...
torch::Tensor forward_cuda(torch::Tensor A, torch::Tensor B) {
int M = A.size(0);
int K = A.size(1);
int N = B.size(1);
auto C = torch::empty({M, N}, A.options());
dim3 block(32, 8);
dim3 grid((M + 127) / 128, (N + 127) / 128);
gemm_kernel<<<grid, block>>>(
A.data_ptr<float>(), B.data_ptr<float>(), C.data_ptr<float>(),
M, K, N);
return C;
}
'''
cpp_source = r'''
torch::Tensor forward_cuda(torch::Tensor A, torch::Tensor B);
'''
functions = ["forward_cuda"]
# In ModelNew:
class ModelNew(nn.Module):
def __init__(self):
super().__init__()
def forward(self, A, B):
return forward_cuda(A, B)
```
Wait, I should also check that A and B are contiguous and on CUDA. The problem doesn't say they're on CUDA by default, but for the benchmark they probably are. In get_inputs, they use `torch.rand` which creates CPU tensors. But the test framework might move them to GPU before calling forward. I should handle this.
Actually, looking at similar problems, the test framework typically moves inputs to the same device as the model parameters. Since ModelNew has no parameters, we should check the device of inputs and handle CPU/CUDA.
If the inputs are on CPU, we can either move them to CUDA and compute there, or fall back to torch.matmul. But the point is to write a CUDA kernel, so we should handle GPU tensors. If inputs are CPU, calling .cuda() on them inside forward might be slow but correct.
Let me make it robust: check if A is on CUDA. If not, move it. Or just do:
```python
if not A.is_cuda:
return torch.matmul(A, B)
```
Actually, for the correctness check, they probably test with the same inputs as get_inputs() which are CPU tensors. But the timing is done on GPU. So I should handle both.
Let me also think about whether my kernel supports arbitrary M, N, K. The grid uses (M+127)/128 and (N+127)/128. In the kernel, I have bounds checks for i < M and j < N. And k < K. Good.
What about K > 32? The kernel handles it with the outer loop. Good.
One thing: `A_smem[TileM][4]` where the 4 is VecK. But what if K < 4? The load loop should guard against it. And in compute, the kk loop should guard.
Let me update the load and compute to handle arbitrary K properly.
For the load:
```cuda
for (int idx = tid; idx < TileM * VecK; idx += BlockDimX * BlockDimY) {
int i = idx / VecK;
int kk = idx % VecK;
if (block_i + i < M && k_base + kk < K) {
A_smem[i][kk] = A[(block_i + i) * K + k_base + kk];
} else {
A_smem[i][kk] = 0.0f;
}
}
```
For B:
```cuda
for (int idx = tid; idx < VecK * TileN; idx += BlockDimX * BlockDimY) {
int kk = idx / TileN;
int j = idx % TileN;
if (k_base + kk < K && block_j + j < N) {
B_smem[kk][j] = B[(k_base + kk) * N + block_j + j];
} else {
B_smem[kk][j] = 0.0f;
}
}
```
For compute: the `if (k_base + kk < K)` check in the inner loop might hurt performance due to divergence. Since K is typically a multiple of 4 and known, we could avoid it. But for correctness I'll add it outside the hottest loop:
```cuda
int valid_k = min(VecK, K - k_base);
#pragma unroll
for (int kk = 0; kk < VecK; ++kk) {
if (kk < valid_k) {
#pragma unroll
for (int ii = 0; ii < ItemsPerThreadM; ++ii) {
int i = ii * BlockDimY + threadIdx.y;
float a = A_smem[i][kk];
#pragma unroll
for (int jj = 0; jj < ItemsPerThreadN; ++jj) {
int j = jj * BlockDimX + threadIdx.x;
float b = B_smem[kk][j];
accum[ii][jj] += a * b;
}
}
}
}
```
Or we can just rely on the loads padding with 0 when out of bounds. If k_base + kk >= K, we loaded 0 into A_smem and B_smem. Then a * b = 0, so accum += 0. No need for the if check in compute! The loads already handled it. This is much better - no divergence in compute.
So the key is: in the load loops, always write to A_smem and B_smem, but use 0 if out of bounds. Then the compute loop is pure FMAs with no branches.
Great! Let me write the final kernel.
Wait, I also need to make sure that A and B are contiguous in memory and on the CUDA device. If not, we can either transpose/copy or fall back.
Actually, let me think about whether to add an early exit for non-contiguous tensors. For the test inputs, they are standard dense tensors so they'll be contiguous.
One more optimization: use `__ldg` for global memory loads. `__ldg` is the read-only data cache load. It can help with caching. Since A and B are read-only in the kernel, we can use `const float* __restrict__` which the compiler might optimize with LDG automatically.
Actually, `__restrict__` is already used. The compiler usually generates LDG for global loads into shared memory.
Let me now draft the full Python module and try to compile it.
I want to make sure the CUDA code has no syntax errors. I'll be very careful.
Also, should I handle both fp32 and fp16? The problem says fp32 precision is required. The reference runs at fp32 (default dtype). So just float is fine.
Wait, `torch.rand` creates float32 by default. Yes.
Let me write it.