reasoning
The first version passed correctness. Now I need to think about performance.
The problem size is:
- B=32
- H=32
- N=512
- D=1024
With my custom kernel:
- BQ=8, BK=16
- 512/8 = 64 query tiles per head
- 512/16 = 32 key tiles per head
- Each block processes BQ=8 queries
- BlockDim=256 (8 warps)
- Grid = (32, 32) = 1024 blocks total
- Each block: D=1024, dt=64, so 16 iterations for dot-product, 16 for V update
- Shared mem: approx 8*(1024+4)*4 = ~32KB for smem_o alone, plus others. Total around 40KB.
This might be OK but there are some concerns:
1. Only 8 warps, meaning occupancy might be OK but threads aren't super utilized
2. The smem_o is 8x1024 = 8192 floats, that's a lot of shared memory. At 1024 blocks total, with shared mem per block ~40KB, that's about 40MB total dynamic, but occupancy is limited to maybe 2 blocks per SM = 2*40KB = 80KB per SM, which is fine (limit is 227KB).
3. But the main issue: the memory access pattern. Each block loads Q once and reuses across all key tiles. K and V are loaded N/BK = 32 times each. Total memory traffic is similar to standard attention.
However, PyTorch's scaled_dot_product_attention is likely using cuDNN/FlashAttention on H100, which is highly optimized with Tensor Cores (even in FP32 it might use TF32 tensor cores or other optimizations). My custom kernel is using only FP32 scalar operations, so I'm unlikely to beat cuDNN unless I also use Tensor Cores.
Actually, wait. The reference uses `torch.nn.functional.scaled_dot_product_attention`. On H100, this will dispatch to FlashAttention if available, cuDNN SDPA, or efficient attention. That will be VERY hard to beat with a simple CUDA kernel, especially without Tensor Cores.
But I can try to use shared memory more effectively and perhaps optimize the kernel for this specific shape. Let me think...
Actually, for fp32 on H100, FlashAttention might not be called (it often uses bf16/fp16). For fp32, it might use the math backend or cuDNN attention. But cuDNN attention for fp32 should still be quite fast.
Wait, the task says "Run faster than the PyTorch reference." It might be possible that the reference is using a non-optimal backend for fp32, or that I can beat it by doing something clever.
Actually, looking at the problem size again:
- Q/K/V size: 32 * 32 * 512 * 1024 * 4 bytes = 2GB
- That's a lot of data
- The attention computation is memory bound
For memory bound operations, I need to minimize memory traffic. Floppy attention minimizes memory traffic by focusing on SRAM. But the problem size (D=1024, N=512) means the standard algorithm might have high memory traffic.
Wait, N=512 is relatively small. D=1024 is large.
- Each Q row has 1024 floats = 4KB
- For all 512 queries in a head: 2MB
- Each head processes this...
Actually, with N=512 and D=1024, for standard attention:
- Q·K^T: each query dot products with all keys.
- Memory reads per head: Q (2MB) + K (2MB) + V (2MB) + output (2MB) = 8MB if we compute exactly once.
- But for the matrix multiply, we read K N/BK times or similar.
For FlashAttention style, the main reduction is in writes to O: only one write per element instead of multiple writes through softmax. But for D=1024, the intermediate activations might be harder to fit.
With my kernel, each block loads:
- Q: 8 rows * 1024 cols = 8K floats = 32KB (loaded once per block)
- K: 16 rows * 1024 cols = 16K floats = 64KB (loaded 32 times) = 2048KB
- V: same as K, 2048KB
- O: written once, 32KB
Total traffic per block = ~4.1MB
1024 blocks * 4.1MB = 4.2GB total memory traffic.
Wait, but K/V are shared across all query groups within a head? Actually no, in my kernel each block only processes 8 queries. So for all 64 query groups in a head, each reads the full K and V. So total K reads per head = 64 * 2MB = 128MB. Similarly V reads = 128MB. Q reads = 2MB once. Output = 2MB once.
Total = 260MB per head * 1024 heads = 260GB? That can't be right.
Let me recalculate:
B=32, H=32, so total heads = 1024.
Per head: Q = 512*1024*4 = 2MB, K = 2MB, V = 2MB, O = 2MB.
In my kernel, per block (8 queries):
- Q: 8*1024*4 = 32KB (loaded once)
- K: 2MB * (512/16) iterations? Actually K is loaded tile by tile: 16*1024*4 = 64KB per tile, and there are 32 tiles, so 2MB total K reads per block.
- V: same 2MB total reads per block.
- O: 32KB written once.
Total per block = 4.1MB.
Total for all blocks (1024) = 4.2GB.
Wait, but K and V are the same across all query blocks! So my kernel is reading K and V way too many times. In a well-optimized attention kernel, you'd want to load K and V once and compute dot products with multiple query tiles, but with D=1024 and only ~227KB shared memory, and 1024 features, you can't fit a whole K/V head in shared memory.
Actually, wait. The standard flash attention loads k tiles and keeps them in shared memory while iterating over query tiles. Wait, no - the standard FA loads Q tiles and keeps them in register, then iterates over K/V tiles. Or it loads Q tiles in SRAM and iterates over K/V tiles.
But optimized implementations use fused MMAs and load K/V through L1 cache, benefiting from cache hits when multiple blocks process different query tiles but load the same K/V tiles. On H100, with:
- 1024 blocks total
- 132 SMs
- 1024/132 ≈ 7.7 blocks per SM
But each block reads the full K and V for its head. Each SM might process blocks from different heads (with different K/V), so L1 cache effectiveness depends on whether blocks from the same head are co-located on the same SM.
Actually, with a grid of (B, H) = (32, 32), blocks are launched in x-major order. All 32 blocks for b=0, h=0 to h=31, then b=1, etc. Or maybe y-major? blockIdx.x is fastest, so blocks (0,0), (1,0), ..., (31,0), (0,1), ..., (31,31). So within a wave of 132 SMs, we process all 32 batches for h=0, then some for h=1. So blocks from the same head may be quite spread out.
Hmm, this is hard to optimize.
Alternative approach: for N=512, D=1024, maybe I can transpose the loops? Let me think about using Tensor Cores...
Actually, H100 has very fast Tensor Cores. Even for FP32, using Tensor Cores properly could be beneficial. But writing a Tensor Core kernel from scratch for attention in C++ CUDA is very complex.
Wait, let me think about the reference. The PyTorch SDPA reference on H100 for FP32 might not even use Tensor Cores. Let me check what backends are available:
- For FP32, SDPA often falls back to the math backend (plain cuBLAS) if `is_causal=False`, or it might use cuDNN attention.
Actually, on H100, anything built into PyTorch would be well optimized. Beating it with a simple hand-written kernel is hard, but not impossible if the reference uses a less optimized backend for some reason.
Let me reconsider the kernel. Maybe I should use a simpler approach: since N=512 and D=1024, I could compute S = Q @ K^T (32*32 blocks of 512x512), then softmax, then S @ V.
For the matmuls I could call cuBLAS from the kernel? No, from cuda source I can call cublas if I link it, but load_inline might not link cublasSgemm by default. Actually, I can specify extra libraries.
But the kernel overhead of separate kernel calls for matmul + softmax + matmul might be high. Still, cuBLAS matmuls are extremely optimized.
Alternative: use a single kernel but write it better. Let me think about using more threads and better memory patterns.
What if I use a warp-level approach where each warp computes one output element? For D=1024, with 32 threads per warp, each thread does 1024/32 = 32 MACs. That's feasible.
Actually, a common approach for efficient attention on small sequence lengths is to have each thread compute a single attention output. But with D=1024, we might need multiple queries per block.
Let me try a different approach: use `mma` via Tensor Cores. But that's very complex for inline CUDA.
Alternatively, I can use `wmma` or the newer `mma.sync.aligned` PTX. With Hopper, there are specialized warpgroup-level matrix multiply instructions. But again, this is extremely complex for a single file submission.
Wait, can I just use cuBLAS within my extension? Yes! I can create a CUDA extension that calls cuBLAS for the matrix multiplications, and write custom CUDA for softmax. This would likely be very fast.
Let me think:
1. Compute S = Q @ K^T for all heads: this is a batched matrix multiply of [BH, M, K] @ [BH, N, K]^T -> [BH, M, N]
2. Apply softmax row-wise
3. Compute O = S @ V: [BH, M, N] @ [BH, N, K] -> [BH, M, K]
For Q, K, V shapes: [32, 32, 512, 1024]
We can reshape to [32*32, 512, 1024] = [1024, 512, 1024]
Then S = [1024, 512, 512]
O = [1024, 512, 1024]
Memory usage:
- Q: 32*32*512*1024*4 = 2GB
- K: 2GB
- V: 2GB
- S: 32*32*512*512*4 = 1GB
Total intermediate = 1GB. That's a lot but H100 has 80GB so it's fine.
The overhead is one extra 1GB read/write. But cuBLAS GEMM is extremely fast on H100, likely saturating memory bandwidth or compute.
However, PyTorch's SDPA probably does the exact same thing: call cuBLAS for the matmuls, but fused with the softmax in between. Actually, SDPA likely uses the "flash attention" algorithm where intermediate S isn't materialized to global memory.
If I materialize S to global memory, I'll be doing 1GB extra read/writes. That might hurt performance.
But wait - on H100, with 3.35 TB/s bandwidth:
- Reading Q+K+V once: 6GB
- Writing O once: 2GB
- Extra S read/write: 2GB
Total: 10GB vs 8GB for fused.
4GB * 1024/3.35TB/s ≈ 1.2 ms difference in bandwidth-bound time. But the compute savings might help.
Actually, PyTorch SDPA for FP32... let me check if it might be using efficient_attention kernel from xFormers or something. On recent PyTorch with H100, `scaled_dot_product_attention` with FP32 and no mask might call cuDNN attention or flash attention.
If it uses cuDNN attention, it's extremely hard to beat.
However, maybe I can write a better fused kernel using proper attention algorithm but with better parallelism.
Let me reconsider my first kernel. The issue is:
1. Only 8 warps per block
2. warp_id = tid/32, so warp_id can be 0..7
3. q_local = warp_id (0..7) which maps to one of the BQ=8 queries
4. lane_id < BK=16 means only 16 threads per warp are used
5. The remaining 16 threads are wasted for all score computations and softmax
6. That's only 50% thread utilization for the critical path
This is a major issue. Let me redesign to use all 32 threads.
New design:
- Process 32 queries per block? Or process 16 queries but use all lanes?
- With D=1024, let's keep BQ and BK smaller.
Alternative approach inspired by FlashAttention:
- Load Q tile into smem (or registers)
- Iterate over KV tiles
- Use all warps to parallelize over queries and keys
What if:
- Each block processes BQ queries (e.g., 16)
- Each warp processes one query? No, 16 warps = 512 threads, which might be too many.
- Or each warp processes part of the score matrix.
Actually, for online softmax + matrix multiply, a common efficient approach:
- Load Q tile: 16 rows x 1024 cols
- Each thread computes partial dot products for one (q, k) pair
- With 256 threads, 16 queries, 512 keys -> 16*512 = 8192 elements in score matrix. 8192/256 = 32 elements per thread. Each element needs 1024 MACs. 32*1024 = 32768 MACs per thread. At 50 TMAC/s on FP32... that's too much for one thread!
Wait, I mean each key tile is 16 keys (BK=16). Each iteration of the outer loop processes one KV tile.
Let me try:
- BQ = 16, BK = 16
- Block = 256 threads = 8 warps
- Each warp handles 2 queries? Actually let's reconsider.
For score computation:
- We need S = Q_tile @ K_tile^T where Q_tile is BQ x D and K_tile is BK x D
- Result is BQ x BK = 256 elements
- With 256 threads, each thread computes exactly one dot product of D=1024... that's 1024 MADs per thread per tile, very slow!
Alternative: each thread computes partial dot product, then reduce.
- Each thread handles 1 element in D dimension
- For a (q, k) pair: 1024 threads would each multiply one pair, then reduce with warp shuffle / shared memory.
With 256 threads, we'd need 1024/256 = 4 accumulations per thread for each (q, k). But there are 256 (q, k) pairs. Total parallel work = 256 pairs * 1024 elements = 262k multiplications per tile. With 256 threads, each thread does 1024 multiplications, then we need reductions.
Wait, here's a better way:
- Load Q_tile and K_tile into shared memory
- For each query in Q_tile:
- Have all threads compute its dot product with all keys in K_tile in parallel
- Specifically, each warp computes scores for 1 query across 32 keys? But BK=16, so only 16 keys.
Actually, let me think differently. Let's transpose: have BQ=64 but each query processed by fewer threads? No...
What about:
- BQ=32 queries
- BK=32 keys
- Block=256 threads
- D=1024
- For dot product, use all 256 threads to compute partial dot products for 32 queries:
- Each thread handles D/256 = 4 features for all queries
- Wait, this doesn't parallelize across K.
Actually, for a block of 32 queries x 32 keys x 1024 dims:
- Total work: 32*32*1024 = 1M multiplications
- 1M / 256 threads = 4096 mults per thread per tile
- Number of key tiles: 512/32 = 16
- Total mults per thread: 65536
- At say 20 Tops FP32 on normal CUDA cores (not tensor cores), that's about 3.3 microseconds of compute per head... negligible? Actually 65536 mult-adds / (67e12 FLOPS / 132 SMs) ≈ very small. The bottleneck is memory.
Hmm wait, 67 TFLOPS is the full chip FP32 throughput. 67e12 ops/sec. 65536 ops is 1 nanosecond. But memory is much slower. So the kernel is definitely memory bound.
For memory-bound attention with D=1024, N=512:
- Total data movement for fused kernel = roughly: load Q once, load K N/BK times... no, load K once per key tile per query tile.
- In Flash Attention, each query tile loads all K tiles. To reduce this, you'd want larger query tiles to amortize K/V reads.
With my kernel:
- BQ=8, so K is read 64 times per head
- Total K reads = 64 * 2MB = 128MB per head
- Total for all heads = 128MB * 1024 = 131GB
- Similarly for V
- Q reads: 2MB * 1024 = 2GB
- O writes: 2GB
- Total traffic ≈ 264GB
- At 3.35TB/s = 3350GB/s, time = 79ms. That's way too high!
Wait, but that's not right because:
1. K and V are cached between query tiles if multiple blocks from the same head run on the same SM
2. Also, with HBM3 on H100, maybe actual sustained bandwidth is different
3. But 264GB for one forward pass on 2GB inputs seems plausible for naive attention
Let me recalculate more carefully. The standard attention complexity is O(N^2*D). The memory operations if we materialize nothing:
- Q: read once, K: read N/BK times per query block... but in theory, if we have enough cache, each K is read once and stays in L2.
Actually, H100 has 50MB L2 cache. Per head, K = 2MB. So we can fit about 25 heads in L2. With 1024 heads total and ~132 concurrent blocks, if each block is from a different head, L2 caching doesn't help.
But if multiple blocks from the same head are running concurrently on different SMs, they all compete for L2 cache but each reads the same K, so L2 can serve it without going to HBM after the first read.
Actually, in a wave of 132 blocks, if they're all from different heads, then K isn't shared. But the next wave will also be from different heads. So every K read goes to HBM. This is terrible.
PyTorch's FlashAttention would be much better because it:
1. Uses much larger tile sizes, amortizing K/V loads
2. Or uses the online softmax algorithm without materializing S to global memory (which my kernel also does, but PyTorch's does it better)
3. And it uses Tensor Cores, which are much more efficient for the GEMM-like operations
Wait, but if PyTorch SDPA uses cuDNN attention, cuDNN may also materialize some intermediates internally in FP32. Or it might also be highly optimized.
Actually, wait. If PyTorch uses SDPA for FP32, it might actually dispatch to the `math` backend or a less optimized kernel because FP32 flash attention often isn't available. Let me check:
In PyTorch, `scaled_dot_product_attention` has the following backends:
- `flash_attention` (typically FP16/BF16 only, not FP32)
- `mem_eff_attention` (FP16/BF16/FP32, but for FP32 it might be less optimized)
- `cudnn_attention` (available for FP16/BF16 on Ampere+, maybe FP32 too?)
- `math` (always available)
For FP32 input on H100, it will likely try cuDNN first, then memory efficient attention, then math backend. cuDNN attention for FP32 might still be available.
If cuDNN attention is used, it's almost impossible to beat.
But wait, the correctness test passed. Let me think about whether I can do something simpler that actually performs better.
Actually, let me think about this more carefully. The trick in these kinds of challenges is often that the reference implementation doesn't use the best available backend because of some configuration issue, or because FP32 SDPA doesn't dispatch to the fastest path.
For instance, if cuDNN attention requires FP16/BF16 for its optimized path, then the FP32 reference might fall back to `memory_efficient_attention` or even `math`. The `math` backend computes `softmax(Q @ K^T / sqrt(d)) @ V` using explicit PyTorch operations, which involves materializing the N x N attention matrix to memory.
If the reference materializes S = Q @ K^T (1GB), then does softmax, then O = S @ V, that would be:
1. S = matmul(Q, K^T): read Q (2GB) + K (2GB) + write S (1GB) = 5GB
2. softmax(S): read S (1GB) + write S (1GB) = 2GB
3. O = matmul(S, V): read S (1GB) + V (2GB) + write O (2GB) = 5GB
Total: 12GB
My fused kernel avoids materializing S, so even with poor reuse, it might be faster if the reference does materialize S.
But if the reference uses flash attention (even in FP32 with some support), it would be very fast.
Hmm, let me try my kernel and see if it's actually faster. But I can't run timing directly - I can only submit. However, I can reason about it.
Actually, I realize my kernel has a bug in the soft max update with corr_prev. In the kernel I wrote:
```cpp
float p_weight = exp_val * corr_tile;
```
Wait, in online softmax, when updating O accumulator, we should use:
```cpp
float corr_prev = expf(m_prev - m_new);
float p_weight = expf(my_s - m_new);
```
Currently I compute:
```cpp
float p_weight = exp_val * corr_tile;
```
where `exp_val = expf(my_s - m_val)` and `corr_tile = expf(m_val - m_new)`. So `p_weight = expf(my_s - m_new)`. That's correct! And `corr_prev = expf(m_prev - m_new)`. And O is updated with `corr_prev * old + sum(p_weight * V)`. That's correct.
OK, the kernel is numerically correct. Now about performance...
Actually, I just realized the kernel submission passed correctness, but it's probably quite slow. Let me try to optimize it.
Main optimization ideas:
1. Use all 32 threads in each warp effectively
2. Increase tile sizes to reduce key tile count (fewer K/V loads)
3. Consider coalescing and vectorized loads
For N=512, D=1024:
- If we can load 32 keys at a time (BK=32), that's only 16 key tiles instead of 32. This halves K/V reads.
- But shared memory for K: 32*1024*4 = 128KB, which fits in shared memory (227KB limit).
- Q tile: 8*1024*4 = 32KB.
- O tile: 8*1024*4 = 32KB.
- p tile: 8*32*4 = 1KB.
- Total: ~200KB, under the 227KB limit.
Let's try BK=32.
With BK=32 and block=256 threads (8 warps):
- Each warp should process some queries. With BQ=8, we have exactly 8 warps, so each warp handles 1 query.
- For score computation with BK=32, only lane_id < 32... wait, all 32 lanes are used since BK=32. Great!
- For the dot product, each thread handles part of D=1024. With DT=64 and 2 parts (lane_id/BK = lane_id/32, that's 0 or not... wait, BK=32 means lane_id/32 = 0 for all lanes. So we need a different way to split D.
Actually, with BQ=8, BK=32, block=256:
- There are 8 warps.
- Each warp handles 1 query, and all 32 lanes handle the 32 keys.
- For dot product over D=1024:
- We need to split the 1024 features among threads.
- One approach: each thread computes 1024/32 = 32 features for its 1 query x 1 key pair? No, there are 32 keys, but each lane only handles 1 key.
Wait, each warp handles 1 query (q_local = warp_id) and 32 keys (k_local = lane_id). So each warp computes 32 dot products, one per lane. Each dot product is D=1024.
To parallelize the dot product, we can have each lane process only part of D, then shuffle reduce. For example:
- load Q for query q_local into shared memory
- load K for all 32 keys into shared memory
- each lane computes its own dot product for key k_local
- To parallelize within lane: each lane does partial dot product over a chunk of D
- But a lane is a single thread. We can load D in chunks of DT=64.
- Lane k_local accesses K[k_local][d] and multiplies by Q_shared[q_local][d].
- This is just a scalar loop over D.
Wait, but each lane is computing its own dot product serially over D=1024. That's 1024 iterations, each doing one multiply-add. With ~32 FMAs per clock per thread on H100, and ~2GHz clock, that's about 16 cycles per FMA (throughput limited). Actually FP32 FMA throughput is 1 per thread per clock on Ampere/Hopper for simple FMAs (when not using tensor cores). At 2GHz, 1024 FMAs take ~500 ns. That's quite fast.
Then with 16 key tiles (512/32), each warp/ thread does 16*1024 = 16384 FMAs. At 2GHz, ~8 microseconds. Not too bad.
For the V update:
- Each query needs O = P @ V where P is 512 keys, V is 512x1024.
- With BK=32, we process 32 keys at a time.
- For each V tile, we need to multiply P_tile (1x32) by V_tile (32x1024) and add to accumulator.
- Our warp is computing for 1 query.
- Lane k handles key k. P[q][k] is loaded.
- Then for each output dimension d, we accumulate P[q][k] * V[k][d].
- With 32 lanes and D=1024, we need all lanes to help. But a single warp only has 32 threads.
Hmm, we could have lane d (for d < 32) accumulate over k=0..31 for V[k][d] and then process multiple chunks of D. For D=1024, each lane handles 1024/32 = 32 dims. So lane i accumulates for dim i + 32*j where j=0..31.
Wait, that's a lot of unrolling. Let me think:
- Load P[q][k] for k=0..31 (in shared memory from earlier step)
- Load V[0..31][0..1023] into shared memory (or directly access)
- Each of the 32 lanes computes a partial sum for a block of 32 dims.
- Specifically, lane d computes sum over k of P[q][k] * V[k][d_base + d] for a d_base.
With D=1024 and 32 lanes, we can process 32 dims at a time.
- iteration 0: dims 0..31, lane d handles dim d
- iteration 1: dims 32..63, lane d handles dim d+32
- ...
- Total: 1024/32 = 32 iterations.
In each iteration, each lane does 32 multiply-adds (over k=0..31). So 32*32 = 1024 FMAs in 32 iterations.
Again, 16 key tiles * 1024 FMAs = 16384 FMAs. Plus score computation: 16384 FMAs. Total ~32768 FMAs per warp per head.
8 warps * 32768 = 262k FMAs per block.
Total blocks = 32*32*64 = 65536 blocks? Wait: B=32, H=32, N/BQ = 512/8 = 64. Total blocks = 32*32*64 = 65536.
65536 blocks * 262k FMAs = 17.2e9 FMAs total.
At 67 TFLOPS, that's 257ms. But memory bandwidth might dominate.
Actually, 65536 blocks is a lot. At 256 threads each, that's 16.7M threads. With 132 SMs at ~2000 concurrent threads each, occupancy is fine.
But 65536 blocks means each K tile is loaded 65536 times. Hmm, no: for each head, there are 64 query blocks, and each loads all 16 key tiles. So per head: 64*16 = 1024 K tile loads. Total heads = 1024. Total K tile loads = 1,048,576. Each tile is 32*1024*4 = 131KB. Total K traffic alone = 1,048,576 * 131KB = 137GB. That's massive.
137GB / 3.35TB/s = 41ms just for K reads. That's actually not terrible, but there's also V reads and Q reads. So total time could be 80-100ms.
Wait, but the reference might be doing something similar or worse.
Actually, let me reconsider whether PyTorch's scaled_dot_product_attention uses FlashAttention for FP32. Looking at PyTorch code, FlashAttention in PyTorch typically only supports FP16 and BF16. For FP32, it likely falls back to either memory-efficient attention or the math backend.
Memory-efficient attention (the xFormers-style one) does fuse the operations to avoid materializing S, but might still use a similar algorithm to mine, possibly worse.
Let me actually check: can I write a kernel that's competitive?
Given the constraints, maybe the simplest winning strategy is:
1. Keep the current fused kernel approach (no materialized S)
2. Optimize thread utilization and tile sizes
3. Try to make it faster
But there's another trick: since D=1024 and N=512, I can fit the entire S matrix for a query tile in shared memory or registers? S for BQ queries and all N keys would be BQ x N = 8 x 512 = 4096 floats = 16KB. That's small.
Actually, what if each thread block processes ALL queries for a head at once?
- BQ=512, but we process them in chunks for shared memory.
- Q for all 512 queries: 2MB. Too big for smem.
- But we can load Q to registers one query at a time.
Alternative: Each thread computes one complete attention output vector.
- There are B*H*N = 32*32*512 = 524,288 output vectors.
- Each output vector is D=1024 floats.
- Total threads can work in parallel.
- Each thread computes softmax attention weights for one query.
- Steps for one thread (one query):
1. Load Q vector (1024 floats) into registers/local memory.
2. Iterate over all keys (512 keys):
- Load K vector (1024 floats).
- Compute dot product Q·K (1024 mult-adds).
- Track max and exp-sum for softmax online.
3. Iterate over all keys again:
- Load K and V vectors.
- Compute weight = exp(Q·K - max) / sum
- Accumulate weight * V to local output buffer.
4. Write output.
This requires loading K twice and V once, but avoids any inter-thread synchronization.
- It also avoids shared memory entirely.
Each thread does:
- Read Q: 1024 floats
- Read K: 512 * 1024 floats
- Read V: 512 * 1024 floats
- Write O: 1024 floats
Total per thread = 1025K floats = 4.1MB
Total for all threads = 524288 * 4.1MB = 2.1TB. That's way too much global memory traffic.
But wait, when multiple threads access the same K/V from global memory, the GPU coalesces and caches. If all threads in a warp access the same key (for different queries or different features), that would be broadcast. But if they access different keys, it's coalesced if they're contiguous in memory.
Actually, let's organize threads differently:
- A warp of 32 threads processes 32 different queries, and they all load the SAME key to compute dot products with their respective Q vectors.
- For key k, warp loads K[k][:1024] (broadcast to all threads).
- Thread t has Q[t][:1024] in registers or local mem.
- Each thread computes Q[t]·K[k] (1024 FMAs).
- After 32 queries x 512 keys, each thread has its attention scores.
- Then for V accumulation: load V[k][:] for key k, each thread multiplies by its softmax weight and accumulates.
This way:
- K is loaded once per key per warp (32 queries), and broadcast.
- V is loaded once per key per warp.
- Q is loaded once per query per warp.
- Memory traffic per 32 queries: Q(32*1024) + K(512*1024) + V(512*1024) + O(32*1024) floats = ~1056K floats = 4.2MB
- Total warps = 524288/32 = 16384
- Total traffic = 16384 * 4.2MB = 68.8GB
That's much better than my shared-memory kernel! 68.8GB is very reasonable. At 3.35TB/s = 3350GB/s, that's about 20.5ms just for memory bandwidth. Add compute, maybe 25-30ms.
And every K/V access is perfectly broadcast across a warp, maximizing cache utilization. Each key is loaded once from HBM and then broadcast to 32 threads via L1/L2.
But 32 threads computing 32 dot products simultaneously, each against the SAME K vector:
- Each thread processes one query.
- K is broadcast to all threads.
- But Q must be in registers for each thread.
- D=1024, so each thread needs 1024 registers for Q... that's 1024*32 = 32k floats, which is 128KB for the warp. Each SM has 64K registers total, so that's already 2x too much per warp! We can't keep all of Q in registers across all 1024 dims if we need registers for other things too.
We can process Q in chunks. Load Q[d:d+step] and K[d:d+step], compute partial dot product. Then reload for next chunk.
- Each thread keeps an accumulator for the score.
- Register usage per thread: Q_partial (step floats) + score (1 float) + max,sum (2 floats) + V_partial (step floats) + O_partial (step floats)...
- If step=64: Q_partial = 64, V_acc = 64, plus a few scalars ≈ 130 floats ≈ 130 registers. That's reasonable.
So the algorithm per thread block:
- BlockDim = warp_size * num_warps
- In my earlier calculation, one warp handles 32 queries. We can use multiple warps per block.
- Say 8 warps per block = 256 threads, handling 32*8 = 256 queries.
- Batch size = 32, Heads = 32, N = 512. Total queries = 524288.
- 524288/256 = 2048 blocks. Grid = (2048,) or split as (32, 32, 16) for clarity.
For 256 threads (8 warps) in a block:
- Each warp handles 32 queries.
- Warp 0: queries [0..31]
- Warp 1: queries [32..63]
- etc.
Synchronization:
- No need for __syncthreads() within a block if we just want atomic writes or separate outputs. Each warp computes independently.
- But we need to make sure warps can load K and V cooperatively if we want to use shared memory for K/V to reduce global memory traffic.
Wait, with 8 warps loading the same K/V but for different queries, if each warp loads K/V independently, the cache will serve it. On H100, L1 cache is quite good for broadcast. But we could also put K/V in shared memory to guarantee one load per key per block.
Shared memory approach for K/V:
- Key tile size = say 1 key (since we process one key at a time for all warps in the block).
- Load K[k][0:1024] and V[k][0:1024] to shared memory.
- 1024 floats * 2 = 8KB per key.
- With shared memory 227KB, we can load multiple keys.
- But 8 warps * 32 queries each = 256 queries. All processing same key tile.
- If we load 8 keys at a time: 8 * 8KB = 64KB of shared memory.
- 512 keys / 8 = 64 key tiles.
Let's design it:
- Process BK=8 keys at a time.
- Load K[0:8][0:1023] to smem_k[8][1024]
- Load V[0:8][0:1023] to smem_v[8][1024]
- Each warp handles 32 queries.
- For key j in 0..7:
- Load Q for each query into registers (or process in chunks).
Actually, let's process D in chunks of DT=64 to keep register pressure low.
For each query, we do:
1. Score accumulation:
For each key tile (8 keys):
For each chunk d in 0..1023 step 64:
Load Q[q][d:d+64] into registers
Load K[k][d:d+64] from shared memory
Compute partial dot product
Sum partials via shuffle to get score
2. softmax
3. V accumulation:
For each key tile (8 keys):
For each chunk d in 0..1023 step 64:
Load V[k][d:d+64] from smem
Multiply by exp(score-max)/sum
Accumulate into O[q][d:d+64]
But to do this, we need two passes (one for scores, one for V accumulation), which doubles Q reads.
Flash style online softmax avoids this by doing it in one pass. But with warps processing different queries and sharing K/V, online softmax is harder because each thread needs its own running max and sum.
With two passes though:
- Pass 1: compute all scores for all 512 keys, per query. Store 512 scores in shared memory or local memory. 512*4 = 2KB per query. For 256 queries per block = 512KB, too much.
- We can do it in tiles: compute scores for BK keys, do softmax update, compute V contribution for those keys, discard scores. This is the online softmax approach.
Online softmax per query in one pass:
- For each key tile:
- Compute scores and find tile max
- Update running global max
- Update running L and O accumulators
This works perfectly! Each thread maintains its own max, sum, and O accumulator in registers.
Let's formalize the per-thread algorithm:
Each thread handles exactly 1 query (1 out of 256 in a block).
Wait, 256 threads = 8 warps * 32 threads. We have 256 queries in a block, so 1 thread per query. That's perfect!
Algorithm per thread (handling query q):
max_val = -inf
sum_exp = 0
o_acc[DT] = {0} but we'll process in chunks
For key_tile = 0 to N-1 step BK:
// compute scores for this key tile
score = 0
For d_chunk = 0 to D-1 step DT:
load Q[q][d_chunk:d_chunk+DT] from global to registers
For k_local in 0..BK-1:
load K[key_tile+k_local][d_chunk:d_chunk+DT] from shared
Compute partial score...
Wait, each thread needs BK scores (one per key in tile).
If we have 1 thread per query, then for BK=8 keys and D=1024:
- Each thread needs to compute 8 dot products.
- Each dot product = 1024 FMAs.
- 8 * 1024 = 8192 FMAs per thread per tile.
- 512/8 = 64 tiles.
- Total FMAs for scores = 8192 * 64 = 524k FMAs per thread.
- Plus V accumulation: BK * D = 8 * 1024 = 8192 FMAs per tile * 64 = 524k.
- Total: ~1M FMAs per thread.
At 67 TFLOPS for 524288 threads? Actually, the GPU parallelizes across threads. Total work = 524288 threads * 1M FMAs = 5.24e11 FMAs = 524B operations.
At 67 TFLOPS = 67e12 ops/sec, time = 7.8ms just for compute. But we're memory bound so it'll be slower.
Memory traffic:
- Q: Each query loaded once = 524288 * 1024 * 4 = 2GB. Actually Q is 32*32*512*1024*4 = 2GB.
- K: Each key loaded once per block? Wait:
- 256 queries per block, all processing all 512 keys.
- With BK=8 keys per tile, 64 tiles.
- K is loaded 64 * 8 = 512 keys per block, which is the full K.
- 2048 blocks total, but each block loads its head's K.
- Total K traffic = 32*32*512*1024*4 = 2GB.
- V: same 2GB.
- O: 2GB.
Total: 8GB. That's the absolute minimum!
At 3.35TB/s, 8GB = 2.4ms. Plus compute overhead. This could be quite fast!
But wait, with 1 thread per query, and 256 threads per block, how do we load K and V into shared memory cooperatively?
Block organization:
- blockDim.x = 256 threads
- Each block processes 256 queries from one head
- There are 512 queries per head, so 2 blocks per head.
- Total blocks for one head = 2
- For B=32, H=32, total blocks = 32*32*2 = 2048 blocks.
Wait, each head has 512 queries. With 256 queries per block, we need exactly 2 blocks per head.
But each block loads all 512 keys of its head's K and V. So K is loaded twice per head (once per block). That doubles K traffic.
K traffic = 2GB * 2 = 4GB. Similarly V = 4GB. Q = 2GB. O = 2GB. Total = 12GB.
Still quite good! 12GB / 3.35TB/s = 3.6ms.
For 2048 blocks and 132 SMs, that's 15.5 blocks per SM. With 256 threads per block, that's about 4000 threads per SM, which is fine.
But wait, can we fit 256 threads per block? Yes, up to 1024. Can we schedule 15 blocks per SM concurrently? Max is 32 blocks per SM. So we're fine.
However, each block processes 256 queries but only 1 at a time per thread. Let me think about how to structure the kernel:
```cpp
__global__ void attention_kernel(
const float* Q, const float* K, const float* V, float* O,
int B, int H, int N, int D
) {
int tid = threadIdx.x;
int q_per_block = 256; // blockDim.x
int q_block_id = blockIdx.z; // 0 or 1 for N=512
int q = q_block_id * q_per_block + tid;
int b = blockIdx.x;
int h = blockIdx.y;
if (q >= N) return;
// read Q[q][0:D] into registers in chunks
// iterate over all keys
// use shared memory for K and V tiles
int head_offset = ((b * H + h) * N + q) * D;
int head_offset_kv = ((b * H + h) * N) * D;
float scale = 1.0f / sqrtf((float)D);
float m = -1e20f;
float l = 0.0f;
// O accumulator: process in chunks
// We can't keep 1024 floats in registers per thread.
// We'll process D in chunks of, say, 64.
for (int d_base = 0; d_base < D; d_base += DT) {
float o_acc = 0.0f;
// Wait, we need to interleave keys and D chunks simultaneously.
// Standard approach: inner loop over keys, but that means we need Q for each key.
// Better: outer loop over keys, inner loop over D chunks.
// Load Q[q][d_base:d_base+DT] once, then iterate over key tiles.
// But for online softmax, we need scores for ALL keys before updating O?
// No, we can do online softmax at the key level.
// Actually, if we process D in chunks, we need partial scores per key tile.
// This gets complicated.
// Alternative: keep all O partials in a small shared memory per thread,
// but 256 threads * 1024 floats = 1MB, which is too much.
}
}
```
The problem is that for online softmax, the correction factor for O accumulators depends on the global max, which isn't known until all keys of a tile are processed. But we can process all keys for one D chunk at a time? No, that doesn't work because the dot product requires ALL D for each key.
We need to process per key, over all D. But D=1024 is too big for per-thread registers. We could break D into chunks but then we need to solve the online softmax problem.
Alternative 1: Two-pass algorithm for each query.
Pass 1: compute all 512 scores, store in shared memory.
Pass 2: compute softmax, then accumulate V.
But 512 scores per query, and 256 queries per block = 131k floats = 512KB. That fits in shared memory!
Wait, shared memory is 227KB per SM, but per block it's also limited to 227KB. 512KB > 227KB. So we can't store all scores for all 256 queries in shared memory.
But we can do it in chunks: process 64 queries at a time within a block.
- 64 queries * 512 scores = 32768 floats = 128KB. That fits!
- Process 64 queries, store their 512 scores.
- Then compute softmax and accumulate V for those 64 queries.
- Repeat 4 times (for 256 total).
But then we have 64 queries computing scores, while the block loads K and V. The other 192 threads are idle or could help with loads.
Actually, a better approach: have ALL 256 threads work together to load K/V into shared memory, but only 64 threads compute scores at any given time. That's wasteful.
Wait, each thread handles one query. We could have each thread compute its score and immediately keep track of max, then in a second pass compute exp and accumulate. But that requires storing 512 scores per thread in local memory. Local memory is backed by global memory, which is slow. 512 floats = 2KB per thread. With 256 threads, total local memory = 512KB per block. This spills to global memory and causes 512KB of global memory traffic per block, which is acceptable (2048 blocks * 512KB = 1GB extra). But local memory accesses are slow.
Hmm, local memory: each thread accesses its own private area, so they should be coalesced if laid out nicely. But the latency is higher than registers.
Actually, for CUDA, local memory is just global memory with known addressing, so each access is a regular global load/store. 512 floats accessed sequentially should be coalesced per warp.
But if each thread needs to write 512 floats and then read them back, that's 2 * 512 * 4 = 4KB per thread, and with 256 threads running in parallel... Actually no, local memory is allocated per thread and the accesses go through L1/L2.
Alternative: Use a smaller number of threads per block but more blocks. For instance:
- 1 block = 1 warp = 32 threads, processing 32 queries.
- 512 queries / 32 = 16 blocks per head.
- Total blocks = 32 * 32 * 16 = 16384 blocks.
- Score storage: 32 threads * 512 scores = 16384 floats = 64KB per block. This easily fits in shared memory!
- Wait, 64KB < 227KB, yes!
16384 blocks is actually fine. It gives very good parallelism.
Each block processes 32 queries and loads the full K/V for its head.
- K is loaded 16 times per head.
- V is loaded 16 times per head.
- Total K traffic = 32GB, V = 32GB. That's too much!
Wait: 16 blocks per head * 2MB K per head = 32MB per head. 1024 heads = 32GB. Yes, that repeats K and V 16x.
But with L2 cache of 50MB, if each head's K (2MB) can stay in L2 while its 16 blocks run, then only the first block reads K from HBM and the rest use L2. But the blocks might be scattered across SMs.
Actually, with 16384 blocks, and 132 SMs, there are ~124 blocks per SM. These blocks are from different heads. L2 cache can hold K for 25 heads. The remaining 1000 heads' K won't be cached. So most K/V blocks will be HBM reads.
This is bad. We need to reduce per-head K/V reads.
To reduce K/V reads per head, we need FEWER blocks per head that all load the same K. That means MORE queries per block, which brings us back to the shared memory score storage issue.
What if we have 1 block per head that processes ALL 512 queries?
- Block = 512 threads. Each thread handles 1 query.
- 32*32 = 1024 blocks total.
- Each block loads K and V once. Great!
- 512 scores per thread. We need to store them somewhere.
- If we keep them in local memory, each thread writes and reads 512 floats.
- Total local memory traffic = 1024 blocks * 512 threads * 512 floats * 4 bytes * 2 (read+write) = 2GB. Not too bad!
But can we have 512 threads per block? Yes, up to 1024.
512 threads * 512 scores * 4 bytes = 1MB for score storage if using shared memory. But 1MB exceeds 227KB per block.
If we use local memory (spill to global memory), the compiler handles it.
However, 512 threads per block means we can only have up to 32 blocks per SM (2048 threads / 512 = 4 concurrent blocks). But with 1024 blocks and 132 SMs, we need 7.8 blocks per SM. So we'd process in 2 waves. That's OK.
Wait, with 512 threads per block, total threads = 512 * 1024 = 524288, which is exactly the number of queries. Each thread 1 query, perfect.
Then K and V are loaded once per block. Total traffic:
- Q: 2GB
- K: 2GB
- V: 2GB
- O: 2GB
- Score spill/read-back: 2GB (write + read = 1GB each way? Actually local memory read+write)
Total: 10GB
At 3.35TB/s = 3ms for memory. Plus compute. This could be very fast!
But wait - will local memory be coalesced? With 512 threads, each accessing 512 floats, if they're stored contiguously per thread, the accesses are NOT coalesced within a warp because thread 0 accesses elements 0,1,2,... and thread 1 accesses 512,513,... These are 512*4=2048 bytes apart, which might be coalesced across warps but not within a warp. Wait, actually a warp of 32 threads accessing 32 different 2048-byte-offset elements: each thread accesses a different cache line (assuming 128B lines). So 32 cache lines per warp transaction. That's terrible.
To make score storage coalesced, we need thread t to access element t, t+blockDim, t+2*blockDim, etc. So we need to store scores as transposed[score_idx][thread] instead of [thread][score_idx]. But these scores are in local memory, and the compiler decides the layout. The layout for local arrays is typically row-major with the thread index as the slowest varying dimension, which means consecutive elements in the array are consecutive in memory.
Wait, for `float scores[512];` in a thread, the compiler allocates `512 * 4 = 2048` bytes per thread. Memory layout is thread 0: bytes 0-2047, thread 1: bytes 2048-4095, etc. (Actually it's often interleaved in groups of 128B for coalescing, but for spilled arrays I'm not sure.)
If the layout is thread-interleaved, accessing `scores[i]` by all threads in a warp would be perfectly coalesced! Wait, how does the compiler store local arrays?
I recall NVIDIA docs say that access to local memory is only ever fully coalesced when ALL threads in a warp access the same address. Otherwise, it's strided. Actually, for 1D arrays, the compiler tries to ensure that `array[i]` by thread `t` is coalesced with `array[i]` by thread `t+1`. This means the elements `scores[i]` for all threads in a warp are stored together.
Actually, the standard compiler behavior for `float local[512]` is to allocate them such that `local[i]` for all threads in a warp are contiguous. More precisely, the layout is: for each `i`, `thread0_local[i], thread1_local[i], ..., thread31_local[i], padding to 128B boundary`. This makes per-warp access to `local[i]` coalesced!
But a warp of threads all accessing `local[i]` at the same time would be coalesced. In our algorithm:
Pass 1:
- Thread t computes score for key k over all D, writes to local[k].
- At key k, all threads in a warp write to their local[k]. This IS coalesced!
Pass 2:
- Thread t reads local[k] for key k.
- All threads in a warp read from their own local[k] at the same time, which is also coalesced.
WAIT, but pass 2 also needs ALL scores for a thread to compute softmax. The softmax reads scores sequentially:
```cpp
for (int k=0; k<N; k++) s[k] = scores[k];
```
With `float scores[512]`, and all threads reading `scores[0]`, then `scores[1]`, etc., at each k all threads read their local[k]. This is coalesced.
Then we need to find max over 512 scores per thread. This is done in registers, no memory issue.
Then compute exp(scores[k] - max) and sum. Again, read scores[k] sequentially, coalesced.
Then accumulate V: read scores[k] again for weights.
Wait, do we need to store scores at all? If we do two passes:
- Pass 1: compute all scores, track max.
- We only need max from pass 1. We don't need to store scores!
- Pass 2: recompute scores and accumulate V with softmax weights.
This avoids storing scores entirely! But it recomputes the dot products:
- Pass 1: Q·K for all 512 keys = 512*1024 = 524k FMAs.
- Pass 2: Q·K again for all 512 keys = 524k FMAs.
- Plus V accumulation: 512*1024 = 524k FMAs.
Total: ~1.57M FMAs per query.
Compute total: 524288 queries * 1.57M = 823B FMAs.
At 67 TFLOPS = 12.3ms.
Memory:
- Q: 2GB (loaded twice, once per pass) = 4GB
- K: 2GB (loaded twice) = 4GB
- V: 2GB (loaded once) = 2GB
- O: 2GB
Total: 12GB. At 3.35TB/s = 3.6ms.
The total could be around 12-16ms. But if we're compute bound at 12ms, that's the time.
But PyTorch reference might do something slower, especially if it's not using flash attention.
However, this approach reads Q and K twice each. Can we do better?
Alternative: keep scores in registers or shared memory to avoid recomputation.
With 512 threads per block processing 512 queries in one head:
- We could have shared memory store a tile of K and a tile of scores.
- Actually, instead of each thread storing its own 512 scores, we could process keys in tiles.
- For each key tile of 16 keys:
- All 512 threads compute their scores for these 16 keys.
- For 512 threads, 16 keys each, that's 8192 scores.
- Store in shared memory? 8192 * 4 = 32KB. That's very reasonable!
- Run softmax on the 16 scores + running stats
- Accumulate V for these 16 keys
- Discard scores
Wait, this is EXACTLY the online softmax tile approach, but with the twist that each thread handles one query and all keys in the tile, and scores for a tile are stored in shared memory or registers.
Actually, with 16 keys per tile:
- Each thread computes 16 scores (one per key in tile).
- It can keep them in 16 registers!
- Then do warp-level or thread-level softmax over 16 elements.
- With online softmax, we update running stats and O accumulator.
- No need to store scores in shared memory at all!
This is the key insight. Let me formalize:
Per thread (handles 1 query):
```cpp
float m = -1e20f; // running max
float l = 0.0f; // running sum exp
float o[DT]; // partial output accumulator, for D processed in chunks of DT
for (int kt = 0; kt < N; kt += BK) { // BK = 16
// compute BK scores for this query
float s[BK]; // 16 registers
#pragma unroll
for (int k = 0; k < BK; ++k) {
s[k] = 0;
}
for (int d = 0; d < D; ++d) {
float q_val = Q[q*D + d]; // each thread loads its own Q once per D
#pragma unroll
for (int k = 0; k < BK; ++k) {
float k_val = K[(kt+k)*D + d]; // broadcast or coalesced load from shared
s[k] += q_val * k_val;
}
}
// apply scale
#pragma unroll
for (int k = 0; k < BK; ++k) s[k] *= scale;
// online softmax for this tile
float tile_max = s[0];
#pragma unroll
for (int k = 1; k < BK; ++k) tile_max = fmaxf(tile_max, s[k]);
float new_m = fmaxf(m, tile_max);
float corr_m = expf(m - new_m);
float corr_tile = expf(tile_max - new_m);
float tile_sum = 0;
#pragma unroll
for (int k = 0; k < BK; ++k) tile_sum += expf(s[k] - tile_max);
// update O accumulator
for each d_chunk:
load O accumulators
for k in 0..BK-1:
p = expf(s[k] - new_m);
V_val = V[(kt+k)*D + d];
o[d_chunk] = o[d_chunk] * corr_m + p * V_val; // wait, this isn't right.
```
Actually, the standard online softmax update for O is:
```cpp
O_new = O_old * exp(m_old - m_new) + sum(exp(s_i - m_new) * V_i)
```
So for each key:
```cpp
float p = expf(s[k] - new_m);
#pragma unroll
for (int d = 0; d < DT; ++d) {
O[d] = O[d] * corr_m + p * V[(kt+k)*D + d];
}
l = l * corr_m + tile_sum * corr_tile;
```
But wait, `corr_m` is applied to all elements of O. This is correct.
And `p = expf(s[k] - new_m)` is the weight for each key.
The issue is: we need to iterate over keys first, then D, or D first then keys?
Standard approach: iterate over D chunks, and for each D chunk, iterate over keys to accumulate.
But with 16 keys and DT=64, the inner loop would be:
```cpp
for each D chunk:
for each key k:
p = expf(s[k] - new_m);
for d in chunk:
o[d] = o[d] * corr_m + p * V_val;
```
Wait, `corr_m` should be applied ONCE to all of O, not per key. So:
```cpp
for each D chunk:
#pragma unroll
for (int d = 0; d < DT; ++d) o[d] *= corr_m;
for each key k:
p = expf(s[k] - new_m);
#pragma unroll
for (int d = 0; d < DT; ++d) o[d] += p * V[(kt+k)*D + base + d];
```
This reduces the number of multiplications by factor BK. But we need `corr_m` per tile.
Wait, is this correct? Online softmax update:
At tile t, we have old accumulator `O_old` with old max `m_old`.
We compute tile max `m_t` and new max `m_new = max(m_old, m_t)`.
The full output should be:
`O_new = O_old * exp(m_old - m_new) / Z_new * Z_old + sum_i exp(s_i - m_new) * V_i / Z_new`
But in the standard flash attention formulation, the running O accumulator stores `O = sum(exp(s_i - m_running) * V_i)`, and we don't divide by L until the end.
So when we update:
```cpp
O_new = O_old * exp(m_old - m_new) + sum_i exp(s_i - m_new) * V_i
L_new = L_old * exp(m_old - m_new) + sum_i exp(s_i - m_new)
```
Then at end, `O_final = O_new / L_new`.
So yes, `corr_m = exp(m_old - m_new)` should be applied to the old accumulator once. Then add contributions from the new tile.
But there's a subtlety: if the old accumulator O_old corresponds to previous keys, and then we add new key contributions, each with their own `exp(s_i - m_new) * V_i`, that's correct.
So the algorithm per D chunk:
```cpp
for (int d = 0; d < DT; ++d) o_local[d] *= corr_m;
for (int k = 0; k < BK; ++k) {
float p = expf(s[k] - new_m);
for (int d = 0; d < DT; ++d) {
o_local[d] += p * smem_v[k][d];
}
}
```
With BK=16, we do 16 loads of V and 16 * DT FMAs. If DT=32, that's 512 FMAs. Plus the initial multiply by corr_m (32 FMAs). And the 16 exp() computations.
For K loading, we can put K in shared memory to avoid each thread loading K independently.
- For BK=16 and D=1024, K tile = 16*1024*4 = 64KB. This fits in shared memory!
- Similarly V tile = 64KB.
- Total shared memory = 128KB + a bit for smem_q etc. Under 227KB.
Wait, each thread needs to load Q values. Each thread processes one query. Q for that query can be loaded directly from global memory. Actually, Q is only read once per query, and each thread handles one query, so Q can be in global memory. But to maximize throughput, since each query's D=1024 is sequential, it's fine.
Actually, for the dot product, each thread loads Q[d] and multiplies by K[k][d]. K[k][d] should be in shared memory for good reuse. There are 512 threads per block, all processing queries from the same head but different queries. They all access the SAME K[k][d] for their respective dot products. So K should definitely be in shared memory to avoid 512 redundant global loads.
Let's structure a block to process one head:
- blockDim.x = 512 threads (or 256 if 512 is too many and limits occupancy)
- Each thread handles 1 query.
- 512 queries per head means exactly 1 block per head.
- Grid = (B, H) = (32, 32)
Wait, with 512 threads per block:
- On H100, max threads per SM = 2048 (if using 32 registers, but we need maybe 64-128).
- With 512 threads per block, we get at most 4 blocks per SM.
- Total blocks = 1024. 1024 / 132 = 7.76 blocks per SM.
- So occupancy = 4 blocks concurrently? Actually 7.76 total, so about 2 waves. Fine.
Register usage:
- Each thread needs: q_val (1), score array BK=16 (16), O accumulator DT (say 32), plus some temporaries.
- 16 + 32 + ~10 = ~58 registers. Well within 255 limit.
- With 512 threads * 58 regs = ~30K registers. SM has 64K, so we can fit 2 blocks = 1024 threads comfortably. Actually 64K/58 = 1100 threads per SM. So with 512 threads, we can fit 1 block = 512 threads. With 2 blocks = 1024 threads, regs = 1024*58=59K, which is OK.
Wait, I also need to consider shared memory:
- K tile: 64KB
- V tile: 64KB
- Total: 128KB
- Max shared memory per SM = 227KB
- With 1 block per SM: 128KB < 227KB, OK.
- With 2 blocks per SM: 256KB > 227KB. Oops.
So we can only have 1 block per SM concurrently because of shared memory. Total concurrent blocks = 132. But we have 1024 blocks. So it takes 1024/132 = 7.76 waves. That's about 8 waves. If each wave takes ~5ms, total = 40ms. Not great.
But wait, can't we use less shared memory? If we have 2 blocks per SM, we need shared memory < 113KB per block.
- If BK=8 instead of 16: K=32KB, V=32KB, total=64KB. Then 2 blocks = 128KB < 227KB. OK.
- Or we could use a different approach.
With BK=8:
- Number of key tiles = 512/8 = 64.
- Shared memory per block = 64KB.
- 2 blocks per SM = 128KB. OK.
- Total concurrent blocks = 264.
- Waves = 1024/264 = 3.88. Better.
But BK smaller means more tiles, more softmax iterations, more overhead.
What if we use 256 threads per block, processing 256 queries per head, with 2 blocks per head?
- Block = 256 threads.
- BK=8. K=32KB, V=32KB. Total smem=64KB.
- 2 blocks per head = 2 blocks loading K and V independently.
- 2 blocks per SM = 128KB smem used. OK.
- 2 blocks * 256 threads = 512 threads per SM concurrent.
- Total blocks = 32*32*2 = 2048.
- Waves = 2048 / (132*2) = 7.76. Same issue as 512 threads.
Hmm, the issue is the total number of head-blocks and how many fit concurrently.
What if we process multiple heads per block? No, because queries are independent.
What if we reduce the number of heads per grid launch by merging the head dimension into the block?
Actually, let me reconsider: why do we need K and V in shared memory at all?
If each of the 512 threads loads K[k][d] independently from global memory, the cache will handle it! Specifically:
- All 512 threads want K[k][d] at the same time (for different queries).
- H100 L1 cache can broadcast a single value across a warp, but 512 threads = 16 warps.
- L2 cache is shared across all SMs. If a block has 512 threads on one SM, they share the same L1. A single L1 cache line can serve at most... hmm, I don't think L1 broadcasts to multiple warps within the same SM simultaneously.
Actually, L1 cache on Ampere/Hopper can serve multiple requests per cycle but not necessarily broadcast to 512 threads. If all 512 threads request the same address, it might take multiple cycles to broadcast.
Shared memory, on the other hand, CAN broadcast to all threads in a block in a single cycle! Wait, can shared memory broadcast to all threads? Shared memory has multiple banks. For a broadcast, all threads access the same bank address, which causes bank conflicts but some hardware might handle it specially. Actually, reading the same shared memory address from all threads in a warp might be broadcast (no conflicts), but from 16 warps, you need 16 access ports or time-multiplexing.
Wait, shared memory banks: 32 banks, each 32 bits wide. To broadcast, multiple threads access the same 32-bit address. A broadcast occurs when all threads in a warp read the same address - this is handled with no bank conflicts. But if 16 warps all try to read the same shared memory address simultaneously, they can't all be served in one cycle because shared memory only has 32 banks/ports. It would take multiple cycles, but still much faster than multiple L1/L2 accesses.
So yes, shared memory is good for data that's read by many threads.
But shared memory limits occupancy. Can we use L1 caching effectively without explicit shared memory?
For 512 threads per block processing one head:
- When processing key k, all threads load K[k][d] and V[k][d].
- If K and V are accessed with a stride pattern, are they cached?
- The first warp loads K[k][d] for d=0. The next warp loads K[k][d] for d=0 as well (since all queries use the same key). All warps load the SAME address. This is not good for cache but might be handled.
Actually, if we process D as the inner loop:
```cpp
for (int d = 0; d < D; ++d) {
float k_val = K[(kt + k) * D + d];
for (int k_local = 0; k_local < BK; ++k_local) {
...
}
}
```
Wait, no. Let me think about the memory pattern.
If we don't use shared memory for K/V, and each thread does:
```cpp
for each key k:
float score = 0;
for each d:
score += Q[q*D + d] * K[k*D + d];
```
Then all 512 threads (for 512 queries) are executing:
```
tid 0: load Q[0*D+d], K[k*D+d]
tid 1: load Q[1*D+d], K[k*D+d]
...
tid 511: load Q[511*D+d], K[k*D+d]
```
For a given d and k, all 512 threads load the SAME K[k*D+d] address, but different Q addresses.
The Q loads are perfectly coalesced (thread t loads Q[t*D+d], consecutive threads load consecutive addresses if D is large... wait, no. Q is laid out as [N][D], so Q[0][d], Q[1][d] are D floats apart. That's a stride of D*4 = 4096 bytes. That's terrible for coalescing!
Wait, the layout is [B, H, N, D]. So for a fixed b and h, the data for query q starts at offset q*D. Query q+1 is at (q+1)*D. So Q[q][d] and Q[q+1][d] are separated by D elements = D*4 bytes = 4096 bytes. A warp of 32 threads accessing Q[0][d], Q[1][d], ..., Q[31][d] has a stride of 4096 bytes between consecutive threads. That's NOT coalesced. Each thread in the warp accesses a different 128B cache line (4096/128 = 32 lines apart). Each warp transaction needs to fetch 32 cache lines = 4096 bytes. That's the worst possible pattern.
This is a CRITICAL realization. For attention where each thread handles one query, we need the Q, K, V access pattern to be coalesced.
To access Q coalesced, threads in a warp should access the same query's consecutive elements, or different elements of the SAME query.
If a warp of 32 threads processes ONE query together, then they access Q[q][d+tid], which is perfectly coalesced. That's the flash attention style: one warp or block of warps processes a tile of queries, and threads within the tile access different D elements of the SAME query or queries.
In my original kernel, I had a block processing BQ=8 queries, and threads collectively loaded Q[0..7][d] with coalesced accesses. That part was actually better.
So the single-thread-per-query approach fails on Q coalescing. We need thread cooperation within queries.
Back to the original kernel style, but optimized:
- A block of threads cooperates to process BQ queries.
- Threads in a warp access the same query or different D elements.
Actually, the original kernel had:
```cpp
for (int i = tid; i < BQ * DT; i += blockDim.x) {
int q = i / DT;
int d = i % DT;
smem_q[q][d] = Q[head_offset + q_idx * D + dt + d];
}
```
This is coalesced because consecutive threads load consecutive memory addresses (outer loop over BQ*DT in thread-index-linear order, and memory is laid out as [N][D], so q_idx*D + dt + d are consecutive when d increments).
So the original kernel's memory access for Q loading was actually perfect.
The main performance issues with the original kernel were:
1. Wasted threads (only 16 of 32 lanes used for BK=16)
2. Small BQ (8 queries), meaning K and V are loaded many times.
Can we increase BQ to amortize K/V loads?
With shared memory = 227KB and D=1024:
- Q tile: BQ * DT (but DT doesn't have to be all of D)
- Actually, we load Q in chunks of DT into shared memory. So smem_q size = BQ * DT * 4 bytes.
- K tile: BK * DT * 4 bytes.
- O accumulator: BQ * D * 4 bytes.
With D=1024, O accumulator is BQ * 4096 bytes. If BQ=16, O = 64KB. If BQ=32, O = 128KB.
Q tile (for DT=64): 16*64*4 = 4KB.
K tile: BK*64*4. If BK=32, K=8KB.
P tile: BQ*BK*4 = 16*32*4 = 2KB.
M/L: BQ*8 = 128 bytes.
Total: 64 + 4 + 8 + 2 + ... = ~80KB for BQ=16, BK=32. That's fine!
So we can use BQ=16 and BK=32.
- Blocks per head = 512/16 = 32.
- Total blocks = 32*32*32 = 32768 blocks.
- K read per head = 32 blocks * (512/32) tiles * 32*1024*4 bytes? Wait.
- Each block loads K 512/BK = 16 times.
- Total K reads per head = 32 blocks * 16 tiles * 32*1024*4 bytes = 32 * 16 * 128KB = 64MB per head.
- Total for 1024 heads = 64GB.
- Plus V = 64GB. Plus Q = 2GB. O = 2GB.
- Total = 132GB. That's... a lot, but the original BQ=8 would be 128GB of K/V + 2GB + 2GB = 264GB.
Ah, I see! Doubling BQ halves the number of blocks per head, but each block still loads the full K and V. Wait, no:
- BQ=8: 64 blocks per head, each loads all K and V. Total K reads = 64 * 2MB = 128MB per head.
- BQ=16: 32 blocks per head, each loads all K and V. Total K reads = 32 * 2MB = 64MB per head.
- BQ=32: 16 blocks per head, K reads = 32MB per head.
- BQ=64: 8 blocks per head, K reads = 16MB per head.
- BQ=128: 4 blocks per head, K reads = 8MB per head.
- BQ=256: 2 blocks per head, K reads = 4MB per head.
- BQ=512: 1 block per head, K reads = 2MB per head.
So as BQ increases, K/V reads decrease. But so does parallelism (fewer blocks). With BQ=512, we have 1024 blocks. That's enough for 132 SMs (7.8 blocks per SM). Not bad.
But with BQ=512, the block has 512 threads processing 512 queries. We need a good strategy for this.
Let me design the BQ=512, BK=64 kernel. But BK=64 means K tile shared mem for one tile of 64*1024 = 256KB, which exceeds shared memory. We need to process K in smaller tiles but load them in shared memory chunks.
What if we process D in chunks? For each D chunk of DT=64:
- Load Q chunk: 512*64 = 32K floats = 128KB into shared memory. That's too much! 128KB is on the edge but actually fits.
Wait, BQ=512, DT=64, Q in smem = 512*64*4 = 128KB. K in smem = BK*64*4. If BK=32, K=8KB. O accumulator = we don't need all of O in smem if we process D chunks and update output directly.
Actually, we can do an interesting thing:
Since N=512 and D=1024, and we want to process all queries in a head in one block.
For O dimension, we can write directly to global memory without a separate accumulator tile. Or keep a partial O in registers per thread.
But each thread handling one query would have the strided Q access problem. So we need warps to handle multiple queries with cooperative loading.
Let me try a different approach altogether. What if we keep it simple AND fast:
- Use `cublasGemmStridedBatched` for both matmuls via `torch`'s operations, but my own custom fused softmax in between.
- Actually, since I can write C++ CUDA code with torch extensions, I can use cublas from within the CUDA extension.
Let me check if I can link cublas with `load_inline`. The `extra_ldflags` parameter allows passing linker flags. And I can include `cublas_v2.h`.
However, two separate cublas calls + a custom softmax kernel would involve multiple kernel launches and global memory round-trips:
1. S = Q @ K^T (writes 1GB)
2. softmax(S) (reads 1GB, writes 1GB)
3. O = S @ V (reads 1GB + 2GB, writes 2GB)
Total extra traffic: 7GB. Not bad, plus the original 4GB for Q/K/V/O = 11GB total.
But if PyTorch's SDPA uses flash attention, it avoids materializing S and would be ~8GB total. The 3GB extra from materializing S might make my approach slower.
However, if PyTorch falls back to the math backend, it does explicitly materialize everything. In that case, using cublas might be just as fast or faster (because of highly optimized GEMM).
Let me investigate: what does PyTorch `scaled_dot_product_attention` do for FP32?
According to PyTorch docs and code:
- It runs a backend selection algorithm:
- If CUDA and cuDNN available and meets conditions, use cuDNN.
- Else if flash attention available and dtype is fp16/bf16, use flash.
- Else if mem_eff attention available and appropriate, use it.
- Else use the math backend.
For FP32, cuDNN SDPA might still be available. But cuDNN SDPA in FP32 might not be as optimized as flash attention. Still, it's likely very fast.
Actually, I just realized: I don't actually know what backend the reference uses. But to beat the reference, I need to be faster than whatever it uses. If it uses flash attention, I probably can't beat it with a simple custom kernel. If it uses math backend, I might be able to beat it.
But given that this is on H100 with the latest PyTorch, I bet it uses something quite optimized (cuDNN or flash attention).
However, maybe there's a winning move: use Tensor Cores. Even for FP32, H100's Tensor Cores can process matrices very fast. But using them in custom CUDA is extremely hard.
Wait, what if I just use a much simpler approach: transpose the dimensions so the GEMMs are more efficient, and use `wmma`? No, too complex.
Let me reconsider my first kernel and think about how to optimize it for this specific hardware.
In the first kernel:
- BQ=8, BK=16
- Lane_id < BK means only 16 lanes used for score storage
- This means 50% thread waste in the critical loops
If I change BK=32, all 32 lanes are used.
- q_local = warp_id (but only warp_id 0..7 with 8 warps)
- For block=256 threads, 8 warps: yes, warp_id 0..7.
With BQ=8 and BK=32:
- Score matrix = 8 x 32 = 256 elements.
- With 8 warps, each warp handles 1 query and all 32 keys. That's perfect! Each lane handles exactly 1 key.
- Dot product parallelization: each lane computes its own dot product for its key.
- But D=1024 is still done as scalar loop. Can we parallelize the dot product?
For each key in the tile, each lane computes Q[q]·K[k]. With D=1024, a lane would do 1024 FMAs serially. For BK=32 keys, all lanes work in parallel. So each warp does 32*1024 = 32K FMAs.
But actually, the loops in my first kernel load Q and K in chunks. That's fine. The issue is: can we load K more efficiently using vector loads?
Vector loads: if we load 4 floats at a time (float4), we reduce memory instructions.
For Q loading: BQ*D = 8*1024 = 8K floats. With float4, 2K loads.
For K loading per tile: BK*D = 32*1024 = 32K floats = 8K float4 loads.
For a block of 256 threads, 8K/256 = 32 loads per thread per K tile. Not bad.
For O accumulator: BQ*D = 8*1024 = 8K floats. Writing these out is 2K float4 stores.
Let me check: what if I increase BQ to 16 and BK to 32? Shared memory:
- smem_q: 16*64*4 = 4096 bytes (DT=64)
- smem_kv: 32*64*4 = 8192 bytes
- smem_p: 16*32*4 = 2048 bytes
- smem_m/l: 16*8 = 128 bytes
- smem_o: 16*1024*4 = 65536 bytes
Total = ~80KB. Under 227KB limit. Great!
Can we have 256 threads (8 warps) handle BQ=16 queries?
- Each warp would handle 2 queries.
- warp 0: queries 0,1
- warp 1: queries 2,3
- etc.
- But the score matrix is 16x32 = 512 elements. 512 elements / 256 threads = 2 elements per thread.
- Each thread would compute 2 scores. Each score is a dot product of D=1024.
For dot product chunk: load Q[0..15][0..63] and K[0..31][0..63].
- Q: 16*64 = 1024 floats. 256 threads load 4 floats each (1 float4). Coalesced!
- K: 32*64 = 2048 floats. 256 threads load 8 floats each (2 float4). Coalesced!
Then for the dot products:
- Each thread handles... hmm, with 256 threads and 16*32=512 scores, each thread does 2 dot products.
- Thread i handles score for query (i/32) and key (i%32)? But i/32 gives warp_id (0..7), which only gets 8 queries. We'd need to map carefully.
Actually, 256 threads = 8 warps * 32 lanes. If BQ=16 and BK=32:
- We could have each warp handle 2 queries * 32 keys = 64 scores.
- 8 warps * 64 = 512 scores. Perfect.
- Warp w handles queries 2w and 2w+1.
- In warp w, lane k handles key k for both queries (or just one query?).
Wait, there are 32 lanes and 32 keys. So each lane can handle exactly 1 key.
But we need to handle 2 queries per warp. So each lane needs to handle 2 dot products: one for query_a, one for query_b.
Dot product:
```cpp
float acc_a = 0, acc_b = 0;
for (int d = 0; d < DT; ++d) {
float q_a = smem_q[2*w][d];
float q_b = smem_q[2*w+1][d];
float k = smem_kv[lane_id][d];
acc_a += q_a * k;
acc_b += q_b * k;
}
```
Each lane does 2 FMAs per iteration. That's actually good for instruction throughput!
After all DT chunks, each lane has acc_a and acc_b for key k.
Then lane k writes acc_a and acc_b to smem_p. But smem_p needs to store for 16 queries x 32 keys.
Wait, instead of storing in shared memory, we could keep scores in registers and do the softmax directly. But then we need all scores for one query to compute its max. Within a warp, for query_a, lanes 0..31 have acc_a for keys 0..31. We can do warp shuffles to find the max and sum for query_a and query_b.
Yes! This is very elegant:
- Each warp handles 2 queries.
- All 32 lanes participate.
- After computing dot products, lane k has score for key k.
- We do warp-level reduction for max and sum for both queries simultaneously.
- Each lane computes its own exp and weight.
- Then for V accumulation, each lane can help compute partial contributions.
For V accumulation with BK=32 keys and D=1024:
- We need O for 2 queries, each O is 1024 floats.
- With 32 lanes, we can have lane d handle dimension d for both queries.
- Wait, D=1024 and 32 lanes means each lane handles 32 dimensions.
- For each key k, we need V[k][:] for both queries.
- Lane d handles dim d, d+32, d+64, ..., d+992 (32 dimensions total).
- For query_a and query_b, each lane computes:
`O_a[dim_j] += p_a[k] * V[k][dim_j]` for j=0..31
This requires loading V[k][d] for the lane's 32 dims. But V is stored as [N][D], so V[k][d] and V[k][d+1] are contiguous. Good.
We can either load V directly from global memory or use shared memory.
If we use shared memory for V tile:
- For a BK=32 key tile and DT=64 D chunk: smem_v[32][64] = 8KB.
- But if each lane handles 32 scattered dims, we can't easily map to a fixed DT tile.
Maybe we process V in chunks of 32 dims:
- For d_base = 0, 32, 64, ..., 992:
- Load V[k][d_base + d] for k=0..31, d=0..31 into registers or shared memory.
- Each lane owns dim d_base + lane_id.
- For both queries, accumulate p_a * V and p_b * V.
- This assumes p is the same for both queries across keys.
Wait, we don't need shared memory for V if we process D=1024 in chunks of 32 and let each lane load its own V values.
Let's formalize this approach:
Block configuration:
- BQ=16, BK=32
- blockDim = 256 threads (8 warps)
- Each warp handles 2 queries.
Algorithm per warp (for query pair q0, q1):
For each key tile kt = 0..N-1 step BK:
// Load K tile for BK=32 keys, D=1024 dims
// We process D in chunks of DT=64 or 128
float s0 = 0, s1 = 0; // scores for q0, q1
for (int dt = 0; dt < D; dt += DT) {
// Load Q for q0, q1: 2*DT floats
// This could be done cooperatively by all warps once per block.
// Actually, all 8 warps need Q for all 16 queries.
// Load K for this warp's 32 keys: 32*DT floats.
// Warp loads its own K tile into registers or shared memory.
// For each d in dt..dt+DT-1:
// k_val = K[kt + lane_id][d] (if K in registers, each lane loads its key)
// q0_val = Q[q0][d]
// q1_val = Q[q1][d]
// s0 += q0_val * k_val;
// s1 += q1_val * k_val;
}
// s0 and s1 are now dot products for key = lane_id.
// scale
s0 *= scale; s1 *= scale;
// warp-level softmax for q0:
// max0 = warp_max(s0)
// sum0 = warp_sum(exp(s0 - max0))
// warp-level softmax for q1:
// max1 = warp_max(s1)
// sum1 = warp_sum(exp(s1 - max1))
// update running m0, l0, m1, l1
// compute p0 = exp(s0 - new_m0), p1 = exp(s1 - new_m1)
// For V accumulation:
// Each lane also handles some D dimensions.
// wt = lane_id % something?
// Actually for V accumulation, all 32 lanes can help for BOTH queries.
// We have 32 lanes and D=1024 dims. Each lane can handle 32 dims.
// Lane d handles dim d, d+32, d+64, ..., d+992.
for (int d_chunk = 0; d_chunk < D; d_chunk += 32) {
float v_val = V[kt + lane_id][d_chunk + something]? No, we need all keys.
// For key k, lane j handles dim d_chunk + j.
// But each key has its own V vector.
// Actually, for a fixed D chunk, we need V[0..31][dim] for all keys.
// If we process one key at a time:
for (int k = 0; k < BK; ++k) {
float p0_k = expf(s0_for_k - new_m0); // But only the lane that computed key k has s0!
// We need to broadcast p0_k to all lanes.
float p0 = __shfl_sync(0xFFFFFFFF, p0_for_some_lane, k);
// and similarly p1.
float v_val = V[kt + k][d_chunk + lane_id]; // coalesced if V is in shared or global with proper layout
O0_acc[d_chunk/32] += p0 * v_val;
O1_acc[d_chunk/32] += p1 * v_val;
}
}
Wait, there's an issue: `s0` is in each lane, but lane k has s0 for key k. To get s0 for key k, we broadcast from lane k: `float s0_k = __shfl_sync(0xFFFFFFFF, s0, k);`. Then we can compute p0_k.
This works! And V[kt+k][d_chunk + lane_id] should be loaded coalesced if we load from global memory with consecutive lanes accessing consecutive dims.
But we also need O0_acc to have 32 registers (for 32 dims). With O1_acc also 32 registers, that's 64 registers just for accumulators. Plus scores (2), plus various temporaries. That's maybe 80 registers per thread. With 256 threads per block = 20K registers. Under 64K per SM, so we can fit 3 blocks. Not bad.
For the V load: `V[kt + k][d_chunk + lane_id]`, with lane_id 0..31, accesses contiguous memory. Perfectly coalesced! But it's a global memory load. The L1 cache might handle this, but it's still going to HBM for V.
Total V traffic: each block reads all V once, but for D=1024 and 16 queries, that's the same as any fused approach.
Let me check: with BQ=16, blocks per head = 32. Each block reads V once = 2MB per head. Total V traffic = 64MB per head * 1024 heads = 64GB. Hmm, that's still high.
If we increase BQ to 64:
- 8 blocks per head.
- V traffic = 16MB per head.
- Total V = 16GB.
Can we do BQ=64 with a 256-thread block?
- 256 threads / 64 queries = 4 threads per query. That's not much.
- We want each warp to handle some queries. 256 threads = 8 warps.
- 64 queries / 8 warps = 8 queries per warp. But a warp has 32 threads, so 4 threads per query.
- 4 threads can compute 4 keys each for one query... or all 4 can cooperate on the same key.
Actually, with 4 threads per query and BK=32, we can have each query handled by 4 threads:
- Thread 0: keys 0..7
- Thread 1: keys 8..15
- Thread 2: keys 16..23
- Thread 3: keys 24..31
Then for dot product, all 4 threads can load Q cooperatively.
But warp-level shuffle reductions become more complex since the 32 lanes are divided into 8 queries × 4 threads.
Hmm, this is getting complicated. Let me think about whether we can do BQ=64 with 2 warps per query:
- 8 warps / 64 queries = 0.125 warps per query. No, that's not even integer.
OK, the issue is that D=1024 is fixed, and we need to amortize K/V loads with larger BQ, but that reduces the threads per query.
What if we use BQ=32 with 8 warps? That's 0.25 warp per query... wait, 32 queries / 8 warps = 4 queries per warp. That's what we had for BQ=16: 2 queries per warp.
With BQ=32 and 8 warps: 4 queries per warp.
Each warp handles 4 queries, and 32 lanes / 4 = 8 keys per query per warp? No, that's not right.
Actually, 4 queries in a warp, each query needs scores for BK=32 keys. That's 128 scores. But a warp only has 32 lanes. So each lane would need to handle 4 scores.
Lane assignments for 4 queries per warp, BK=32 keys:
- Query 0: lanes 0..7 handle keys 0..7, 8..15, 16..23, 24..31? No, if each lane handles 1 key, that's only 8 keys per query.
- Query 0: lanes 0..7 handle keys 0..31, 4 keys each (0,4,8,...,28 by lane 0)?
Maybe we should assign lanes like:
- Lane j handles key j for all 4 queries.
- Query i (in this warp) needs score for key j, computed by lane j.
- So lane j has scores for query 0,1,2,3 all for key j.
- `s[i] = Q[qi]·K[kj]` where i=0..3.
This means each lane computes 4 dot products, one per query. That's 4 FMAs per d iteration:
```cpp
acc0 += q0 * k;
acc1 += q1 * k;
acc2 += q2 * k;
acc3 += q3 * k;
```
With 32 lanes processing 32 keys, and 4 queries per warp, we cover 4 * 32 = 128 scores.
For BQ=32, we need 8 warps. 8 * 128 = 1024 scores total across all warps. But BQ=32 and BK=32 means only 32*32 = 1024 scores total. Perfect!
So the score matrix exactly fits: 1024 scores = 8 warps * 128 scores/warp.
Then for warp-level softmax:
- For query 0, we need all 32 scores. These are in lanes 0..31, but each lane has score for query 0 with DIFFERENT keys.
- We need `warp_max` across lanes for query 0: extract all query-0 scores via shuffle.
- Similarly for queries 1,2,3.
- But a warp only has 32 values total for query 0 (1 per lane). That's a 32-element max reduction. Easy with shuffles.
- Same for sum.
For V accumulation:
- Query 0 needs O of size D=1024.
- With 32 lanes and 4 queries in the warp, we need to divide D among the lanes.
- Total dimensions to compute = 4 * 1024 = 4096.
- 32 lanes can handle 4096 dims by each lane handling 128 dims.
- Lane j handles dim j, j+32, j+64, ..., j+992 for one query?
Actually, we need each lane to load V[k][dim] for all keys k to contribute to the output. The weight for key k and query i is only known to one lane (lane k). We need to broadcast weights.
For query 0, lane k has weight p0k. For dim d_j handled by lane j, we need:
`O0[d_j] += p0k * V[k][d_j]` for all k.
We can have each lane broadcast its p weights one by one:
```cpp
for (int k = 0; k < BK; ++k) {
float p0 = __shfl_sync(0xFFFFFFFF, p0_lane, k);
float p1 = __shfl_sync(0xFFFFFFFF, p1_lane, k);
float p2 = __shfl_sync(0xFFFFFFFF, p2_lane, k);
float p3 = __shfl_sync(0xFFFFFFFF, p3_lane, k);
// lane j loads V[k][d_start + j] for its d_start
}
```
But if each lane handles 128 dims, for a single key, each lane loads 128 elements of V. That's a lot of loads.
Wait, can we process V accumulation in chunks too? Like DT=32 per iteration?
If each lane handles dim `d = lane_id` in chunk `d_base`, then for all 32 lanes, they load V[k][d_base + 0..31] which is coalesced. Each lane gets one float.
But with 4 queries, after loading the V chunk, each lane needs to update 4 output accumulators. And it needs to receive 4 p weights (p0k, p1k, p2k, p3k) for each key k.
For BK=32 keys and DT=32 dims:
- 32 iterations of k.
- In each k: 4 broadcasts + load V[k][d_base+lane_id] (coalesced).
- 4 FMAs per lane (one per query).
- Then move to next k.
This is quite efficient on memory (32 floats loaded per k, all coalesced) and compute.
For one V tile (BK=32 keys), each lane processes 32 dims (DT=32). There are 1024/32 = 32 D tiles.
So per key tile: 32 * 32 = 1024 iterations of k+d. Actually 32 keys * 32 D tiles = 1024 operations per lane.
With 32 lanes * 4 queries = 128 accumulators total. Each lane keeps 4 * 32 (DT registers)? No, we process one D tile at a time and reuse accumulators across D tiles. But we need to keep the full O across all D tiles. Since D=1024 and DT=32, we process 32 D tiles and for each tile we have 4 query accumulators. Total registers = 4 * 32 = 128 floats.
Wait, NO. For online softmax, the O accumulator needs to be corrected (`* corr_m`) for each new key tile. The correction is per query, not per dimension. So:
```cpp
for each key tile:
compute new_m, corr_m
for each D tile:
for each query i in warp:
O_i[d] *= corr_m;
for k in 0..BK-1:
p_i_k = ...
O_i[d] += p_i_k * V[k][d];
```
For 4 queries and DT=32, O accumulators = 4 * 32 = 128 registers. That's a lot. Actually, we'd process one D tile at a time and move to the next, accumulating O for that tile. But we need O from previous key tiles to persist across D tile processing? No, we can process a D tile for ALL key tiles before moving to next D tile:
```cpp
for d_base in 0..D-1 step DT:
initialize O0[d_base..d_base+DT-1] = 0
m = -inf, l = 0
for each key tile:
compute scores (using full D, not just d_base)
update m, l
for k in 0..BK-1:
for d in 0..DT-1:
O0[d] = O0[d] * corr_m + p0k * V[k][d_base+d];
```
Wait, this works! Because the score computation always uses the FULL D dimension. We just process output dimensions in chunks.
So the algorithm is:
1. For each key tile, compute scores (dot products using all D=1024).
2. Do online softmax update.
3. For each D chunk (32 dims), update O accumulators using weights from this key tile.
4. Move to next key tile.
Register usage:
- Scores: 4 floats (one per query in warp).
- O accumulation: 4 queries * 32 dims = 128 floats.
- Q loads: 4 * 32 (DT for Q chunk, but Q for score uses full D, so Q is loaded in chunks and accumulated).
Actually, for full D dot products, we don't need all Q in registers simultaneously. We load Q in chunks:
```cpp
float s0 = 0, s1 = 0, s2 = 0, s3 = 0;
for (int dt = 0; dt < D; dt += DT) {
load Q[q0][dt]..Q[q0][dt+DT-1] into q0_vals[DT]?
...
load K[k][dt]..K[k][dt+DT-1] into k_val[DT]?
}
```
Wait, if each lane handles 4 queries, and each query's score is sum over D of Q*K, then:
- We load Q for the 4 queries into registers.
- We load K for one key into registers.
- Multiply and accumulate.
If DT=32:
- q0[0..31], q1[0..31], q2[0..31], q3[0..31] = 128 floats.
- k[0..31] = 32 floats.
- Total: 160 floats. That's a lot of registers per thread.
If DT=64: 256 + 64 = 320 floats. Too much.
So DT needs to be smaller. What if DT=16?
- 4*16 = 64 + 16 = 80 floats for Q+K.
- Plus scores (4), plus O accumulators (4 * 32 = 128 if we process D output chunks of 32).
- Total: 212 floats. That's a lot but might still be OK (255 limit).
Wait, O accumulators (128 floats) are only needed during V accumulation, not during score computation. So the compiler might reuse registers.
Let me try to write this kernel more concretely:
```cpp
#define BQ 32
#define BK 32
#define WQ 4 // queries per warp
#define DT_Q 64
#define DT_O 32
__global__ void attention_kernel(...) {
int tid = threadIdx.x;
int warp_id = tid / 32;
int lane_id = tid % 32;
int b = blockIdx.x;
int h = blockIdx.y;
// this block processes BQ=32 queries
int qg = 0; // Actually with 1 block per head and BQ=32, we need 16 blocks per head? No, grid=(B,H,16) for 512/32=16 query groups per head.
int qg = blockIdx.z; // 0..15 for N=512/BQ=32
int qid = warp_id * WQ; // 0, 4, 8, ..., 28 for warp 0..7
// queries handled by this warp: qg*BQ + qid + i, for i=0..WQ-1
int q0 = ...;
int q1 = q0 + 1;
int q2 = q0 + 2;
int q3 = q0 + 3;
float m0=-inf, m1=-inf, m2=-inf, m3=-inf;
float l0=0, l1=0, l2=0, l3=0;
// O accumulators for DT_O dims
float o0[DT_O], o1[DT_O], o2[DT_O], o3[DT_O];
// initialize to 0
for (int kt = 0; kt < N; kt += BK) {
// compute scores for this key tile
float s0 = 0, s1 = 0, s2 = 0, s3 = 0;
for (int dt = 0; dt < D; dt += DT_Q) {
// load Q for 4 queries and K for 1 key (lane_id's key)
float q0v[DT_Q/WQ]; // wait, we need all DT_Q dims.
}
}
}
```
This is getting complex, and I need to be careful about shared memory, coalescing, and register spilling.
Let me reconsider and go for a simpler but well-optimized kernel.
Given the complexity and time constraints, maybe I should try a kernel with:
- BQ=16, BK=32, DT=64
- 8 warps in a block
- Each warp handles 2 queries (16/8)
- All 32 lanes are active
- Process D=1024 in DT=64 chunks
- Shared memory: Q[16][64]=4096B, K/V[32][64]=8192B, P[16][32]=2048B
- O accumulator stays in shared memory: 16*1024*4=65536B
- Total smem = ~80KB
Wait, can each warp handle 2 queries with 32 lanes? For scores:
- Lane j computes score for key j for BOTH queries.
- So dot products: `acc0 += q0 * k_j`, `acc1 += q1 * k_j`.
- Two FMAs per iteration.
For V accumulation:
- Lane j loads V[k][d_chunk + j] for some d_chunk?
- But we need to update O for both queries.
- Each lane can handle one dimension. For D=1024, lanes 0..31 handle dim 0..31 for the first chunk, then 32..63, etc.
- For each key k, lane j gets p0k and p1k via shuffle.
- Then lane j updates o0[j] += p0k * v[k][j], and o1[j] += p1k * v[k][j].
- Wait, o0 needs 1024 dims. We can keep just the current chunk in registers, but we need to persist across key tiles.
With online softmax, O must persist across key tiles. So we need the full O[16][1024] somewhere. Options:
1. Shared memory (80KB for O alone, as noted above).
2. Global memory (slow).
3. Split D processing: process all key tiles for a D chunk, then write that chunk and move to next. But that means processing scores (which need all D) for every D chunk, which is terrible.
Shared memory for O seems unavoidable if we want to avoid global round-trips. With 80KB total smem, we can only fit 1-2 blocks per SM (SM has 227KB). That's acceptable. 1024 blocks / (132 SMs * 1 block) = 7.76 waves. Each wave processes all SMs for one iteration.
But wait, we can further reduce O size by processing fewer queries per block or using smaller accumulators. Actually, 80KB is not that large.
Let's calculate total time with this kernel:
- Memory traffic ≈ 132GB (from earlier calculation? Actually with BQ=16: K=64GB, V=64GB, Q=2GB, O=2GB = 132GB)
- 132GB / 3.35TB/s ≈ 39ms
- Plus compute.
This seems too slow. But I think my traffic calculation was wrong because with 1024 blocks and only 132 SMs and 80KB smem, the actual number of concurrent blocks per SM is limited by shared memory to maybe 2 (160KB), but with warps, it's limited by registers too.
Actually, at 132GB of traffic, the kernel is clearly bandwidth bound. Is the reference doing significantly less traffic?
If the reference uses flash attention with the optimal algorithm, it also loads Q once, K once, V once per head (if it can fit everything in tiles effectively). But with D=1024, a flash attention tile for Q might be small.
Wait, in the standard flash attention algorithm, each block loads Q tiles of size BQ x DT into SRAM, and iterates over KV tiles. The total K/V reads per head is (N/BQ) * D * N / BK? No, each query block loads all K and V. So total K reads = (N/BQ) * (N * D). Similarly for V.
With BQ=8: K reads = 64 * 512 * 1024 = 33.5M floats = 134MB per head? Wait no:
N/BQ = 512/8 = 64 query blocks per head.
Each query block loads K in tiles of BK. Total K loaded per query block = N * D = 512 * 1024 floats.
Total K reads per head = 64 * 512 * 1024 floats = 33.5M floats = 134MB.
Total for all heads = 134MB * 1024 = 137GB.
Similarly V = 137GB.
Q = 2GB.
O = 2GB.
Total = 278GB.
Ah! I was calculating 2MB per head earlier, but I forgot that each query block loads the FULL K of 2MB. Yes, 64 query blocks * 2MB = 128MB per head. That's correct.
So with my first kernel (BQ=8, BK=16), total traffic is roughly:
- K: 128MB/head * 1024 = 131GB
- V: 131GB
- Q: 2GB
- O: 2GB
Total: ~266GB
At 3.35TB/s, that's 79ms. My kernel could never beat a flash attention kernel that does only 8GB.
Wait, but does PyTorch SDPA use the optimal flash attention? For FP32 on H100, I think it might not even use flash attention (since flash attention typically only supports FP16/BF16).
If it uses cuDNN attention, cuDNN might also use a similar tiling algorithm. But cuDNN is highly optimized and might use larger tiles or Tensor Cores.
If the reference uses the `math` backend, it explicitly computes: `softmax(Q @ K^T / scale) @ V`
This involves:
- Q @ K^T: [1024, 512, 1024] batch GEMM. cuBLAS handles this very efficiently and might even use Tensor Cores.
- Writes S = [1024, 512, 512] = 1GB
- Softmax on S: 1GB read + 1GB write
- S @ V: [1024, 512, 512] @ [1024, 512, 1024] = 1GB read + 2GB read + 2GB write
Total: Q(2GB) + K(2GB) + S(2GB) + V(2GB) + O(2GB) = 10GB
10GB / 3.35TB/s = 3ms for memory. Plus compute for GEMMs. cuBLAS GEMM at these sizes... [1024, 512, 1024] is actually [1024 batches, M=512, K=1024, N=512]. This is a fairly small GEMM. cuBLAS might not achieve peak performance for such sizes.
The compute for Q @ K^T is B*H*N*N*D = 32*32*512*512*1024 = 274B FMAs.
274B / 67 TFLOPS = 4ms.
Plus the second GEMM: same 274B FMAs = 4ms.
Plus softmax: negligible.
Total reference = ~11ms if compute bound, or ~10ms if memory bound.
If I write a custom fused kernel with traffic ~266GB, it would take ~80ms. That would be MUCH SLOWER than the reference.
So the reference is very likely NOT using the math backend. It's probably using cuDNN attention or some other optimized path that avoids materializing S.
Hmm, but if the reference uses the math backend, my custom kernel with 266GB traffic would still lose.
Wait, maybe I should write a better custom kernel that materializes S and uses cuBLAS for GEMMs.
Actually, hold on. Can I just use `torch.matmul` inside `ModelNew.forward`? The challenge says to use a custom CUDA kernel. But I can write a kernel that uses `cublas` for the GEMMs and a custom kernel for softmax. That might be fast.
Or even simpler: maybe I can use the same approach as the math backend but optimize the softmax.
Wait, let's reconsider. Maybe PyTorch FP32 SDPA actually dispatches to something unoptimized on this platform. I need to check what the reference timing would likely be.
Actually, I can look at the problem from a different angle: the evaluator says "run faster than the PyTorch reference". If the reference takes X ms, I need < X ms. The main variable I control is the algorithm.
Given that custom kernels for attention without Tensor Cores are usually slower than FlashAttention/cuDNN, but if the reference doesn't use those...
Actually, wait. For FP32, SDPA might fall back to `math` because:
- FlashAttention doesn't support FP32
- Memory-efficient attention might not support FP32 either (or it's disabled for some shapes)
- cuDNN attention requires certain conditions
If the reference falls back to `math`, it computes three separate operations: matmul, softmax, matmul. My custom fused kernel would likely be faster if I can make it reasonably efficient, even with 266GB traffic? No, 266GB is way too much.
Unless... wait. My traffic calculation was for the worst case where K/V aren't cached. But with the custom kernel running on 132 SMs, and 1024 blocks, the K might stay in L2 if multiple blocks from the same head run on the same SM.
Actually, let's recalculate more carefully. A head has 64 blocks (for BQ=8). Each block reads 2MB of K. 64 blocks need to run. H100 has 132 SMs. So it's entirely possible that, say, 10 blocks from the same head run on 10 different SMs, none sharing cache. Then each reads K from HBM. But maybe the scheduler puts blocks from the same head together?
Actually, CUDA block scheduling for 2D grid (B, H) usually assigns blocks in order to SMs. If the grid is (32, 32) and each block processes multiple query groups... in my first kernel there was no z dimension, so grid=(32,32), each block processed 1 query group. The blocks would be scheduled as (0,0), (1,0), ..., (31,0), (0,1), etc. So consecutive blocks in x are from the same head. This means blocks from the same head will indeed run on the same SMs, potentially sharing L2 cache.
But L2 cache is only 50MB. K for one head is 2MB. We can fit 25 heads in L2. With 32 blocks per head all wanting K, and they might not all fit in L2. But K is 2MB, and if we read it 64 times (for 64 query groups), L2 can serve most of those reads if the working set is 2MB. Actually, L2 has much higher bandwidth than HBM (maybe 10-20 TB/s). So if K stays in L2 after the first read, subsequent reads come from L2. Then total K traffic from HBM is only 2MB per head, not 128MB!
Is this realistic? For a single head, all 64 blocks access the same K. The first block reads K from HBM into L2. The next 63 blocks read from L2. So K HBM traffic is indeed 2MB per head, not 128MB!
Similarly for V. So my total traffic calculation was overly pessimistic.
But what about Q? Each block reads a different Q tile. Those don't share across blocks.
So total traffic with L2 caching of K/V:
- Q: 2GB
- K: 2GB (from HBM to L2, then served from L2)
- V: 2GB
- O: 2GB
Total: 8GB
At 3.35TB/s, that's 2.4ms. Plus some L2 cache overhead. This is excellent!
But wait, is L2 cache big enough? Per head, K+V = 4MB. With 132 SMs and blocks racing through, if 132 blocks are from different heads, L2 needs 132*4MB = 528MB of L2 cache, but we only have 50MB. So most K/V will be evicted.
However, let's think about the wave structure. 1024 blocks total. First wave of ~132 blocks. If these are from 132 different x-coordinates (batches) and the same y-coordinate (head), then L2 needs only 4MB (one head). But if they're from different heads, L2 needs 132*4MB.
Given grid=(B,H)=(32,32), blocks are launched in x-major order. Block sequence:
(0,0), (1,0), ..., (31,0), (0,1), (1,1), ..., (31,31)
So the first 32 blocks are all head=0, batch=0..31. They all share the same head! They all need K/V for head 0. So L2 only needs 4MB for the first 32 blocks.
Then blocks (0,1) to (31,1) need head 1. K for head 0 might be evicted.
But within the first 32 blocks, all use head 0's K/V. If they all run concurrently across 32 SMs (out of 132), then L2 bandwidth is shared. The remaining 100 SMs work on other blocks from other heads.
Actually, the block scheduler usually assigns blocks to SMs in a round-robin or similar. The exact order isn't strictly x-major across the chip, but within a cluster or GPC it might be.
On H100, with 132 SMs, and 1024 blocks, a wave of 132 blocks would include roughly 4 blocks from head 0 (batches 0,1,2,3), 4 from head 1, etc. Actually no, if scheduler takes blocks in sequence: (0,0), (1,0), ..., (131,0) but batch only goes up to 31. So first 32 blocks: (0..31, 0). Next 32: (0..31, 1). Next 32: (0..31, 2). Etc.
For a wave of 132 blocks, we'd get:
- 32 blocks from head 0 (all batches)
- 32 blocks from head 1
- 32 blocks from head 2
- 32 blocks from head 3
- 4 blocks from head 4
So 5 different heads are active. K+V per head = 4MB. 5 heads = 20MB. L2 is 50MB. So K/V for these heads might entirely fit in L2!
As blocks from head 0 complete, new blocks from head 4, 5, etc. join the active set. If the scheduler maintains a pool of ~132 blocks and replaces completed ones, we might have ~10-20 heads active at a time. 20 heads * 4MB = 80MB > 50MB L2. So some eviction occurs.
But still, K/V traffic from HBM is much less than the worst case. Maybe 10-20GB total for K+V instead of 256GB.
This makes fused attention much more viable. The original kernel might actually be decent if L2 caching works well.
Now, for performance optimization, I should:
1. Increase BQ to reduce the number of query blocks per head (reducing K/V reads, even from L2, and improving L2 hit rates).
2. Use all threads efficiently.
3. Use shared memory for K/V tiles to amortize the cost of loading from L2/L1.
With BQ=64:
- 8 blocks per head.
- A wave of 132 blocks would include ~4 heads fully + some partial.
- L2 pressure = ~4 heads * 4MB = 16MB. Excellent caching.
- Per head, K loaded 8 times from L2, but only once from HBM (if cached).
- Total K HBM traffic = 2GB. V = 2GB.
So BQ=64 seems sweet. But can we fit the kernel with BQ=64?
If each block processes 64 queries:
- With 512 threads per block (16 warps), each warp handles 4 queries.
- Warp w handles queries 4w, 4w+1, 4w+2, 4w+3.
- Each warp has 32 lanes.
- Lane j handles key j for all 4 queries? No, BK would be <=32.
- If BK=32, lane j handles key j for all 4 queries.
- Score per lane: 4 dot products of D=1024 each.
This is similar to my earlier calculation. For V accumulation:
- Lane j handles D dimensions: j, j+32, j+64, ... for all 4 queries.
- With D=1024 and 32 lanes, each lane handles 32 dims.
- For BK=32 keys and all D, each lane does 32*32 = 1024 FMAs per query? No, per 4 queries.
Actually, let's think about the total compute intensity.
Total FLOPs = 2 * B*H*N*N*D = 2 * 32*32*512*512*1024 = 549B FMAs.
At 67 TFLOPS = 8.2ms.
If memory traffic is ~8GB (Q+K+V+O), that's 2.4ms at 3.35TB/s.
So the kernel is mostly compute bound. The 8ms compute suggests we need high FMA throughput.
With a custom kernel using scalar FMAs on CUDA cores, can we achieve 67 TFLOPS?
- 67 TFLOPS requires nearly 100% utilization of FP32 units.
- A warp can do 1 FMA per thread per clock. At 2GHz, a warp does 64 GFLOPS. An SM with 4 warps (128 threads) does 256 GFLOPS. 132 SMs = 33.8 TFLOPS. Wait, that's half of 67 TFLOPS.
Actually, modern GPUs have TWO FP32 units per SM per scheduler? Let me think.
H100 has 128 FP32 CUDA cores per SM. At 2GHz, that's 256 GFLOPS per SM (128 * 2). 132 SMs = 33.8 TFLOPS. But the spec said 67 TFLOPS. How?
Ah, H100 has FP32 Tensor Cores too? No, the 67 TFLOPS is peak FP32 throughput for the whole chip. That implies some SM-level parallelism I'm not accounting for. Actually, H100 has 128 FP32 cores per SM, each can do 1 FMA per cycle. 1 FMA = 2 FLOPs. So per SM: 128 * 2 * 2GHz = 512 GFLOPS. 132 * 512 GFLOPS = 67.5 TFLOPS. Yes!
So an SM with 128 active threads (4 warps) can achieve peak if all threads issue FMAs every cycle.
My kernel needs to issue enough FMAs to hide memory latency.
With BQ=64 and BK=32, and each block having 512 threads (16 warps):
Each warp processes 4 queries, 32 keys.
Score computation: 4 queries * 32 keys * 1024 D = 131k FMAs per warp.
V accumulation: 4 queries * 32 keys * 1024 D = 131k FMAs per warp.
Total: 262k FMAs per warp.
At 512 GFLOPS per SM, 16 warps on an SM = 16*32=512 threads. If all issue FMAs, that's peak utilization.
262k FMAs / (512 threads * 1 FMA/cycle/thread) = 262k cycles. At 2GHz = 0.13ms per block.
1024 blocks / 132 SMs = 7.76 waves. 7.76 * 0.13ms = 1.0ms.
Hmm, that seems too optimistic. But memory will slow it down.
Actually, I realize I'm overcomplicating this. Let me just try to write a solid BQ=64 kernel with BK=32, 16 warps (512 threads), and see if it compiles and works.
But wait, I can't easily test performance without submitting. Let me try a more conservative approach: write a kernel that's known to be well-structured, use a moderate BQ like 32 or 64, and submit.
Actually, I have another idea. Since the problem size N=512 is relatively small, and D=1024 is large, what if I use a different parallelism strategy:
Each block = 1 warp = 32 threads.
Grid = (num_heads * num_batches * num_queries, 1) = 32 * 32 * 512 = 524288 blocks.
But CUDA has a limit on grid dimensions: 2^31-1 for x, so 524288 is fine.
Each block/thread handles 1 query.
Algorithm per thread:
```cpp
float q_val = Q[qid * D + threadIdx.x]; // Load only 1 element of Q? No.
```
Wait, no. A single thread can't efficiently compute the full dot product of D=1024, but we can use warp-level reductions.
Per warp (32 threads), for 1 query:
1. Load all of Q: 1024 floats. With 32 threads, each thread loads 32 floats. Coalesced (thread t loads Q[t*32:(t+1)*32]).
2. Iterate over all keys (512 keys):
a. Load K[k][t*32:(t+1)*32] for thread t. Coalesced.
b. Compute partial dot product per thread: thread t computes dot(Q_chunk, K_chunk).
c. Warp shuffle reduce to get score.
d. Track max for softmax.
3. Iterate over keys again for V accumulation:
a. Load V[k][t*32:(t+1)*32].
b. Broadcast weights via shuffle.
c. Accumulate to O.
This is actually a very clean approach! Let me formalize:
Each warp handles 1 query. 1 block = 32 warps = 32 queries.
Grid = (32 heads * 32 batches * 16 query_groups) = 16384 blocks.
Wait, 512 queries / 32 queries per block = 16 blocks per head.
Total blocks = 32 * 32 * 16 = 16384.
Algorithm per warp (query q):
```cpp
// Load Q into shared or local memory
float q_local[32]; // 32 floats per thread
for (int i = 0; i < 32; ++i) {
q_local[i] = Q[q * D + threadIdx.x * 32 + i]; // thread t loads elements [t*32, (t+1)*32)
}
// Wait, threadIdx.x is 0..31 for warp 0? But a block has 32 warps.
// threadIdx.x % 32 = lane_id
// We can load: q_local[i] = Q[q*D + lane_id*32 + i];
```
For score computation:
For each key k = 0..511:
float dot = 0;
for (int i = 0; i < 32; ++i) {
float k_val = K[k*D + lane_id*32 + i]; // same pattern as Q
dot += q_local[i] * k_val;
}
// warp sum
#pragma unroll
for (int offset = 16; offset > 0; offset /= 2) dot += __shfl_xor_sync(0xFFFFFFFF, dot, offset);
dot *= scale;
// update max and sum
// store dot in local array? We need all 512 scores.
```
For 512 scores, storing in local memory: 512 floats = 2KB per thread. Local memory is global-backed and coalesced if laid out nicely. But 512 scores * 32 threads per warp = 64KB per warp. For 32 warps = 2MB per block. That's the same order as O output, so local memory spill might be OK.
But actually, we can do online softmax WITHOUT storing all scores!
The trick: process keys in groups, compute tile max, update running max, compute exp relative to tile max, accumulate to sum.
For online softmax, we need to know the max to compute exp. If we process one key at a time:
```cpp
float m = -inf;
float l = 0;
float o_local[32] = {0};
for (int k = 0; k < N; ++k) {
float s = 0;
for (int i = 0; i < 32; ++i) {
s += q_local[i] * K[k*D + lane_id*32 + i];
}
warp_reduce(s);
s *= scale;
float new_m = fmaxf(m, s);
float corr = expf(m - new_m);
float p = expf(s - new_m);
l = l * corr + p;
for (int i = 0; i < 32; ++i) {
o_local[i] = o_local[i] * corr + p * V[k*D + lane_id*32 + i];
}
m = new_m;
}
```
Wait! This is BRILLIANT! With warp-level online softmax, each warp processes ONE query. All 32 threads in the warp cooperate:
- Each thread handles a 32-dim chunk of Q, K, V.
- The score is computed as a warp reduction.
- The output is accumulated in each thread's 32-dim local accumulator.
- No shared memory needed!
- No need to store scores!
Let me verify the algorithm carefully:
1. Thread t loads Q[q][t*32:(t+1)*32] into 32 registers.
2. For each key k:
a. Thread t loads K[k][t*32:(t+1)*32] into 32 registers (or directly use).
b. Thread t computes partial dot product: sum_i Q[i] * K[i] for its 32 dims.
c. Warp reduce across 32 threads to get full score.
d. Apply scale.
e. `new_m = max(m, score)`
f. `corr = exp(m - new_m)`
g. `p = exp(score - new_m)`
h. `l = l * corr + p` (but we don't actually need l if we keep unnormalized O)
Wait, we DO need l for the final normalization. Or do we?
Actually, with the standard online softmax formulation:
```
O_new = O_old * exp(m_old - m_new) + exp(s - m_new) * V
l_new = l_old * exp(m_old - m_new) + exp(s - m_new)
```
At the end, O_final = O_new / l_new.
This is correct.
i. Thread t loads V[k][t*32:(t+1)*32].
j. Thread t updates each output dim:
`o_local[i] = o_local[i] * corr + p * v_local[i]`
3. After all keys:
`o_local[i] /= l;` // wait, l is the sum of exp(s_i - m_final)
Actually, O_new already has the exp(s_i - m_new) weights. We just need to divide by l_new.
4. Thread t writes O[q][t*32:(t+1)*32] = o_local.
But wait, there's a subtle issue: `corr = exp(m - new_m)` is the same for all threads in the warp. We should broadcast it. And `p = exp(s - new_m)` is the same for all threads (since s is the same for the whole warp after reduction). We can broadcast p too.
Let's trace through carefully:
```cpp
float q_local[32];
// load Q
for (int i = 0; i < 32; ++i) q_local[i] = Q[q*D + lane*32 + i];
float m = -1e20f;
float l = 0.0f;
float o_local[32] = {0};
for (int k = 0; k < N; ++k) {
float dot = 0.0f;
// load K chunk into k_local? We can load directly.
for (int i = 0; i < 32; ++i) {
float k_val = K[k*D + lane*32 + i];
dot += q_local[i] * k_val;
}
// reduce dot across warp
#pragma unroll
for (int offset = 16; offset > 0; offset /= 2)
dot += __shfl_xor_sync(0xFFFFFFFF, dot, offset);
dot *= scale;
// broadcast m and l? Actually all threads have the same m and l.
// Wait, m and l are scalars, each thread keeps its own copy.
// But they all update identically because dot is the same.
float new_m = fmaxf(m, dot);
float corr = expf(m - new_m);
float p = expf(dot - new_m);
l = l * corr + p;
// update O accumulators
for (int i = 0; i < 32; ++i) {
float v_val = V[k*D + lane*32 + i];
o_local[i] = o_local[i] * corr + p * v_val;
}
m = new_m;
}
// normalize
for (int i = 0; i < 32; ++i) {
o_local[i] /= l;
O[q*D + lane*32 + i] = o_local[i];
}
```
But wait, `m` and `l` are updated identically on all threads, but without synchronization. Is that OK? Let's see:
- `dot` is the same on all threads after the warp shuffle reduction.
- `new_m = fmaxf(m, dot)` where `m` started the same on all threads (and remains the same).
- `corr = expf(m - new_m)` is the same on all threads.
- `p = expf(dot - new_m)` is the same.
- `l = l * corr + p` updates identically.
- `m = new_m` updates identically.
But what about divergence? All threads in the warp execute the same instructions because the data is identical after the shuffle. Actually, no divergence because `m`, `dot`, `new_m` are the same for all lanes. Even though each thread computes its own `dot` locally before the shuffle, after the shuffle all threads have the same `dot` value (because `__shfl_xor` exchanges and adds... wait, no. Let me re-examine.
`dot` is the partial sum for each thread's chunk. After the warp reduction `dot += __shfl_xor(...)`, all threads in the warp hold the TOTAL sum. Yes, because `__shfl_xor` does a butterfly reduction. After all shuffles, each thread has the full sum. So yes, all threads have the same `dot` after reduction.
Then the rest of the computation is identical for all threads. So no divergence!
But there's a register issue: each thread needs:
- q_local[32]
- o_local[32]
- k_val and v_val can be loaded directly without an array, or in a loop.
- So registers ≈ 64 + ~10 temp = 74. That's fine.
Actually, inside the loop, if we load one K element at a time:
```cpp
for (int i = 0; i < 32; ++i) {
float k_val = K[k*D + lane*32 + i];
dot += q_local[i] * k_val;
}
```
`q_local` needs 32 registers. `k_val` is a temp. `dot` is one register. So total ~34 registers.
Then for V accumulation:
```cpp
for (int i = 0; i < 32; ++i) {
float v_val = V[k*D + lane*32 + i];
o_local[i] = o_local[i] * corr + p * v_val;
}
```
`o_local` needs 32 registers. This is fine.
Total registers = 32 (q) + 32 (o) + m + l + corr + p + dot + k_val/v_val + loop counters.
Max ≈ 80 registers. That's well within 255.
Now, memory access pattern:
- Q load: thread lane loads Q[q*D + lane*32 + i] for i=0..31. These are contiguous within the thread, but strided across threads by 32 floats = 128 bytes.
- Thread 0 loads bytes 0..127
- Thread 1 loads bytes 128..255
- ...
- Thread 31 loads bytes 3968..4095
- This is perfectly coalesced! A warp issues one transaction per cache line (128B), covering the full 4096B with 32 threads.
Same for K and V loads.
Same for O writes.
This access pattern is EXCELLENT.
Simplicity: very simple. No shared memory. No __syncthreads.
Performance analysis:
- 1 block = 32 warps = 32 queries.
- Total blocks = 32 * 32 * 16 = 16384.
- 16384 blocks on 132 SMs = 124 blocks per SM if perfectly distributed.
- Each block processes 32 queries. Per query: 512 keys.
- Per warp: for each key, 32 iterations of dot product + 32 iterations of V accumulate.
- Each dot product iteration: 1 FMA + 1 load (K) + 1 load (q_local, already in register).
- Each V iteration: 1 FMA + 2 multiplies + 1 load (V).
- Plus warp shuffle for score.
- Compute per warp: 512 keys * 32 dims * ~5 ops = 81920 ops.
- Memory per warp: Q (32*4B) + K (512*32*4B) + V (512*32*4B) + O (32*4B)
= 128B + 65536B + 65536B + 128B = 131KB per query.
Wait, per warp (1 query): K loaded 512 times * 32 floats * 4 bytes = 65536 bytes. V same. Q and O negligible.
- Total memory per block (32 queries): 32 * 131KB = 4.2MB.
- K and V are the same across all 32 warps in a block! But each warp loads K and V independently.
- However, with L1 cache, the first warp loads K[0]...K[511] and V[0]...V[511] into L1. The other 31 warps read from L1!
- Wait, can L1 serve 31 other warps? The L1 cache on H100 is ~128KB or ~256KB per SM. K+V for 512 keys * 32 dims... no, that's for all keys but only 32 dims per thread.
Actually, K for one key is 4096 bytes. For 512 keys, K = 2MB. That's too big for L1. But each warp only accesses K[k] sequentially. When warp 0 loads K[0], that goes to L1. When warp 1 also loads K[0] (for its query), it might still be in L1. But with 32 warps all loading K sequentially, it's likely that:
- Warp 0 loads K[0] to L1.
- Warp 1 loads K[0] from L1 (or L2).
- Warp 2 loads K[0] from L1 (or L2).
- ...
- This is effectively broadcasting K[0] to all warps once, then they move on to K[1].
The L1 cache might not have enough bandwidth to serve all 32 warps simultaneously for the same address, but it can probably serve a few. They might serialize, but the cache hit means it's much faster than HBM.
However, if 32 warps in a block all try to access the same K[k] at the same time, there will be contention. `__syncwarp()` doesn't help.
Wait, the warps in a block execute independently and might not be perfectly synchronized. They'll naturally stagger.
But here's a key optimization: let ALL warps in a block process the SAME key at roughly the same time. This happens naturally if they iterate over k in lockstep (though not guaranteed). But even without lockstep, the L2 cache can serve K and V.
L2 bandwidth on H100 is much higher than HBM. If K and V are in L2, loading is fast.
For this kernel, total HBM traffic:
- Q: each query loaded once = 2GB
- K: each key loaded once per block? No, each block loads K 32 times (once per warp). 16384 blocks * 2MB = 32GB
- V: same 32GB
- O: 2GB
Total = 68GB
At 3.35TB/s = 20ms. But with L2 caching, most K and V loads come from L2, so HBM traffic is just 2GB for K and 2GB for V (loaded once per key into L2). So effective HBM traffic might be closer to 8GB = 2.4ms.
The 68GB total includes L2 traffic. L2 bandwidth is ~10 TB/s, so 68GB / 10TB/s = 6.8ms.
Compute: each warp does 512 keys * 32 dims * ~8 ops = 131k ops per query. 524288 queries * 131k = 68.7B ops. At 67 TFLOPS = 1ms. But wait, that's the total FLOPs, not just this part. Actually full attention FLOPs is 549B. 549B / 67T = 8.2ms.
Hmm, my per-query calculation might be wrong. Let me recalculate:
For one query:
- Scores: 512 keys * 1024 dims = 524k FMAs.
- V accumulation: 512 keys * 1024 dims = 524k FMAs.
- Total: ~1M FMAs = 2M FLOPs.
- 524288 queries * 2M = 1.05 T FLOPs.
- At 67 TFLOPS = 15.7ms.
My earlier 549B came from 2*B*H*N^2*D = 2*32*32*512*512*1024 = 549B. Yes, that's correct.
Each warp has 32 threads. The partial dot product for 32 dims per thread requires 32 FMAs. Then warp shuffle reduction (negligible). So 32 FMAs per thread per key for scores + 32 FMAs per thread per key for V = 64 FMAs.
Total FMAs per thread = 512 keys * 64 FMAs = 32768 FMAs.
Total threads = 16384 blocks * 32 warps * 32 threads = 16.7M threads.
Total FMAs = 16.7M * 32768 = 549B. Matches.
At 67 TFLOPS: 549B / 67T = 8.2ms.
At L2-bound 68GB / 10TB = 6.8ms.
So total time might be around 10-15ms.
Is the reference faster or slower?
If reference uses cuBLAS for two GEMMs, each GEMM [1024, 512, 1024, 512] is actually:
- GEMM 1: [1024 batches, M=512, N=512, K=1024]
- This is large enough for cuBLAS to achieve good efficiency. cuBLAS might get 80% of peak on H100.
- But each GEMM only has 512x512 output, which is relatively small.
- cuBLAS batch GEMM: 274B FMAs each = 4ms each at 67 TFLOPS.
- Plus softmax: ~1ms.
- Plus materializing S: extra memory traffic.
Actually, for a batch size of 1024, a [512, 512, 1024] GEMM total FLOPs = 1024 * 512 * 512 * 1024 * 2 = 549B. cuBLAS on H100 might achieve say 50-100 TFLOPS for this size? 549B / 80TF = 6.8ms per GEMM. Two GEMMs = 13.6ms. Plus memory for S = maybe +2ms. Total ~16ms.
If I write a fused kernel taking ~10-15ms, I might be competitive. But can I do better than 10ms?
What if I use shared memory to explicitly prefetch K and V tiles, avoiding redundant loads entirely within a block?
With 32 warps in a block all processing the same K and V:
- K tile: 32 keys * 1024 dims = 128KB. Wait, BK=32, D=1024, that's 128KB.
- With BK=32, we process 512/32 = 16 key tiles.
- For each key tile, all 32 warps in the block need the same K[kt:kt+31][:] and V[kt:kt+31][:].
- Store K tile and V tile in shared memory: 2 * 128KB = 256KB. Exceeds 227KB limit.
What if BK=16? K tile = 16*1024*4 = 64KB. V tile = 64KB. Total = 128KB. Plus smem_q for 32 queries? But in this warp design, each warp handles 1 query and loads its own Q into registers. So no need for Q in smem.
Total smem = 128KB. Under 227KB. Good!
With BK=16:
- Number of key tiles = 512/16 = 32.
- Each block loads each K tile once into smem (instead of 32 times via L2).
- This reduces L2 traffic significantly, and ensures perfect reuse.
Wait, how do we load K tile into smem?
- Block has 1024 threads. K tile is 16*1024 = 16K floats.
- Each thread loads 16 floats into smem. Perfectly coalesced.
- __syncthreads() to ensure all K is loaded.
- Then all 32 warps compute their dot products using smem_k.
- Then all 32 warps do V accumulation using smem_v.
- __syncthreads() before loading next tile.
This is a much better use of shared memory! Let me design this:
Block = 1024 threads? Actually 32 warps = 1024 threads. That's the max.
But wait, with 32 warps and 1024 threads per block, and 2 blocks per SM (2048 threads), that's fine. Occupancy might be limited.
Actually, for Ampere/Hopper, max threads per SM = 2048. So with 1024 threads per block, we can have 2 concurrent blocks per SM. That gives good occupancy.
But 1024 threads means register pressure matters. 1024 threads * 80 registers = 81K registers per SM. But an SM only has 64K registers. So we'd be limited to 1 block per SM if each thread uses 64 registers. With 80 registers, 1024*80 = 81K > 64K, so the compiler might spill or we'd only fit 1 block = 1024 threads per SM.
What if we use 512 threads per block (16 warps)?
- 512 threads * 80 regs = 40K registers. We can fit 1 block (40K) comfortably, and maybe 1.5 but effectively 1 block per SM.
- Or if we use 64 registers: 512*64 = 32K, so 2 blocks = 64K, perfect!
- But can we keep register usage under 64?
Let's count:
- q_local: 32 floats
- o_local: 32 floats
- loop temp: 3 floats
- m, l, corr, p, s, dot: 6 floats
Total = 73 floats. That's a lot. The compiler might keep q_local and o_local in registers, but might spill some.
Actually, q_local and o_local are arrays. The compiler stores them in registers if possible. 64 registers is tight.
Alternative: process Q, K, V in chunks of 16 instead of 32:
- q_local[16], o_local[16]
- Total registers: 16 + 16 + 6 = 38. Very comfortable!
- But we need 1024/16 = 64 chunks. More loop iterations.
- D chunk = 16 means each thread handles 16 dims.
- 512 keys * 64 chunks * 2 FMAs = 65536 FMAs. Same as before.
With D chunk = 16, and 512 threads per block:
- Q load: 512 threads * 16 floats = 8192 floats per block. Coalesced.
- Block processes 16 queries? (16 warps, 1 query per warp)
- Wait, 512 threads = 16 warps. Each warp handles 1 query. So block handles 16 queries.
- Total blocks = 32 * 32 * (512/16) = 32*32*32 = 32768 blocks.
This is a lot of blocks. 32768 / 132 = 248 blocks per SM. Each block does 65536 FMAs per thread. 512 threads * 65536 = 33.5M FMAs per block. Total FMAs = 1.1T? Wait, total work should be 549B FMAs. 32768 blocks * 512 threads * 65536 FMAs = 1.1T FMAs. That's twice as much.
Oh wait, I miscalculated. For one query, the total FMAs is 512 keys * 1024 dims * 2 (score + V) = 1,048,576 FMAs.
With chunk size 16: 512 keys * 64 chunks * 16 dims * 2 = same.
Per thread: for one key, chunk loop does 16 iterations. 16 FMAs for score + 16 for V = 32 FMAs.
For 512 keys: 512 * 32 = 16384 FMAs per thread.
Total threads = 32768 blocks * 512 threads = 16.7M threads.
Total FMAs = 16.7M * 16384 = 274B. That's only half of 549B.
Why? Because 16384 FMAs per thread * 16.7M threads = 274B. But we need 549B.
Ah, I see. The score computation alone is 512 keys * 1024 = 524k FMAs. But if 512 threads process 512 keys... no, in my design, one thread processes ONE query, not one key. So for one query with 512 keys and D=1024:
- Score: 512 * 1024 = 524k FMAs.
- V: 512 * 1024 = 524k FMAs.
- Total: ~1M FMAs.
With chunk size 16:
- Score: 512 keys * 64 chunks * 16 FMAs = 524k FMAs.
- V: same.
Per warp (1 query): total FMAs = ~1M.
Per block (16 warps/queries): 16M FMAs.
Total blocks = 32*32*32 = 32768.
Total FMAs = 32768 * 16M = 524B. Yes, that matches.
Per thread FMAs: 1M FMAs / 32 threads per warp = 32768 FMAs per thread. (Actually, 1M / 32 = 32k)
16.7M threads * 32k FMAs = 535B. Close enough.
OK so D chunk size doesn't change total work. It just changes the loop structure.
Now, with shared memory for K and V:
- Block = 16 warps = 512 threads, processing 16 queries.
- K tile: BK=16, D=1024, so smem_k = 16*1024*4 = 64KB.
- V tile: same 64KB.
- Total: 128KB.
- With 227KB limit per block, and 64K registers per SM, we might fit 1 block per SM (128KB smem, ~40K registers).
- Or if we use fewer threads per block...
With 512 threads per block and 128KB smem:
- Max concurrent per SM: limited by smem to 1 block (128KB) or 2 blocks (256KB > 227KB, no). So 1 block per SM.
- Total concurrent blocks = 132.
- Waves = 32768 / 132 = 248 waves.
- That's a lot of waves! If each block takes a short time, the overhead of launching all those waves adds up.
What if we reduce the number of blocks by increasing queries per block?
Block = 64 warps = 2048 threads? Exceeds max 1024. Not possible.
Block = 32 warps = 1024 threads. Then block handles 32 queries.
Total blocks = 32*32*16 = 16384. Waves = 16384/132 = 124.
Still many waves, but better.
Block = 64 queries? Not possible with 1 thread per query (max 1024 threads). Unless we have 1 block per 64 queries with 16 threads per query. But 64*16=1024 threads. Then each query is handled by 16 threads cooperatively. That changes things.
Actually, the cooperative approach with 1 warp per query and multiple warps per block seems clean. Let's commit to:
- 1 block = 32 warps = 1024 threads
- Block handles 32 queries
- Total blocks = 32 * 32 * 16 = 16384
- Uses shared memory for K and V tiles if we want.
But register pressure with 1024 threads is a concern. If each thread needs 64 registers, total = 65536. An SM has 65536 registers. So exactly 1 block per SM. This works!
With 1 block per SM = 1024 threads per SM. H100 SMs can handle up to 2048 threads, so occupancy is 50% in terms of threads, but 100% in terms of blocks (1 block is the max due to smem). Actually, 50% thread occupancy is fine.
16384 blocks / 132 SMs = 124 waves. Each wave of 132 blocks. At per-block time of, say, 0.5ms, total = 62ms. That's too long.
Wait, 0.5ms per block is way too high. Let's estimate:
Each warp does ~1M FMAs. 1024 threads = 32M FMAs per block. At 512 GFLOPS per SM: 32M / 512G = 62 microseconds. 124 waves * 62us = 7.7ms.
That seems reasonable! Plus memory loads. 7.7ms is competitive.
But can we reduce waves? With 1024 threads per block processing 32 queries, it's 16384 blocks. If we could process 64 queries per block, total blocks = 8192, waves = 62. But we'd need 2048 threads, which exceeds the 1024 limit.
Alternatively, process 32 queries with 512 threads (16 warps). Then 2 blocks per SM concurrently if registers allow.
If 512 threads * 64 regs = 32K registers, then 2 blocks per SM = 64K registers. OK.
And smem per block = 128KB. 2 blocks = 256KB > 227KB. Not OK. So we'd need to reduce smem.
If we don't use shared memory and instead rely on L2 cache, we can have 2 blocks per SM with 512 threads each. Total concurrent = 264. Waves = 32768 / 264 = 124. Same issue.
If we don't use shared memory and use 256 threads (8 warps) processing 8 queries:
- Total blocks = 32*32*64 = 65536.
- Waves = 65536 / (132 * 8) ??? With 256 threads * 64 regs = 16K regs. Up to 4 blocks per SM = 64K regs. So 4 blocks = 1024 threads.
- Concurrent blocks = 132 * 4 = 528.
- Waves = 65536 / 528 = 124.
It seems like regardless of block size, if we process all queries with 1 thread per query, we get ~124 waves. The key is that total parallel work is 524288 queries, and SMs can run ~1000-2000 threads. With 132 SMs * ~1500 threads = 198k threads concurrent. 524k threads / 198k = 2.6 waves.
Wait, I think my wave calculation was wrong. Total threads = 524288.
Concurrent threads across all SMs = 132 SMs * threads_per_SM.
If threads_per_SM = 1024 (1 block of 1024), total concurrent = 132 * 1024 = 135k.
Waves = 524k / 135k = 3.88 waves.
If threads_per_SM = 2048 (2 blocks of 1024 or 4 blocks of 512), total concurrent = 132 * 2048 = 270k.
Waves = 524k / 270k = 1.94 waves.
So we're looking at 2-4 waves. That's very good! I was dividing blocks by SMs before, but it should be total threads by concurrent threads.
For shared memory limited to 1 block per SM (1024 threads), waves = 524k / (132*1024) = 3.88.
For 2 blocks per SM (if possible), waves = 2.
For 512 threads with 2 blocks = same 1024 threads per SM. Waves = 3.88.
Actually 512 threads * 2 blocks = 1024 threads per SM. Same as 1 block of 1024.
To get to 2048 threads per SM:
- 2 blocks of 1024: each block would need < 113KB smem. With BK=16 and K+V=128KB, we can't.
- 4 blocks of 256: each block needs < 56KB smem. We can if we don't use smem or use very small tiles.
- 4 blocks of 512: need < 56KB smem, and 512*128 regs? 512*64=32K, 4 blocks=128K > 64K. Not possible.
So max concurrency is likely 1 block of 1024 threads, or 2 blocks of 512 with no shared memory. Either way, ~1024-2048 threads per SM. Let's say 1024 for safety. Waves = 3.88.
Total runtime ≈ 4 * block_time.
Block_time for 32 warps processing 32 queries:
- Each warp: 1M FMAs. At 50% peak (due to memory), ~0.1ms.
- 4 waves * 0.1ms = 0.4ms. That's very fast!
But this assumes perfect utilization. Let's be more conservative: 2ms.
Now, how does the reference do? The reference with cuBLAS + softmax might take 5-15ms depending on backend.
So if I can write a clean fused kernel with no shared memory and good register usage, I might actually beat the reference.
Let me think about whether to use shared memory or not.
- Without shared memory: relies on L2 cache for K/V reuse among warps. Clean, no sync overhead.
- With shared memory: explicit reuse, but limits occupancy.
For N=512, D=1024, K and V for BK=16 keys = 128KB shared memory. This is fine if we use 512 threads per block and only 1 block per SM.
Actually, I just realized: with shared memory for K and V, we MUST synchronize threads. But each warp handles one query, and the score computation for each query needs ALL K dimensions. So we can load a chunk of K[0:BK][0:DT] into smem, then all warps compute their partial scores using this chunk. But partial scores must be accumulated across chunks. This means we need a per-query score accumulator. But the score is per query per key, not per query per chunk. Wait, no: the score is the dot product of Q and K, which is sum over D. If we compute partial dot products for a chunk of D, we accumulate them.
Actually, for BK=16 and D chunk = 64:
- Load K[0:16][0:64] into smem: 16*64 = 1024 floats = 4KB.
- All 32 warps load their Q[q][0:64] into registers.
- Warp q computes partial scores for q and keys 0..15 using chunk 0..63.
- Then load next chunk of K and Q.
- After all chunks, we have accumulated the full scores.
This doesn't require a score array. Each warp just accumulates s_k for each key k:
```cpp
float s[16]; // registers
for key_tile:
for d_chunk:
load K chunk
for k=0..15:
s[k] += q[d] * K[k][d];
```
But `s[16]` per warp = 16 registers. And we need 32 lanes computing. If each lane handles one key, then lane k handles s[k] for ONE query. So each lane needs only 1 register per query. But we have 32 lanes and 16 keys. Some lanes are idle.
If each warp handles 1 query and BK=16:
- 32 lanes, 16 keys. Each lane handles 2 keys? Or one key with 2 lanes?
- Actually, 16 keys and 32 lanes means each key gets 2 lanes. 2 lanes can split D.
- Lane 0,1 handle key 0, splitting D in half.
- This uses all threads and parallelizes dot product!
But then we need shuffle to reduce across the 2 lanes for each key.
This is getting complicated. The no-shared-memory approach with warp-level reductions and 32 dims per thread is much simpler and likely performance-sufficient.
Let me finalize on this approach: one-warp-per-query, no shared memory.
Wait, I think I should try both the no-shared-memory version and a shared-memory version. But I need to write the code first.
Actually, let me reconsider: with the no-shared-memory version, each warp loads K independently. L2 bandwidth is used for every warp. L2 bandwidth on H100 is around 10-20 TB/s. HBM is 3.35 TB/s.
If K and V don't fit in L2 for reuse (because many different heads compete), traffic goes to HBM and we're very slow.
With shared memory, K is loaded explicitly once per block and reused by all 32 warps. This guarantees we read K from HBM only once per block, and then all 32 warps use it from shared memory.
For V, same thing.
So shared memory is safer and gives explicit control over memory traffic.
I'll design a kernel with:
- 1 warp per query
- Block = 32 warps = 1024 threads
- Block processes 32 queries
- Shared memory for K and V tiles
- BK=16 keys per tile
- D processed in chunks, but wait...
Actually, for the dot product with shared memory, let me reconsider the structure:
If BK=16 and each warp handles 1 query:
- We want to compute 16 scores per query per tile.
- A warp has 32 lanes. We can have each lane compute 1 score using half the lanes, or have lanes cooperate (2 lanes per key).
With 2 lanes per key and BK=16: 16*2 = 32 lanes. Perfect!
- Lane pair (2k, 2k+1) computes score for key k.
- They split D in half. Lane 2k handles D[0:512], lane 2k+1 handles D[512:1024].
- After partial dot products, shuffle reduce lane 2k and 2k+1.
- Now all 32 lanes have all 16 scores? No, lane 2k has score for key k.
To have all lanes know all scores for softmax, we can shuffle. But with only 32 lanes and 16 keys, we can store scores in 16 lanes and have the other 16 lanes idle during softmax? Or we can distribute scores and do warp-level operations.
Actually, with 16 scores, we can have:
- lane 0..15 store scores for keys 0..15.
- lanes 16..31 are idle during score storage.
Then for softmax:
- Max over 16 lanes: `__shfl_xor` can reduce 16 values if they're in lanes 0..15. But `__shfl_xor` operates on all 32 lanes.
- We can do: lanes 0..15 have scores, lanes 16..31 have -inf.
- Warp max = fmaxf across all 32 lanes = max of the 16 scores.
- Similarly for sum.
Then for V accumulation: all 32 lanes can participate, each handling 1024/32 = 32 dims of output.
- For key k (0..15), we need p_k for that key.
- p_k is stored in lane k.
- Lane k broadcasts p_k to all lanes.
- Then all lanes update their 32 dims for this key.
This works! And all lanes are used efficiently for V accumulation.
Register usage per lane:
- Q: 1024/2 = 512 floats if we split D... wait, we don't need to keep all Q in registers.
We can process D in chunks of DT (say 64), and load Q and K for each chunk.
```cpp
float q_local[32]; // DT=64, split across 2 lanes = 32 elements per lane.
float s0 = 0; // score for key handled by this lane
```
Actually, with 2 lanes per key, DT could be the full 64 of a chunk.
- For D chunk 0..63:
- All lanes load Q[0..63]? No, each lane loads only its portion.
- Wait, Q for the warp's query needs to be accessible to all lanes.
- Since we split D across 2 lanes, each lane loads a sub-chunk.
For chunk dt (64 elements):
- Lane 2k loads Q[dt:dt+32]
- Lane 2k+1 loads Q[dt+32:dt+64]
- K[k][dt:dt+32] and K[k][dt+32:dt+64] loaded from shared memory.
Actually, K is in shared memory, so all lanes can access it. We'll load K[0:16][dt:dt+64] into smem.
For simplicity, let me process D in chunks of 32 to match the warp size:
- DT = 32.
- Each lane loads exactly one float of Q and one float of K per chunk? No, we want each thread to process a contiguous chunk.
Wait, instead of this complexity, let me just go with the no-shared-memory warp design I described earlier. It's simpler and avoids all these tile complications. And with online softmax, it's correctly fused.
Final design: one-warp-per-query, no shared memory.
- Block = 32 warps = 1024 threads.
- Grid = (32*32, 16) where x = batch*head, y = query_group (512/32 = 16).
- Each warp handles one query. 32 warps handle 32 queries in a block.
Actually, in CUDA, blocks are limited to 1024 threads, and grids are up to 2^31-1 in x, 65535 in y and z.
I can use grid = (32*32*16, 1) = 16384 blocks, each with 32 warps.
Or to be cleaner: grid = (32, 32, 16), block = (32, 32, 1)?
Wait, block can have 3D thread indexing. blockDim.z is up to 64. blockDim.x * blockDim.y * blockDim.z <= 1024.
If blockDim = (32, 32, 1) = 1024 threads.
Then grid = (32, 32, 16).
Within block:
- threadIdx.x = lane_id (0..31)
- threadIdx.y = warp_id (0..31), representing query offset within the block
So query_id = blockIdx.z * 32 + threadIdx.y.
head_id = blockIdx.y.
batch_id = blockIdx.x.
Wait, lane_id = threadIdx.x, and warp_id = threadIdx.y. But threadIdx.x=32 and threadIdx.y=32 means it's 2D. CUDA will create 32 warps, each of 32 threads. But the warp grouping is by threadIdx.x being the fastest varying dimension, so a warp consists of threads (0..31, y) for fixed y. Yes, that's correct.
So within warp y, all threads have threadIdx.x = 0..31. This is one warp.
And blockIdx.z tells which group of 32 queries in the sequence.
In this configuration:
- threadIdx.x determines which slice of D the thread handles.
- Each thread handles D/32 = 32 dimensions (since D=1024 and 32 threads).
For accessing Q: Q[b][h][q][d]
- d start = threadIdx.x * 32.
- Thread loads Q[q][d_start + i] for i=0..31. Coalesced.
For K and V: same pattern.
For dot product and accumulation, we use warp-level shuffle for the score.
Let me write this kernel out carefully. I need to be careful with numerical precision and correctness.
```cpp
__global__ void attention_kernel(
const float* __restrict__ Q,
const float* __restrict__ K,
const float* __restrict__ V,
float* __restrict__ O,
int B, int H, int N, int D)
{
int batch = blockIdx.x;
int head = blockIdx.y;
int q_group = blockIdx.z;
int lane = threadIdx.x;
int warp = threadIdx.y;
int q = q_group * 32 + warp;
if (q >= N) return;
int head_offset = ((batch * H + head) * N) * D;
int q_offset = head_offset + q * D;
float scale = 1.0f / sqrtf((float)D);
// Load Q: 32 floats per thread
float q_local[32];
#pragma unroll
for (int i = 0; i < 32; ++i) {
q_local[i] = Q[q_offset + lane * 32 + i];
}
float o_local[32] = {0.0f};
float m = -1e20f;
float l = 0.0f;
// Iterate over all keys
for (int k = 0; k < N; ++k) {
int k_offset = head_offset + k * D;
// Compute dot product Q·K
float dot = 0.0f;
#pragma unroll
for (int i = 0; i < 32; ++i) {
float k_val = K[k_offset + lane * 32 + i];
dot += q_local[i] * k_val;
}
// Warp-level reduction
#pragma unroll
for (int offset = 16; offset > 0; offset /= 2) {
dot += __shfl_xor_sync(0xFFFFFFFF, dot, offset);
}
dot *= scale;
// Online softmax update
float new_m = fmaxf(m, dot);
float corr = expf(m - new_m);
float p = expf(dot - new_m);
// Load V and update output
#pragma unroll
for (int i = 0; i < 32; ++i) {
float v_val = V[k_offset + lane * 32 + i];
o_local[i] = o_local[i] * corr + p * v_val;
}
l = l * corr + p;
m = new_m;
}
// Normalize and write output
float norm = 1.0f / l;
#pragma unroll
for (int i = 0; i < 32; ++i) {
O[q_offset + lane * 32 + i] = o_local[i] * norm;
}
}
```
Wait, I need to check something critical about this kernel.
The issue: `l = l * corr + p`. This maintains the sum of exp(s_i - m_current). At the end, m_current = max over all keys, so l = sum exp(s_i - max). And we normalize by l. This is correct.
But what about `corr = exp(m - new_m)`? If m = -inf initially, and new_m = dot for the first key:
- corr = exp(-inf - dot) = 0
- Then o_local = o_local * 0 + p * v = p * v.
- l = 0 + p = exp(dot - dot) = 1.
This is correct.
For the second key, if dot < m:
- new_m = m (old max)
- corr = exp(m - m) = 1
- p = exp(dot - m)
- o_local = o_local + p * v
- l = l + p
Correct.
If dot > m:
- new_m = dot
- corr = exp(m - dot)
- p = exp(dot - dot) = 1
- o_local = o_local * corr + 1 * v
- l = l * corr + 1
Correct.
However, there's a subtle precision issue. When m updates, we re-weight all previous output by `corr`. This is done elementwise with `o_local[i] = o_local[i] * corr + p * v`. This is the standard online softmax formulation. It should be correct.
Now, the interesting point: `l` is only used at the end for normalization. But we could avoid accumulating l by keeping a running sum. Actually, l = sum exp(s_i - m_current). Yes.
Is there any issue with `__shfl_xor_sync` and the loop unrolling? It should be fine.
One concern: warp divergence due to `if (q >= N) return;`. With N=512 and q_group=0..15, q = q_group*32 + warp. For q_group=15, q = 480..511. All valid since N=512. So no divergence for this problem size.
Another concern: memory alignment. Q, K, V are allocated by torch.rand, which uses aligned memory. q_offset + lane*32 + i: lane*32 is 128-byte aligned for lane=0, 128-byte aligned for lane=1? No, lane*32 + i for i=0 gives q_offset + lane*32. Consecutive lanes access q_offset + 0, 32, 64, 96, ... which are 128 bytes apart. This is NOT coalesced for 4-byte loads.
OH NO! I made a critical mistake.
For coalesced access, thread t should access address `base + t` (or `base + t * stride` where stride is contiguous). Here, lane 0 accesses offset 0, lane 1 accesses offset 32. If each thread loads a single float, lane 0 accesses byte 0, lane 1 accesses byte 128. These are 128 bytes apart. A warp of 32 threads issues 32 separate 128-byte cache line requests. That's terrible!
To be coalesced, thread t should access `base + t`. Then we'd load Q[q][lane + i*32] instead of Q[q][lane*32 + i].
Let me think: if I want each thread to load 32 contiguous floats, the layout should be:
Thread t loads: `Q[offset + t], Q[offset + t + 32], Q[offset + t + 64], ...`
Wait, no. For coalescing, at a given instruction, all threads in a warp should access consecutive memory addresses OR the same address.
If I do:
```cpp
for (int i = 0; i < 32; ++i) {
float val = Q[offset + lane + i * 32];
}
```
At i=0: threads access offset+0, offset+1, offset+2, ..., offset+31. Perfectly coalesced (one 128-byte transaction).
At i=1: threads access offset+32, offset+33, ..., offset+63. Perfectly coalesced.
...
At i=31: threads access offset+992, ..., offset+1023. Coalesced.
This is perfect! And each thread ends up with 32 floats, but they're not contiguous in D; they're strided by 32.
Does this matter for the dot product?
```cpp
dot += q_local[i] * k_local[i];
```
where q_local[i] = Q[q][lane + i*32] and k_local[i] = K[k][lane + i*32].
This computes: sum over i of Q[q][lane + i*32] * K[k][lane + i*32].
After warp shuffle reduce, we get sum over all lanes and all i of Q[q][lane + i*32] * K[k][lane + i*32] = sum over all d of Q[q][d] * K[k][d].
Yes, this is correct! The order doesn't matter for the dot product.
And for the output accumulation with V:
```cpp
o_local[i] += p * v_local[i];
```
where v_local[i] = V[k][lane + i*32].
This updates O[q][lane + i*32].
At the end, we write O[q][lane + i*32] back.
Again, the writes are coalesced because at each i, consecutive threads write consecutive addresses.
Perfect! This is the correct memory pattern.
Let me rewrite the kernel with this pattern. The performance should be much better.
One more optimization: float4 loads and stores. Since each warp load is perfectly coalesced and covers 128 bytes (32 floats), we can use float4 to load 4 floats per instruction.
```cpp
float4 q_local[8]; // 8 * 4 = 32 floats
for (int i = 0; i < 8; ++i) {
q_local[i] = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 8];
}
```
Wait: `reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 8]` means:
- Base pointer: Q + q_offset, cast to float4*
- Index: lane + i * 8
- float4 at address = base + (lane + i*8) * sizeof(float4) = base + (lane + i*8) * 16 bytes.
- For lane=0..31 and i=0..7:
- i=0: offsets 0..31 in float4 units = bytes 0..496 (32 * 16)
- i=1: offsets 32..63 = bytes 512..1008
- ...
- i=7: offsets 224..255 = bytes 3584..4088
- Total: 256 float4s = 1024 floats = 4096 bytes.
- Consecutive lanes access consecutive float4s at each i. Perfectly coalesced.
This is equivalent to thread t accessing floats at indices t, t+32, t+64, ..., t+992.
Because float4 index j corresponds to floats 4j, 4j+1, 4j+2, 4j+3.
If lane loads float4 at index `lane + i*8`, that covers floats:
4*(lane + i*8) through 4*(lane + i*8)+3
= 4*lane + 32*i through 4*lane + 32*i + 3
Wait, thread lane loads a float4 that contains 4 floats, but we don't want that if each thread should handle individual floats for the dot product. We'd need to unpack the float4.
Actually, using float4 is fine if we unpack:
```cpp
float4 q4 = ...;
q_local[i*4+0] = q4.x;
q_local[i*4+1] = q4.y;
q_local[i*4+2] = q4.z;
q_local[i*4+3] = q4.w;
```
But this creates extra move instructions. It could still be beneficial because memory bandwidth improves.
Alternatively, we can restructure the dot product to use float4 for FMAs too. But then we'd need to do 4 FMAs per float4, which is fine.
Actually, let me just keep it simple with float loads first, and switch to float4 if needed. The compiler might vectorize automatically.
Wait, D=1024 and we process 32 chunks of 32 floats per thread. That's 32 iterations. Using float would be 32 memory instructions. Using float4 reduces to 8 memory instructions. This is a significant reduction.
Let me use float4:
```cpp
#define DIM_PER_THREAD 32
#define FLOAT4_PER_THREAD (DIM_PER_THREAD / 4)
float4 q4[FLOAT4_PER_THREAD];
#pragma unroll
for (int i = 0; i < FLOAT4_PER_THREAD; ++i) {
q4[i] = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
}
float q_local[32];
#pragma unroll
for (int i = 0; i < FLOAT4_PER_THREAD; ++i) {
q_local[i*4+0] = q4[i].x;
q_local[i*4+1] = q4[i].y;
q_local[i*4+2] = q4[i].z;
q_local[i*4+3] = q4[i].w;
}
```
Wait, `lane + i * 32` for float4 indexing:
- `lane` ranges from 0 to 31.
- `i` ranges from 0 to 7.
- Index = lane + i*32.
- For i=0: float4 indices 0..31. These are consecutive float4s.
- For i=1: float4 indices 32..63.
- ...
- Total float4s accessed per warp = 32 * 8 = 256 float4s = 1024 floats.
- Memory layout: Q is float array. Cast to float4*. Accessing index j gets bytes j*16 through j*16+15.
- Thread lane accesses index lane + i*32. Since 32 threads access indices 0..31 for i=0, they access bytes 0..511. Perfectly coalesced.
- For i=1, indices 32..63, bytes 512..1023. Coalesced.
- And so on.
This is correct. And it loads 8 float4s per thread.
Now for the inner loop over K and V:
```cpp
for (int k = 0; k < N; ++k) {
float4 dot4 = make_float4(0,0,0,0);
// We need to compute partial dot product using float4
#pragma unroll
for (int i = 0; i < FLOAT4_PER_THREAD; ++i) {
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot4.x += q_local[i*4+0] * k4.x;
dot4.y += q_local[i*4+1] * k4.y;
dot4.z += q_local[i*4+2] * k4.z;
dot4.w += q_local[i*4+3] * k4.w;
}
float dot = dot4.x + dot4.y + dot4.z + dot4.w;
// warp reduce...
```
Wait, I can keep q in float array, but load K as float4 and unpack. Or load both as float4 and keep in float4 for dot product if we reorganize.
Actually, we can avoid unpacking Q and store it as float4 too, but then the dot product needs element-wise multiply. With float4, we can't directly do `q4[i] * k4[i]` because CUDA doesn't have built-in float4 multiply-add. We'd need to do `.x * .x + .y * .y + ...` manually, which is what I wrote.
Alternatively, we can keep everything as floats and let the compiler handle vector loads. Modern compilers (nvcc with -O3) often automatically vectorize loads when they see coalesced float loads in a loop with known bounds.
Let me try a simpler version first with float loads, and then consider float4 if performance is needed. Actually, for correctness and clarity, I'll use float4 for memory loads but process as floats.
Let me also think about whether I should use shared memory. With float4 loads, memory bandwidth is already good. And using no shared memory means more threads per SM. I think the no-shared-memory version is worth trying. The kernel is extremely simple.
One more optimization: the K and V loads in the inner loop could be vectorized with float4.
Let me write the simplified kernel:
Block shape: (32, 32, 1) = 1024 threads.
Grid shape: (32, 32, 16) = 16384 blocks.
Wait, gridDim.z max is 65535. 16 is fine.
But `blockIdx.z` is used for query groups.
However, I need to map addresses correctly. The input layout is [B, H, N, D] with default strides.
In PyTorch, default stride for [B, H, N, D] is [H*N*D, N*D, D, 1].
So offset = ((b * H + h) * N + n) * D + d.
Let me write the full kernel code.
Actually, I should also consider using `__ldg` for read-only data to help caching.
And for the output write, we can just use regular store.
Let me also think about numerical precision. The reference uses fp32. My kernel uses fp32 FMAs and `expf`. `torch.nn.functional.scaled_dot_product_attention` might use a different softmax algorithm or different precision for intermediate accumulations. But online softmax in fp32 should be accurate enough.
One more thought: `expf` can be expensive. With 512 keys and 32 warps, each warp does 512 `expf` calls. `expf` takes dozens of cycles. This might be a bottleneck.
We have 512 keys * 16384 warps = 8.4M `expf` calls. Each `expf` might take ~20-30 cycles. At 2GHz, that's about 10-15ms across all SMs. Not too bad, but significant.
Can we use fast `__expf`? The `--use_fast_math` flag should enable it. But it might sacrifice precision. Let's try it and see if correctness passes.
Actually, the reference uses standard `exp`. With fast math, the error might be too large. But fp32 tolerance is 1e-4. `__expf` vs `expf` difference is usually around 1e-5 or so. It might be fine.
But let's first try with `expf` and see if we pass correctness. If yes, we can try `__expf` for speed.
Wait, the first kernel already passed correctness. This new kernel uses a different algorithm but should give the same result.
Wait, there's a subtle difference. In standard attention:
- S = Q @ K^T
- P = softmax(S / sqrt(D))
- O = P @ V
In my online softmax formulation:
- I compute dot products and scale.
- Then I apply online softmax to get normalized weights.
- Then O = weighted sum of V.
This should give exactly the same result as the standard formulation (up to floating point associativity differences).
But does PyTorch's `scaled_dot_product_attention` use exactly the same online softmax? No, but with fp32 tolerance of 1e-4, minor differences from associativity should be acceptable.
Let me think about register pressure. With `q_local[32]`, `o_local[32]`, and the loop variables, it's about 68 registers. But the compiler might keep `q_local` in registers. With unrolled loops, it definitely will.
However, 68 registers + some stratch = ~80. For 1024 threads, that's 80K > 64K. The compiler might have to spill.
But we can avoid `q_local[32]` by loading Q on the fly:
```cpp
float q_val = Q[q_offset + lane + i * 32];
dot += q_val * k_val;
```
This way, we only need a scalar `q_val` per iteration, not an array. The compiler might not keep it in a register between iterations, or it might reload from L1 cache.
Actually, with warp-level code and 32 iterations, the compiler will likely keep the Q values in registers or at least keep around. But there's no guarantee. If Q gets spilled to local memory, that's bad.
To help the compiler, we can load Q into a smaller array and explicitly unroll:
```cpp
float q00, q01, q02, ..., q31;
// load
// use in loop
// or use float4
float4 q4[8];
#pragma unroll
for (int i = 0; i < 8; ++i) q4[i] = loadQ(...);
// in K loop:
float4 k4[8];
#pragma unroll
for (int i = 0; i < 8; ++i) {
k4[i] = loadK(...);
dot += q4[i].x * k4[i].x;
dot += q4[i].y * k4[i].y;
dot += q4[i].z * k4[i].z;
dot += q4[i].w * k4[i].w;
}
```
With float4, the compiler has fewer variables to track. But the register usage for 8 float4s = 32 registers.
For the V accumulation, we need `o_local` as float or float4:
```cpp
float4 o4[8] = {0};
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 v4 = loadV(...);
o4[i] = o4[i] * corr + v4 * p;
// but float4 doesn't support * for scalars directly
// we need: o4[i].x = o4[i].x * corr + v4.x * p;
}
```
Hmm, float4 arithmetic with scalars isn't built-in. I'd have to write helper functions or manual arithmetic.
Maybe I should stick with plain float arrays. Let the compiler optimize.
Actually, another option: use `__half2` or vector types? No, we're doing fp32.
I think the simplest correct approach is:
```cpp
float q_reg[32];
for (int i = 0; i < 32; ++i) q_reg[i] = ...;
float o_reg[32] = {0};
for (int k = 0; k < N; ++k) {
float dot = 0;
for (int i = 0; i < 32; ++i) {
dot += q_reg[i] * K[k_offset + lane + i * 32];
}
// reduce
// softmax
for (int i = 0; i < 32; ++i) {
o_reg[i] = o_reg[i] * corr + p * V[k_offset + lane + i * 32];
}
}
```
NVCC will likely keep `q_reg` in registers. With 32 registers and `o_reg` with 32 registers, plus loop temps, that's about 70 registers. It might spill a little, but with unrolled loops it should be OK.
Let me write this kernel and compile it.
One more thing: `blockIdx.z` is the third dimension of the grid. `gridDim.z` must be <= 65535. 16 is fine.
Wait, can I use `blockIdx.z` in `load_inline`? Yes, standard CUDA.
Let me write the complete Python code with the kernel.
Actually, I realize something. With 1024 threads per block and only `threadIdx.x` and `threadIdx.y`, blockDim = (32, 32, 1). A warp consists of threads with the same `threadIdx.y` and consecutive `threadIdx.x` values. The hardware groups threads by `threadIdx.x` first. So threads (0..31, y) form a warp. This is correct.
But note: within a warp, consecutive threads are consecutive in `threadIdx.x`. My memory access uses `lane = threadIdx.x`. So lane 0..31 are consecutive threads in a warp. Memory accesses at the same `i` value, address = base + lane + i*32, are perfectly coalesced. Good.
OK, let me write the kernel. Due to the complexity of float4 with scalar multiply, I'll use plain float and let nvcc handle vector loads. But I'll hint vectorization through pointer casting or pragma unroll.
Actually, I can try float4 to save instructions. Let me write helper macros:
```cpp
inline __device__ float4 operator*(const float4& a, float b) {
return make_float4(a.x * b, a.y * b, a.z * b, a.w * b);
}
inline __device__ float4 operator+(const float4& a, const float4& b) {
return make_float4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
}
inline __device__ void operator+=(float4& a, const float4& b) {
a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w;
}
inline __device__ float4 operator*(float b, const float4& a) {
return a * b;
}
```
However, operators on float4 might already be defined in some CUDA versions. Actually, no, CUDA provides float4 struct but not arithmetic operators by default (except in some headers for specific purposes). I can define them as above.
But `corr` is a scalar float. So `o4[i] = o4[i] * corr + v4 * p;` works with these operators.
For dot product with float4:
```cpp
float dot = 0;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 q4 = q_local[i];
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4.x * k4.x + q4.y * k4.y + q4.z * k4.z + q4.w * k4.w;
}
```
This is cleaner and should reduce instruction count.
For the output with float4:
```cpp
float4 o_local[8] = {};
// ...
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 v4 = reinterpret_cast<const float4*>(V + k_offset)[lane + i * 32];
float4 p4 = make_float4(p, p, p, p);
float4 corr4 = make_float4(corr, corr, corr, corr);
o_local[i] = o_local[i] * corr4 + v4 * p4;
}
```
Actually `make_float4(p,p,p,p)` for each iteration creates extra instructions. I can just use a helper function:
```cpp
__device__ float4 fma4(float4 o, float corr, float4 v, float p) {
return make_float4(o.x * corr + v.x * p, ...);
}
```
Or even simpler, just use 4 scalar FMAs:
```cpp
o_local[i].x = o_local[i].x * corr + v4.x * p;
o_local[i].y = o_local[i].y * corr + v4.y * p;
...
```
The compiler will optimize this. With 8 iterations and 4 statements each, plus unrolling, it should generate good code.
Actually, wait. With `q_local` as float4, I can't easily write `o_local[i].x = ...` if o_local is float4 array. It works, but it's verbose.
Let me try the float4 approach with manual element access, and see if it compiles.
Now, let me think about one potential issue: when `lane + i * 32` reaches 256 (for i=8), that's out of bounds for a single float4 row of 1024/4 = 256 elements. Wait:
- float4 indices for one row of 1024 floats: 0 to 255.
- lane ranges 0..31.
- i ranges 0..7.
- index = lane + i * 32.
- i=0: 0..31
- i=1: 32..63
- ...
- i=7: 224..255
- Max index = 255. Correct! And 256 * 4 = 1024 floats.
Good, no out-of-bounds.
Now, what about the alignment requirement for float4? The base pointer Q + q_offset must be 16-byte aligned. Since torch tensors are aligned to at least 256 bytes typically, and q_offset is a multiple of D*4 = 4096 bytes (since q ranges over integers), the alignment is satisfied.
Wait, for K + k_offset, k_offset is also a multiple of D*4 = 4096. So it's also aligned.
OK let me write the code now.
I'll use:
- Block: (32, 32, 1) = 1024 threads
- Grid: (32, 32, 16) = 16384 blocks
- q = blockIdx.z * 32 + threadIdx.y
Wait, grid dimensions: (32, 32, 16) means:
- blockIdx.x: 0..31 (batch)
- blockIdx.y: 0..31 (head)
- blockIdx.z: 0..15 (query group)
But there are only 32 batches and 32 heads. The total output is [B, H, N, D].
For block (b, h, g), it processes queries g*32 .. g*32+31 for batch b and head h.
Now, I should double check: the grid has 32*32*16 = 16384 blocks. Each block has 32 warps, each handling 1 query. Total queries handled = 16384 * 32 = 524288. Total queries in tensor = 32*32*512 = 524288. Perfect match.
Wait, I need to make sure the kernel handles N properly. N=512, groups=16. group=15 handles queries 480..511. All valid.
Edge case: what if N is not divisible by 32? We have `if (q >= N) return;`. In this specific problem N=512, so it's exactly divisible and the check is never true. But it's good practice.
Let me write the complete kernel code now.
Actually, one issue with using `blockIdx.z`: the maximum is quite large, but only 16. Fine.
But I realize something: in my grid=(32,32,16), the launch configuration would be:
```cpp
dim3 grid(32, 32, 16);
dim3 block(32, 32, 1);
```
Wait, block has 32*32 = 1024 threads. Grid has 32*32*16 = 16384 blocks. Total threads = 16.7M.
The kernel code:
```cpp
__global__ void fused_attention(
const float* __restrict__ Q,
const float* __restrict__ K,
const float* __restrict__ V,
float* __restrict__ O,
int B, int H, int N, int D)
{
int batch = blockIdx.x;
int head = blockIdx.y;
int q_group = blockIdx.z;
int lane = threadIdx.x;
int warp = threadIdx.y;
int q = q_group * 32 + warp;
if (q >= N) return;
int stride = H * N * D;
int head_stride = N * D;
int q_offset = batch * stride + head * head_stride + q * D;
float scale = 1.0f / sqrtf((float)D);
// Load Q
float4 q4[8];
#pragma unroll
for (int i = 0; i < 8; ++i) {
q4[i] = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
}
float o_reg[32];
#pragma unroll
for (int i = 0; i < 32; ++i) o_reg[i] = 0.0f;
float m = -1e20f;
float l = 0.0f;
for (int k = 0; k < N; ++k) {
int k_offset = batch * stride + head * head_stride + k * D;
// Compute dot product
float4 k4[8];
float dot = 0.0f;
#pragma unroll
for (int i = 0; i < 8; ++i) {
k4[i] = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4[i].x * k4[i].x;
dot += q4[i].y * k4[i].y;
dot += q4[i].z * k4[i].z;
dot += q4[i].w * k4[i].w;
}
// Warp reduce to get full dot product
#pragma unroll
for (int offset = 16; offset > 0; offset /= 2) {
dot += __shfl_xor_sync(0xFFFFFFFF, dot, offset);
}
dot *= scale;
// Online softmax
float new_m = fmaxf(m, dot);
float corr = expf(m - new_m);
float p = expf(dot - new_m);
// Load V and update output
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 v4 = reinterpret_cast<const float4*>(V + k_offset)[lane + i * 32];
o_reg[i*4+0] = o_reg[i*4+0] * corr + p * v4.x;
o_reg[i*4+1] = o_reg[i*4+1] * corr + p * v4.y;
o_reg[i*4+2] = o_reg[i*4+2] * corr + p * v4.z;
o_reg[i*4+3] = o_reg[i*4+3] * corr + p * v4.w;
}
l = l * corr + p;
m = new_m;
}
float norm = 1.0f / l;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 o4;
o4.x = o_reg[i*4+0] * norm;
o4.y = o_reg[i*4+1] * norm;
o4.z = o_reg[i*4+2] * norm;
o4.w = o_reg[i*4+3] * norm;
reinterpret_cast<float4*>(O + q_offset)[lane + i * 32] = o4;
}
}
```
Wait, I need to check: does `q4[i]` register array consume 8*4 = 32 registers? Yes. `o_reg[32]` = 32 registers. `k4[8]` = 32 registers (temporarily). But with loop unrolling, the compiler might optimize the lifetimes.
Actually, `k4` is only needed inside the K loop. But the loop has N=512 iterations, so K is loaded repeatedly. The compiler might keep `q4` in registers throughout, and `o_reg` too. `k4` might be allocated as temporaries inside the loop body.
Total register pressure: q4(32) + o_reg(32) + m(1) + l(1) + corr(1) + p(1) + dot(1) + k4(32) + temps = ~100 registers. That's too much.
Hmm, 100 registers per thread * 1024 threads = 102400 registers. But an SM only has 65536 registers. So with 1024 threads, we'd be limited to 64 registers per thread. The compiler will definitely spill.
Can we optimize registers?
Alternative: don't load all Q into registers. Load Q on the fly:
```cpp
// No q4 array
for (int k = 0; k < N; ++k) {
float dot = 0.0f;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 q4 = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4.x * k4.x + q4.y * k4.y + q4.z * k4.z + q4.w * k4.w;
}
// ...
}
```
Now q4 is a temporary float4, not an array. It doesn't consume persistent registers. Same for k4 in this loop.
For the output accumulation, we need to load V too:
```cpp
for (int k = 0; k < N; ++k) {
// compute dot...
// update o_reg
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 v4 = reinterpret_cast<const float4*>(V + k_offset)[lane + i * 32];
o_reg[i*4+0] = ...;
}
}
```
q4 is loaded twice (once in dot product, but wait - if we don't cache Q, we load it 512 times! That's bad. Q for this query has 1024 floats = 4KB. Loading it 512 times = 2MB per query. Total Q traffic = 524k * 2MB = 1TB. That's terrible.
So we DO need to keep Q in registers or shared memory.
What if we load Q into shared memory? Each block has 32 queries. Q for 32 queries = 32*1024*4 = 128KB. Shared memory can hold this, but it limits occupancy.
Alternatively, each warp loads Q into SHARED MEMORY (not registers) for its one query:
- smem_q[32][32] (1024 floats for one query).
- But shared memory is per block, so all 32 warps share it.
- Wait, we'd need 32 query worths = 128KB. That's too much.
What if a block only has 8 warps = 256 threads, processing 8 queries?
- Q smem = 8*1024*4 = 32KB.
- Can have 4 blocks per SM = 128KB. But 4*256 = 1024 threads per SM. OK.
- But grid now has more blocks.
Actually, with q4 array of 8 float4s = 32 floats = 32 registers, it's not terrible. The issue is combining q4, o_reg, and k4.
What if o_reg is stored in shared memory too?
- O accumulator for 32 queries: 32*1024*4 = 128KB. Too much for shared.
What if we process outputs in chunks?
- Compute all scores first (using cached Q), then accumulate V in chunks.
- But we said storing scores for 512 keys is 512 floats per query = 2KB per query.
- For 32 warps = 64KB. This fits in shared memory!
- Actually no, 32 queries * 512 scores = 16K floats = 64KB. That's reasonable.
- We could store scores in shared memory, then process V in chunks and update global O directly.
Algorithm:
1. Load Q into registers (8 float4s = 32 regs).
2. For each key, compute score and store in shared memory: score[q][k].
3. Wait, but we need 32*512 = 16K floats = 64KB. BlockDim = 32 warps. If each warp computes 16 scores (one per key for some keys)...
This is getting complicated again. Let me try a simpler register approach with lower block size.
What if block = 8 warps = 256 threads, processing 8 queries.
Grid: (32, 32, 64) = 65536 blocks. (8 queries per block, 512/8 = 64 groups).
Each warp handles 1 query. 8 warps = 8 queries.
Total blocks = 32*32*64 = 65536. Total threads = 65536 * 256 = 16.7M. Same total threads.
With 256 threads per block:
- Register limit = 65536 / 256 = 256 registers per thread! Much more room.
- q4[8] = 32 registers.
- o_reg[32] = 32 registers.
- k4[8] temporary = 32 registers (compiler manages lifetime).
- dot, m, l, corr, p = 5.
- Total peak ~100 registers. Under 256 limit. Plenty of room.
- Also, smem is not used, so occupancy isn't limited by shared memory.
- Concurrent blocks per SM: 65536 / 128 = 512 threads. 2048 / 512 = 4 blocks per SM.
- Actually 2048 threads max / 256 threads per block = 8 blocks per SM.
- But register limit: 8 * 256 * 100 = 204800 > 65536. So max blocks = 65536 / (256 * 100) = 2.5, i.e. 2 blocks.
- 2 blocks * 256 = 512 threads per SM. That's fine.
- Total concurrent blocks = 132 * 2 = 264.
- Waves = 65536 / 264 = 248 waves? No, total threads = 16.7M.
- Concurrent threads = 132 SMs * 512 threads = 67.5k threads.
- Waves = 16.7M / 67.5k = 248 waves.
Wait, that's a lot of waves. Each wave takes some time.
Block time: each block processes 8 queries * 512 keys. Per warp: 512 keys * 1024 FMAs = 524k FMAs. At 50% peak, maybe ~15 microseconds per warp. 8 warps execute concurrently within a block (they're in the same SM). So block time ~15 microseconds.
248 waves * 15 microseconds = 3.7ms. That's not bad!
Actually, with 512 threads per SM and 2 blocks, the 2 blocks from the same SM might execute interleaved. But the wave model still roughly applies.
For 1024 threads per block:
- Concurrent threads per SM: 65536 / (1024 * 100) = 0.64 blocks due to registers. So 1 block = 1024 threads.
- Total concurrent threads = 132 * 1024 = 135k.
- Waves = 16.7M / 135k = 124 waves.
- Block time: 32 warps * 512 keys. But they're all on the same SM. Latency hiding between warps. Time might still be ~15us (or a bit more due to less parallelism per block? No, 32 warps is plenty).
- Total = 124 * 15us = 1.86ms.
1024-thread blocks seem better for overall runtime because fewer waves.
But register pressure is the issue. Can we reduce registers?
With q4[8] and o_reg[32], that's 64 persistent registers. If we can reduce o_reg...
What if instead of 32 floats for o_reg, we process V in chunks and keep only 4 or 8 outputs at a time?
Wait, no. For online softmax, we need ALL output dimensions to be updated at once because the `corr` factor applies to all dims. We could process chunks:
```cpp
for (int d_chunk = 0; d_chunk < D; d_chunk += 32) {
float o_local[32] = {0};
m = -inf; l = 0;
for (int k = 0; k < N; ++k) {
float dot = 0;
// compute full dot product using all Q
// load V[k][d_chunk:d_chunk+32]
// update o_local
}
// write o_local / l
}
```
This requires recomputing the dot product for each d_chunk! That's 32x more compute. Not worth it.
Alternative: compute and store scores in shared memory, then do V accumulation in a separate pass.
- Load Q.
- For each key, compute score and store in smem_score[warp][k] (512 floats per warp = 2KB, times 32 warps = 64KB total).
- Then find max and compute softmax weights in a separate step.
- Then for each key, broadcast weight, load V, accumulate to o_reg[32].
- o_reg[32] still needs 32 registers. But q4 is freed after score computation.
This doesn't save registers.
Another idea: use half the registers for q4 by loading Q in two passes. For D=1024, process i=0..3 first, then i=4..7. But then we load Q twice. Q is small (4KB per query), loading twice is 8KB. From L1 cache, it's fast and doesn't increase HBM traffic. And it halves Q register usage!
```cpp
for (int k = 0; k < N; ++k) {
float dot = 0.0f;
// First half of Q
float4 q4[4];
#pragma unroll
for (int i = 0; i < 4; ++i) {
q4[i] = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
}
#pragma unroll
for (int i = 0; i < 4; ++i) {
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4[i].x * k4.x + q4[i].y * k4.y + q4[i].z * k4.z + q4[i].w * k4.w;
}
// Second half of Q
#pragma unroll
for (int i = 4; i < 8; ++i) {
q4[i-4] = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
}
#pragma unroll
for (int i = 4; i < 8; ++i) {
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4[i-4].x * k4.x + ...;
}
```
This uses only 4 float4s for Q at a time = 16 floats persistent. But `q4` is local to the loop, the compiler might keep it in registers or reload from L1. If it reloads from L1, the extra Q loads are fast (L1 cache hit).
Actually, since we're loading the same Q over and over across 512 iterations, Q will definitely be in L1 cache. Re-loading is cheap.
So we don't need to keep Q in registers at all!
```cpp
for (int k = 0; k < N; ++k) {
float dot = 0.0f;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 q4 = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4.x * k4.x + q4.y * k4.y + q4.z * k4.z + q4.w * k4.w;
}
// ...
}
```
Now `q4` is just one float4 per iteration. The compiler uses 4 registers for it, reusing across iterations. `k4` similarly.
After the dot product, we compute softmax update, then load V and update `o_reg[32]`.
`o_reg` still needs 32 registers. But 32 registers + a few temporaries is very manageable! Total ~50 registers.
Wait, this is great! Q is loaded fresh each key iteration. But since all warps in the block repeatedly access the SAME Q addresses (different queries, but each warp's Q is fixed), and L1 cache on H100 is per-SM, all threads in the block will hit in L1 for Q loads after the first miss.
Actually, each warp loads Q for its own query. Query data is 4KB. L1 cache per SM is 128KB or 256KB depending on configuration. 32 warps * 4KB = 128KB. This might just barely fit in L1, or some of it might spill to L2. But even L2 bandwidth is high.
And the dot product loads Q and K. K is 4KB per key? No, one key is 4096 bytes. For 512 keys = 2MB. K doesn't stay in L1.
So per key iteration, we load:
- Q: 4KB (for one query) from L1 or L2
- K: 4KB (for one key) from L1 or L2
- V: 4KB for accumulation
Total per block per key = 32 queries * 4KB Q + 1 key * 4KB K + 1 key * 4KB V?
Actually, all 32 warps in a block are loading Q for different queries simultaneously. So total concurrent load = 32 * 4KB = 128KB of Q. Plus 32 * 4KB K? No, all warps load the same K[k][0:1023]. So K is broadcast across all warps. The hardware fetches one 4KB cache line set for K and broadcasts it.
Wait, is K broadcast across warps? L1 cache serves requests from all warps in the SM. If 32 warps all request K[k][0:1023], and the data is in L1, it can be broadcast. But L1 might not have enough read ports to serve 32 warps simultaneously.
In practice, warps have some serialization in accessing L1. But with 32 warps on an SM, and each warp issuing loads every few cycles, it should be OK. The L1 cache has high bandwidth.
The key insight is that we only need 32 persistent registers for `o_reg`, plus a few temporaries. This allows us to use 1024 threads per block comfortably without register spilling.
Let me build this kernel:
```cpp
__global__ void fused_attention_kernel(
const float* __restrict__ Q,
const float* __restrict__ K,
const float* __restrict__ V,
float* __restrict__ O,
int B, int H, int N, int D)
{
int batch = blockIdx.x;
int head = blockIdx.y;
int q_group = blockIdx.z;
int lane = threadIdx.x;
int warp = threadIdx.y;
int q = q_group * 32 + warp;
if (q >= N) return;
int stride = H * N * D;
int head_stride = N * D;
int q_offset = batch * stride + head * head_stride + q * D;
float scale = 1.0f / sqrtf((float)D);
float o_reg[32];
#pragma unroll
for (int i = 0; i < 32; ++i) o_reg[i] = 0.0f;
float m = -1e20f;
float l = 0.0f;
for (int k = 0; k < N; ++k) {
int k_offset = batch * stride + head * head_stride + k * D;
float dot = 0.0f;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 q4 = reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32];
float4 k4 = reinterpret_cast<const float4*>(K + k_offset)[lane + i * 32];
dot += q4.x * k4.x;
dot += q4.y * k4.y;
dot += q4.z * k4.z;
dot += q4.w * k4.w;
}
#pragma unroll
for (int offset = 16; offset > 0; offset /= 2) {
dot += __shfl_xor_sync(0xFFFFFFFF, dot, offset);
}
dot *= scale;
float new_m = fmaxf(m, dot);
float corr = expf(m - new_m);
float p = expf(dot - new_m);
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 v4 = reinterpret_cast<const float4*>(V + k_offset)[lane + i * 32];
o_reg[i*4+0] = o_reg[i*4+0] * corr + p * v4.x;
o_reg[i*4+1] = o_reg[i*4+1] * corr + p * v4.y;
o_reg[i*4+2] = o_reg[i*4+2] * corr + p * v4.z;
o_reg[i*4+3] = o_reg[i*4+3] * corr + p * v4.w;
}
l = l * corr + p;
m = new_m;
}
float inv_l = 1.0f / l;
#pragma unroll
for (int i = 0; i < 8; ++i) {
float4 o4;
o4.x = o_reg[i*4+0] * inv_l;
o4.y = o_reg[i*4+1] * inv_l;
o4.z = o_reg[i*4+2] * inv_l;
o4.w = o_reg[i*4+3] * inv_l;
reinterpret_cast<float4*>(O + q_offset)[lane + i * 32] = o4;
}
}
```
Wait, I realize there's a crucial issue: Q is loaded inside the K loop. The loads `reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32]` will be executed for every key iteration (512 times). This is 512 * 4KB = 2MB per query. But Q is only 4KB. So we're reading Q 512x more than necessary.
However, with L1 caching, subsequent reads of Q will be L1 hits. The extra reads don't go to HBM, but they consume L1 bandwidth and instructions.
Can the compiler optimize this and keep Q in registers? It can't keep all of it, but with the loop unrolling and 8 float4 loads, it could keep each float4 in registers for the duration of the inner loop body, then reload in the next key iteration.
The compiler will NOT cache Q across outer loop iterations. But since Q is invariant across the K loop, shouldn't the compiler optimize this? Actually, the compiler is not allowed to move loads out of loops in general if there could be aliasing. But with `__restrict__` and const, it might optimize. However, in my experience, nvcc doesn't always move global loads out of loops.
To help the compiler, we could add `__ldg` (load global with cache) to hint that the data is readonly and cacheable. Or we could manually move Q loads before the loop.
But to manually move Q loads before the loop, we need at least 32 registers. That causes the pressure issue again.
Alternatively, use shared memory to explicitly cache Q. But then we're back to shared memory limits.
Wait, what about L1 automatic caching? On modern GPUs, global loads go through L1. When a warp loads Q in iteration k=0, the data goes into L1. In iteration k=1, the same addresses are loaded from L1. The latency is just a few cycles. The instruction overhead is still there, but it's fast.
With 512 iterations and 8 float4 loads per iteration = 4096 load instructions per warp. At 1 load per cycle per SM, this is negligible compared to the compute.
Actually, the real issue is whether the loads are redundant with respect to the instruction pipeline. Issuing redundant load instructions consumes issue slots.
A warp executing the K loop does:
- 8 loads of Q
- 8 loads of K
- 32 FMAs
- 8 loads of V
- 8*4 = 32 FMAs for V update
- shuffle reduce
- 2 expf
- Total: ~16 loads + ~64 FMAs + 2 expf + some ALU.
The loads can be overlapped with FMAs due to latency hiding. Instruction issue is probably not the bottleneck.
Actually, let me calculate more carefully. A warp issues about 1-2 instructions per cycle to the execution units. For 512 iterations, if each iteration takes ~40 cycles, total is ~20000 cycles = 10 microseconds. For 32 warps in a block, blocks take ~10 microseconds. Total waves = 124. Total time = 1.24ms.
This is really fast! Much faster than I expected. And much faster than cuBLAS approach.
But I realize my estimate might be too optimistic. Let me be more conservative. Say each key iteration takes 100 cycles due to expf latency and memory. 512 * 100 = 51200 cycles = 25.6 microseconds per warp. 32 warps with good latency hiding = 25.6 microseconds per block. 124 waves = 3.2ms.
If the reference math backend takes ~15ms, then 3.2ms is much faster. Even if it's 5-6ms, we might still win.
One concern: `expf` is expensive. We call it twice per key per warp. For N=512, that's 1024 expf calls per warp. Each expf takes ~20-30 cycles. Two expf = 40-60 cycles per iteration. At 2GHz, that's 20-30 nanoseconds per iteration. 512 * 25ns = 12.8 microseconds. For 32 warps executing concurrently with instruction-level parallelism, blocks might take ~15 microseconds. 124 waves = 1.86ms.
Actually, if `expf` is so dominant, the whole kernel might be faster than I thought. But wait - is the expf throughput really the bottleneck? H100 has special function units (SFU) for transcendental ops. Each SM might have limited SFU throughput.
SFUs on Ampere/Hopper: each SM has 4 SFUs (I think). Each SFU can compute 1 `expf` per clock. So an SM can do 4 expf per cycle. With 32 warps on an SM (1024 threads), and each warp needing 2 expf per iteration... well, warps issue expf only 512 times over the entire kernel. But the latency means the SM is not fully utilized on SFUs.
If each warp issues 512 * 2 = 1024 expf instructions, and each takes say 20 cycles to retire. A block has 32 warps, so 32768 expf instructions total. If an SM can retire 4 expf per cycle, it takes 8192 cycles = 4.1 microseconds for expfs on one SM. That's not too bad.
But we're running multiple SMs. With 132 SMs, total expf capacity = 132 * 4 = 528 per cycle. Total expfs = 16384 blocks * 32 warps * 1024 expfs = 536M expfs. 536M / 528 = 1M cycles = 0.5ms. So expf is not a major bottleneck at the chip level.
OK, I think the kernel is promising. Let me try compiling it and checking correctness.
But wait, I just thought of another issue. The loads `reinterpret_cast<const float4*>(Q + q_offset)[lane + i * 32]` where q_offset = ... + q*D. If q_offset is passed in correctly, and Q is a float pointer, then Q + q_offset is float-aligned. But for float4, we need 16-byte alignment. q_offset might not be 16-byte aligned if q*D is not divisible by 4. But D=1024, and q_offset = batch*stride + head*head_stride + q*D. Since D=1024 is divisible by 4, q*D is divisible by 4. stride and head_stride are also multiples of D and thus divisible by 4. So the alignment is actually 4-byte for float, which means for float4 we need 16-byte alignment. Since everything is divisible by 1024, and 1024 / 4 = 256, so offsets are multiples of 256 float4s, which is definitely 16-byte aligned. Good.
One more edge case: when total D is not divisible by 128 (bytes per warp load), but D=1024 is divisible by 128. Good.
Now let me think about whether I should use float4 for loads. On Hopper, 128-byte loads are very efficient. Using float4 issues 128-bit loads, and 8 float4 loads cover the full 1024 bytes. Alternatively, I could use `__ldcs` or `__ldg` for uncached/cached loads. For Q which is reused across iterations, L1 caching is ideal. The default `ld.global` goes through L1 on Hopper (unlike Ampere where it bypasses L1 by default, but `ld.global.L1::128` is available). Actually, on Ampere and later, `ld.global` can use L1. The compiler usually does the right thing.
I'll use `__ldg` (load through read-only cache) for Q, K, V. `__ldg` is a hint for frequent reuse.
For Q: since it's invarient, `__ldg` is perfect.
For K and V: they are sequential and each element is read by all warps in a block (though with 1024 threads, it's 32 warps). `__ldg` might help because the read-only cache can broadcast.
```cpp
float4 q4 = __ldg(reinterpret_cast<const float4*>(Q + q_offset) + lane + i * 32);
```
Actually, I need to cast properly.
`__ldg` takes a pointer. So:
```cpp
const float4* q_ptr = reinterpret_cast<const float4*>(Q + q_offset);
float4 q4 = __ldg(q_ptr + lane + i * 32);
```
Let me check if `__ldg` works on `float4`. Yes, it works on any type.
But wait, when using `__ldg`, the address must be 16-byte aligned. Our addresses are. Good.
One more concern: `__ldg` might bypass L1 on some architectures and use the texture cache / read-only cache. For data that's reused within a warp/thread, this might not be the best. But for data read by all warps, the L1 cache and texture cache both handle it.
Actually, for K which is read by all 32 warps in a block, any cache (L1 or read-only) is good. Q is read repeatedly by one warp. L1 is ideal because Q stays hot in L1.
I'll just use default loads for now. If I want, I can use `__ldg` for all. The compiler will generate `LDG` instructions.
Now let me write the python module and test.
Wait, one more thing. I should check if the problem's `torch.nn.functional.scaled_dot_product_attention` might have different behavior for edge cases. The input tensors are `torch.rand(batch_size, num_heads, sequence_length, embedding_dimension)`. So shape is [32, 32, 512, 1024].
In my kernel, I use the standard formula: scale = 1/sqrt(D).
PyTorch's SDPA does the same: `softmax(Q @ K^T / sqrt(D)) @ V`.
Actually, I need to confirm: PyTorch fp32 SDPA uses `1/sqrt(E)` where E is the embedding dimension per head. Here D=1024.
My kernel does the same.
Let me now implement and test this kernel. I'll use the compile_kernel and run_correctness functions.
One more optimization thought: can I handle the dot product and V update more efficiently?
For the dot product, I do 8 float4 loads and 4 multiply-adds per float4. That's 32 operations. I can combine into:
```cpp
dot += q4.x * k4.x + q4.y * k4.y + q4.z * k4.z + q4.w * k4.w;
```
This is a multiply-add chain. The compiler will generate 4 FMADs (actually FMAs).
For V update:
```cpp
o_reg[i*4+0] = fmaf(p, v4.x, o_reg[i*4+0] * corr);
```
Wait, `o * corr + p * v` = `fmaf(p, v, o * corr)`. The compiler might use FMA if it can. But first it computes `o * corr`, then `fma(p, v, o*corr)`. This uses 2 multiplies and 1 add per element. An FMA instruction does 1 multiply-add. So ideally we'd want `o * corr + p * v` to be a single FMA if `corr` is 1, but it's not. Actually, `fma(p, v, o * corr)` still requires a separate multiply for `o * corr`.
But we can compute `o_reg[i] *= corr;` first, then `o_reg[i] += p * v;` using FMA for the second step. This is actually 2 instructions vs a combined instruction that might be the same. Let's just write it clearly and let nvcc optimize.
Actually, for corr and p being scalars, if we do:
```cpp
#pragma unroll
for (int i = 0; i < 32; ++i) {
o_reg[i] *= corr;
}
#pragma unroll
for (int i = 0; i < 32; ++i) {
o_reg[i] += p * v_local[i];
}
```
This might allow better instruction scheduling because the multiply and the FMA are separated. But it requires 64 instructions. If we fuse them, it's still 64 instructions.
However, with floats, the second form uses FMA: `o += p * v` (after separate multiply). But `o *= corr` is just MUL.
The first form: `o = o * corr + p * v`. If `corr` is not 1, this is not a pure FMA. It requires two multiplications. However, on modern GPUs, there's a `DFMA` or similar? No, standard FMA is `a * b + c`. Here we have `a * b + c * d`, which requires 2 multiplies and 1 add.
So there's no way to reduce it below 2 multiplies + 1 add per output element. The separated form might be better for instruction pipelining.
Alright, let me just write the kernel. I'll keep the code as above: one statement per element.
One more thing: should I use `rpowf` or `fdiv_rn` for `1.0f / sqrtf(D)`? `sqrtf` and division are done once per kernel. Negligible.
OK, let me code this up.
Actually, I realize there might be an simpler issue: the CUDA warp shuffle with `__shfl_xor_sync(0xFFFFFFFF, dot, offset)` works across the warp. But my threads are arranged as `threadIdx.x=lane, threadIdx.y=warp`. Within a warp, `threadIdx.y` is constant and `threadIdx.x` varies. The warp consists of threads with consecutive `threadIdx.x`. So the shuffle operates correctly within each warp.
Wait, is that right? In CUDA, threads in a block are linearized as:
`linear_tid = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y`
Warps consist of 32 consecutive threads in this linearization. So threads (0..31, constant y, constant z) are in the same warp. Yes!
But the warp shuffle operates on threads in the same warp. And within a warp, the lane ID is `linear_tid % 32 = threadIdx.x % 32`. Since `threadIdx.x` ranges 0..31, lane ID is exactly `threadIdx.x`. So `__shfl_xor_sync(0xFFFFFFFF, dot, offset)` correctly exchanges between threads with different `threadIdx.x` values within the same warp.
And since `dot` is computed identically by all threads in the warp (they load different q and k but end up with different partial sums), the shuffle adds them up correctly. Perfect.
But I need to make sure that threads in different warps don't participate in each other's shuffles. The mask `0xFFFFFFFF` means all threads in the warp participate. And since shuffle is warp-local, threads in warp 0 don't see warp 1's values. Good.
Now, let me also think about whether blockDim=(32,32,1) could cause issues with the number of warps. blockDim.x=32, blockDim.y=32, blockDim.z=1. Total threads=1024. Maximum allowed. This is fine.
I'll write the code and compile.
Wait, one more concern: `reinterpret_cast<const float4*>(Q + q_offset)`. Q is `const float*`, so Q + q_offset is `const float*`. Casting to `const float4*` is fine if aligned. But C++ syntax: `reinterpret_cast<const float4*>(Q + q_offset)`.
When indexing: `[lane + i * 32]`. The base pointer is float4*, so this accesses the (lane + i*32)th float4. That's correct.
Let me compile this.
Actually, wait. `threadIdx.y` ranges from 0 to 31. And `warp` = `threadIdx.y`. So 32 warps. But CUDA limits to max 32 warps per block. 32*32=1024 threads. This is exactly at the limit. OK.
What about blockIdx.z? The grid dimension is dim3(32, 32, 16). But maximum gridDim.z is 65535. 16 is fine.
Wait, dim3 constructor: `dim3 grid(32, 32, 16)` or `dim3 grid(X, Y, Z)` where Z defaults to 1. So `dim3 grid(32, 32, 16)` is fine.
Actually, I realize that `load_inline` might have limitations on what CUDA features work, but it should compile standard CUDA code.
Let me write the full python module and compile it.