reasoning
This is a scene graph relation update operation. Let me analyze the reference:
1. `nodes` shape: (B, N, D)
2. `pair_feats` shape: (B, N, N, P)
3. `edge_mask` shape: (B, N, N)
Where B=4, N=24, D=96, P=16, num_heads=4, head_dim=24.
Operations:
1. LayerNorm on nodes: (B, N, D)
2. Q, K, V projections: Linear(D -> D) each
3. Reshape to (B, num_heads, N, head_dim)
4. Q*K^T * scale: (B, H, N, N)
5. rel_mlp on pair_feats: Linear(P->P), SiLU, Linear(P->H) -> (B, N, N, H) -> permute to (B, H, N, N)
6. Add rel_bias to logits
7. Masked fill with edge_mask
8. Softmax over last dim
9. Attn * V -> (B, H, N, head_dim) -> transpose/reshape to (B, N, D)
10. Output projection Linear(D->D)
11. Residual add
This is essentially a self-attention with pairwise edge bias. N=24 is quite small, so we have small matrices. B=4, H=4.
For an H100 with 3.35 TB/s bandwidth, the main focus should be on:
1. Fusing operations to minimize memory bandwidth
2. Using shared memory effectively for the small N x N matrices
3. Avoiding multiple kernel launches
Given N=24 is small and N^2=576, we can potentially launch one block per batch and head (or one block per batch handling all heads).
Let me think about the kernel design:
- Each thread block could handle one (batch, head) pair
- N=24, head_dim=24, so N*head_dim=576 which fits well in shared memory
- The full attention matrix (N x N = 576) also fits easily
- pair_dim=16, and rel_mlp is: Linear(16->16), SiLU, Linear(16->4)
Actually wait - let me re-read. `rel_mlp` outputs `num_heads` features. So pair_feats goes through:
1. Linear(P -> P) = Linear(16 -> 16)
2. SiLU
3. Linear(P -> num_heads) = Linear(16 -> 4)
So for each edge (i,j), we get a num_heads-dimensional vector.
I can fuse the entire operation into one kernel! Since N is small (24), I can:
1. Load nodes, do LayerNorm, QKV projections in shared memory or registers
2. Do rel_mlp on pair_feats (this is independent per edge)
3. Compute attention scores with bias
4. Softmax
5. Attend to values
6. Output projection and residual add
But wait, the weights are in nn.Linear modules, so we need to use them. In the CUDA kernel, we can access them via `.weight` and `.bias` tensors.
Let me think about the approach. I'll write a fused kernel that:
1. Takes all inputs
2. For each (b, h) block:
- Computes Q, K, V for that head by indexing into the weight matrices
- Computes relation bias via the MLP
- Computes attention
- Writes output
Actually, since N=24 is small, I could even have each thread handle multiple elements.
A simpler but still efficient approach: write a kernel that fuses all the compute for a single (batch, head). With B=4, H=4, that's 16 blocks total. Each block has N=24 threads (or more).
Let me think about memory layout and what we need to do:
For LayerNorm:
- weight shape: (D,)
- bias shape: (D,)
For Q projection:
- weight shape: (D, D)
- No bias
For K projection:
- weight shape: (D, D)
- No bias
For V projection:
- weight shape: (D, D)
- No bias
For rel_mlp[0]:
- weight shape: (P, P) = (16, 16)
- bias shape: (P,)
For rel_mlp[2]:
- weight shape: (num_heads, P) = (4, 16)
- bias shape: (num_heads,)
For out_proj:
- weight shape: (D, D)
- No bias
The challenge is that in CUDA we need to handle all these linear algebra operations. For small matrices, doing them directly with loops might be okay.
Since N=24, D=96, P=16, H=4, head_dim=24:
For each batch b and head h:
- Compute normalized nodes (LayerNorm): For each node i, compute mean and var over D=96 dims, then normalize.
Actually we can precompute LayerNorm for all nodes in a separate step or in the same kernel.
But launching a separate kernel for LayerNorm + QKV projection might be okay too, then a fused attention kernel.
Wait, the evaluation is about speedup over PyTorch reference. The reference already uses torch.matmul which is highly optimized. The key wins would come from fusing the rel_mlp and the attention loop to reduce memory bandwidth.
Actually, for this small size, PyTorch might be memory bandwidth bound. Fusing should help.
Let me think about what to fuse:
Option 1: Fuse everything into one big kernel
- Pros: minimal memory traffic
- Cons: complex kernel, might have occupancy issues
Option 2: Two kernels:
- Kernel 1: LayerNorm + QKV projection
- Kernel 2: Fused attention (rel_mlp + matmul + bias + mask + softmax + attn*V + out_proj + residual)
Actually, rel_mlp operates on pair_feats which depends on pair_dim and num_heads. If we do it in the attention kernel, each block would compute the rel_bias for its head.
Wait, a better approach: rel_mlp processes pair_feats to produce a (B, H, N, N) tensor. We could compute the first layer of rel_mlp outside the kernel (since it's shared across heads), and only do the last linear layer inside.
rel_mlp = Sequential(
Linear(P, P), # output: (B, N, N, P)
SiLU(), # output: (B, N, N, P)
Linear(P, H), # output: (B, N, N, H)
)
If we do the first two layers outside, we get (B, N, N, P). Then inside the kernel, for each head h, we do the last linear layer to get (N, N). But since P=16 and H=4, doing it all inside might be fine too.
Actually for maximum fusion, let's think about how to write a single attention kernel:
For each (b, h) pair:
1. Read all nodes for batch b into shared memory (or compute Q, K, V on the fly)
2. Compute Q[b, h, i, :] and K[b, h, j, :], V[b, h, j, :] by doing the linear projections
3. For each edge (i, j): compute rel_bias from pair_feats
4. Compute attention scores, softmax
5. Compute context
6. Do output projection and sum to output
But doing the full Q, K, V projections per thread block means each block loads the same weight matrices. With B=4, H=4, and small sizes, this might be okay.
Actually, a cleaner approach: first do fused LayerNorm+Q+K+V in a single kernel, producing Q, K, V in a compact format. Then do a second kernel for attention + output.
But that requires writing intermediate Q, K, V to global memory.
Hmm, let me think about register and shared memory usage.
For each block handling one (b, h):
- Q, K rows: each is head_dim=24 floats. We can store Q for all N=24 nodes as 24*24 = 576 floats = 2.3 KB.
- Similarly K: 2.3 KB
- V: 2.3 KB
- Attention matrix: N*N = 576 floats = 2.3 KB
- Shared memory: ~10 KB, well within limits.
Actually wait, if a block handles one (b, h), we need:
- For computing Q: we need the normalized nodes (N x D = 24 x 96 = 2304 floats = 9.2 KB)
- Or we can normalize on the fly and keep normalized nodes in shared memory
Hmm, or each thread computes one node's normalized values and stores them, then we can compute Q, K, V.
Wait, let me think about this more carefully. We have:
- 16 blocks total (B=4, H=4)
- Each block computes for one (b, h) pair.
The Q projection for head h uses columns h*head_dim to (h+1)*head_dim of the Q weight matrix.
If we store normalized nodes in shared memory (N x D = 2304 floats), then all threads in the block can compute Q, K, V for their assigned rows.
But 2304 floats * 4 bytes = 9216 bytes. With Q, K, V for the head (each N x head_dim = 576 floats = 2304 bytes), total would be ~16 KB. That's fine.
Actually, we don't need to store full Q, K, V for all nodes before computing attention. We can compute dot products on the fly too.
Let me think about a different granularity: each warp handles a subset of queries, each thread handles a query.
With N=24 and head_dim=24, for a single head:
- Q is (24, 24), K is (24, 24), V is (24, 24)
- We can have 24 threads in a block, each thread handling one query row.
But a block size of 24 might not saturate an H100. Better to use more threads.
Alternative: each block handles all heads for one batch. Then we have 4 blocks.
- Block size could be 128 threads or more.
- Each thread handles some (h, i) pairs.
Actually let me calculate: for one batch, total work is:
- Compute Q, K, V: 3 * N * D * D = 3 * 24 * 96 * 96 = 663,552 FMAs
- rel_mlp: B * N * N * (P*P + P*H) but P is small
- Attention: B * H * (N * N * head_dim + N * N) = 4 * 4 * (24*24*24 + 24*24) = ~221k + something
- Output proj: N * D * D = 24 * 96 * 96 = 221k FMAs
So total compute per batch is roughly 1M FMAs. For B=4, 4M FMAs.
At 67 TFLOPS FP32, this is ~60 microseconds if fully compute bound. But it's likely memory bound because of all the separate PyTorch ops.
With 3.35 TB/s bandwidth, reading nodes: 4*24*96*4 = 36 KB, writing output: 36 KB, plus all intermediates. The reference reads/writes a lot.
For kernel design, let me try fusing as much as possible. I'll write a kernel where:
Option A: One big kernel
For each batch b:
- Compute LayerNorm for all nodes
- Compute Q, K, V for all heads
- Compute rel_mlp output
- For each head: compute attention, context
- Compute output projection
- Add residual
But this requires a lot of shared memory or registers.
Alternatively, let me think about a tile-based approach. Since N is only 24, I can process everything in registers for small tiles.
Actually, for simplicity and correctness, let me write two kernels:
Kernel 1: Fused LayerNorm + QKV
- Input: nodes (B, N, D)
- Output: q, k, v (B, H, N, head_dim)
- Fused because LayerNorm and QKV are both memory-bound
Kernel 2: Fused Attention + Output
- Input: q, k, v (B, H, N, head_dim), pair_feats (B, N, N, P), edge_mask (B, N, N)
- Output: updated nodes (B, N, D)
- Fuses: rel_mlp(last layer), attention scores, softmax, attn@v, out_proj, residual add
Wait, actually rel_mlp includes two linear layers. If we do the first Linear+SiLU outside, it's shared. But maybe we can also compute it on the fly in the attention kernel? The first layer maps P->P=(16,16), which per edge is a matvec of size 16x16. With N*N=576 edges per batch, that's 576*16*16 = 147k FMAs per batch, or 588k total for B=4.
The second layer maps P->H=(16,4), which is 576*16*4 = 36k FMAs per batch per something.
Actually, if I compute rel_mlp on the fly inside the attention kernel, each block that handles one (b, h) needs to read pair_feats for its batch (which is N*N*P=576*16=9k floats = 36 KB) and compute all P->P and P->H. But wait, the P->P part is independent of h! So multiple blocks for the same batch would redo the same work.
Better approach: precompute the first layer (Linear + SiLU) of rel_mlp outside, getting a (B, N, N, P) tensor. Then in the attention kernel, each block only does the final Linear(P, H) for its head.
Or even better: write a dedicated kernel that computes:
1. LayerNorm + QKV
2. First layer of rel_mlp + SiLU
And another kernel for attention.
Hmm, but PyTorch's sequential operations might already do a decent job for rel_mlp since P=16 is tiny.
Let me think about the main bottleneck. In the reference:
1. `self.norm(nodes)` - one kernel
2. `q_proj`, `k_proj`, `v_proj` - three separate kernels (or fused by PyTorch if it sees the pattern)
3. `torch.matmul(q, k.transpose)` - one kernel
4. `self.rel_mlp(pair_feats)` - multiple kernels
5. `logits + rel_bias` - one kernel
6. `masked_fill` - one kernel
7. `softmax` - one kernel
8. `matmul(attn, v)` - one kernel
9. `transpose + reshape` - one kernel
10. `out_proj` - one kernel
11. `nodes + output` - one kernel
That's 12+ kernel launches! Even if each is fast, the launch overhead adds up for small tensors.
So the biggest win will be from fusing these into fewer kernels.
Let me design a two-kernel approach:
Kernel 1: `fused_norm_qkv`
- Computes LayerNorm on nodes
- Computes q, k, v projections
- Writes q, k, v to global memory
Kernel 2: `fused_attention`
- Computes attention scores: q @ k^T * scale + rel_bias
- Applies mask and softmax
- Computes attn @ v
- Applies out_proj
- Adds residual
For rel_bias, we can compute the full rel_mlp inside kernel 2... but again, the first layer is shared across heads.
Actually for simplicity and because P is small, let me compute rel_mlp inside kernel 2 but only the final layer if we precompute the intermediate. Or compute it all - with B=4 and N=24, the total work is tiny and we might not notice the redundancy.
Wait, I want to minimize global memory traffic. Let's see what we'd need to pass between kernels if we do a 2-kernel approach:
Kernel 1 output: q, k, v as (B, H, N, head_dim) = 4*4*24*24*4 = 36,864 bytes.
Kernel 2 also needs: pair_feats (B, N, N, P) and edge_mask (B, N, N) and the original nodes.
Actually, if we precompute rel_mlp first layer outside, that's another global memory read/write.
Let me think... What if I just fuse the entire thing into ONE kernel? With N=24, D=96, P=16, H=4, the shared memory requirements are manageable.
For ONE batch element in shared memory:
- Normalized nodes: N*D = 24*96 = 2304 floats = 9.2 KB
- Q, K, V for all heads: B*H*N*head_dim = 4*4*24*24 = but wait, per block
If a block handles one batch element:
- Norm nodes: 2304 floats (can be computed and stored)
- Q, K, V: each H*N*head_dim = 4*24*24 = 2304 floats, total 6912 floats = 27.6 KB
- Attention matrix per head: H*N*N = 4*576 = 2304 floats = 9.2 KB
- Context before out_proj: H*N*head_dim = 2304 floats = 9.2 KB
- For rel_mlp: pair_feats is N*N*P = 576*16 = 9216 floats = 36.8 KB, which is large!
So storing pair_feats in shared memory is too much for one block.
But we don't need to store all pair_feats. We just need rel_bias for each head, and rel_bias is H*N*N = 4*576 = 2304 floats = 9.2 KB.
So we can compute rel_bias from pair_feats on the fly! Or compute it once and store.
If we compute rel_bias per edge on the fly:
- Each edge needs: read P=16 floats, do Linear(16,16)+SiLU+Linear(16,4) to get 4 floats
- Over 576 edges, this is 576*16 = 9216 reads per batch, done once per batch.
Actually, if a block handles one batch, the block can compute rel_bias for all edges and store in shared memory. That's 9.2 KB. Then proceed with attention.
So here's my plan for a block handling one batch element (B=4 blocks total):
Shared memory needs:
1. norm_nodes: N*D = 2304 floats = 9.2 KB
2. qkv: could be separate arrays. Actually q, k, v for all heads: 3*H*N*head_dim = 3*2304 = 6912 floats = 27.6 KB
3. rel_bias: H*N*N = 2304 floats = 9.2 KB
4. attn_scores: H*N*N = 2304 floats... but we can reuse rel_bias memory if we're careful. Or compute attention per head.
5. output_ctx: H*N*head_dim = 2304 floats = 9.2 KB... but if we do head by head, we only need N*head_dim per head
6. output_before_residual: N*D = 2304 floats = 9.2 KB
Total if we do everything: ~70+ KB, which is close to the 227KB limit but okay.
But wait, we don't need all of these simultaneously! We can stream through heads.
Let me redesign: for each batch element, process heads one by one or in groups.
For a block handling one batch element:
1. Compute norm_nodes, store in SMEM (9.2 KB)
2. Compute Q, K, V for all heads (or head by head). If all heads: 27.6 KB.
3. For each head h:
a. Compute rel_bias for this head (or compute all rel_bias up front)
b. Compute attention scores for this head
c. Softmax
d. Compute context
e. Add to output accumulator
Actually, if we process head by head:
- We need norm_nodes (9.2 KB)
- Q for current head: N*head_dim = 576 floats = 2.3 KB
- K for current head: 2.3 KB
- V for current head: 2.3 KB
- rel_bias for current head: N*N = 576 floats = 2.3 KB
- attn_scores: 576 floats = 2.3 KB
- output accumulator: N*D = 2304 floats = 9.2 KB
Total: ~32 KB. Very reasonable!
Wait, but generating Q, K, V for each head means each thread does a lot of work.
Actually, maybe 4 blocks (one per batch) is too few to saturate the GPU. Let me use more blocks. I could have:
- One block per (batch, head) = 16 blocks
- Or tile the work further
With 16 blocks, each block handles all N queries for one head. That's N=24 threads minimum, but we can use more. 128 threads per block * 16 blocks = 2048 threads. An H100 has many SMs, so this might not fully saturate. But with the work per thread being significant, it should be okay.
Wait, actually 16 blocks on an H100 (which has ~132 SMs) means most SMs are idle. We need more parallelism.
Alternative: each block handles multiple queries for one (batch, head), or we have each block handle one batch with multiple heads.
Let me think about thread-level parallelism. For a single head:
- Computing Q for one query: D dot product with weight matrix column. Q projection: D * head_dim = 96 * 24 = 2304 FMAs for one query.
- Computing K: same
- Computing V: same
- Computing one attention score: dot product of head_dim=24
- Computing all scores for one query: N * head_dim = 576 FMAs
- Softmax for one query: N exponentials and reductions
- Computing context for one query: N * head_dim = 576 FMAs
- Output projection for one query: D * D = 9216 FMAs (wait, no - output is D, and context per head is head_dim*H, so it's D*D actually)
Wait, the output projection takes (B, N, D) and projects to (B, N, D). For one query: D * D = 96*96 = 9216 FMAs.
So output projection is actually the heaviest part! But in our fused kernel, if we process one head at a time, we'd need to accumulate the context for all heads and then do the output projection. Or we can do the output projection once at the end.
Let me reconsider. Maybe I should:
1. Compute Q, K, V for all heads
2. Compute attention and get context for all heads -> reshape to (B, N, D)
3. Do output projection and residual add
This is essentially what PyTorch does, but we fuse it into one kernel.
For the kernel, let me use a block size that handles multiple query nodes.
Actually, I think the cleanest approach is:
CUDA kernel that handles one batch element per block, with many threads cooperating.
Block size: maybe 256 threads (or 512).
For batch b:
1. Threads cooperatively compute norm_nodes for all N nodes by block-level reduction.
2. Threads cooperatively compute Q, K, V for all H heads.
3. Compute rel_bias for all heads. We can have each thread compute a subset of edges.
4. Then for each head, compute attention and context.
5. Finally, do output projection and residual.
With N=24 and 256 threads, threads would be underutilized for some operations but that's okay. The main thing is fusing.
Actually, a better approach: process multiple batch elements per block? No, that gets complicated.
Let me think about warp-level efficiency. A warp is 32 threads.
For computing norm_nodes: each node has D=96 elements. Two warps can handle one node (32 threads per 48 elements, or 96 threads per node, etc.). With N=24 nodes, we need 24*96/32 = 72 warps. At 256 threads/block (8 warps), we'd need ~9 blocks per batch element to parallelize this perfectly.
But the total work is so small that we don't need perfect parallelism. The main bottleneck is memory traffic from all the separate kernels.
Let me try a simpler design:
Kernel 1: `fused_norm_qkv`
- Launch: B blocks, each with H warps or threads
- Each block handles one batch.
- Within block: compute norm for all nodes, then compute QKV for all heads.
- Store Q, K, V in output tensor.
Kernel 2: `fused_attention_output`
- Launch: B * H blocks (16 blocks)
- Each block handles one (batch, head)
- Within block: load Q, K for the head.
- Compute rel_bias on-the-fly or load it.
- Compute attention scores, softmax.
- Compute context.
- Atomic-add context into an accumulator (B, N, D) in global memory.
- After all heads done, do output projection and residual.
But this requires synchronization between heads for the output projection... Unless we do the output projection inside kernel 2 with atomic adds, which might be expensive.
Alternative for kernel 2:
- Compute context for each head
- Apply output projection weight per head: wait, out_proj is (D, D), applied to concatenated heads. We can't apply it per head easily unless we split the weight.
Actually, out_proj.weight has shape (D, D). If context has shape (B, H, N, head_dim) which is (B, N, D) when flattened, then out_proj is just a matrix multiply with the weight.
So in kernel 2, if each block handles one head, it can compute partial output: for each node n, out_partial[n, :] = context[n, :] @ W_out[head_start:head_end, :]
Where W_out[head_start:head_end, :] is a (head_dim, D) slice of the output weight. Then we can atomic-add these partial outputs together.
Actually, since out_proj has no bias, the output for node n is:
output[n, d] = sum_h sum_hd context[n, h, hd] * W_out[h*head_dim + hd, d]
So each head block computes a (N, D) partial output and atomically adds to a (B, N, D) accumulator.
Hmm, atomic adds on float might have performance issues but with only 16 blocks atomically adding to (N*D = 2304) elements, contention should be minimal.
Wait, actually if each block computes the full (N, D) partial output and then atomically adds, we're doing 16 atomic operations per output element, which is fine.
Actually, we can avoid atomics by having each block write its partial output to a (B, H, N, D) tensor and then reduce + residual in a third kernel. Or have the first block of each batch do the final sum.
Let me think... with atomics, each output element gets 4 atomic adds (one per head). That's very manageable.
So kernel 2 design:
- One block per (batch, head) = B*H = 16 blocks
- Each block has say 128 threads
- Threads compute:
1. Load Q[b, h, :, :] and K[b, h, :, :] to shared memory (N * head_dim each)
2. For each query i (distributed across threads):
a. Compute attention scores with all keys: dot(Q[i], K[j]) for j=0..N-1
b. Add rel_bias[i, j]
c. Apply mask
d. Softmax
e. Compute context[i] = sum_j attn[j] * V[j]
3. After computing all contexts for this head:
context is (N, head_dim)
Compute out_partial[n, d] = sum_hd context[n, hd] * W_out[h*head_dim + hd, d]
Atomic-add out_partial to output[n, d]
4. Threads also compute rel_bias on the fly
Wait, but computing rel_bias on the fly means each block does the full rel_mlp for all N*N edges but only extracts the bias for its head. That's wasteful since the first layer is shared.
Let me think about this differently. The first layer of rel_mlp:
- Linear(16, 16) + bias
- SiLU
This maps pair_feats to an intermediate. If we precompute this in a small kernel, we get (B, N, N, 16). Then each attention block reads P=16 floats per edge and does a Linear(16, 1) for its head.
Actually, we can compute the first layer in kernel 1 or a separate small kernel.
For kernel 2, if we have precomputed intermediate of shape (B, N, N, 16):
- Each block loads the intermediate for edges (i, j) and does a matvec with W_rel2[:, h] to get the bias for head h.
Since N=24, N*N=576, and P=16, reading 16 floats per edge is 576*16*4 = 36 KB per block. That's not too bad.
But we can also compute the rel_mlp entirely on-the-fly! The weights are small:
- W_rel1: (16, 16) = 256 floats = 1 KB
- b_rel1: (16) = 64 bytes
- W_rel2: (4, 16) = 64 floats = 256 bytes
- b_rel2: (4) = 16 bytes
We can load these weights into shared memory or registers. Then each thread, when computing attention scores for edge (i, j), does:
1. Read pair_feats[b, i, j, :] (16 floats)
2. Compute intermediate = W_rel1 * pair_feats + b_rel1
3. Apply SiLU
4. bias = W_rel2[h, :] dot intermediate + b_rel2[h]
Wait, but each thread would read pair_feats 16 floats for each edge. If a block has T threads and each thread processes N/T queries, each thread accesses N pairs (one per key). Actually, if we have 24 threads each handling one query, and there are 24 keys, each thread reads 24 edges * 16 floats = 384 floats. That's 1.5 KB per thread. With 24 threads, total read is 36 KB. Not too bad for global memory.
But doing the rel_mlp computation per edge per thread adds overhead. Since the matrices are tiny, it's all registers basically.
Actually, let me compute the cost. Per edge:
- Linear1: 16 multiplies, 15 adds = 31 ops
- SiLU: a few ops
- Linear2: 16 multiplies, 15 adds = 31 ops
Total: ~70 ops per edge
For one query, 24 edges: 1680 ops
For all 24 queries: 40k ops
Q*K^T per query: 24 dots of 24 dims = 24*24*2 = 1152 ops
For all queries: 27k ops
So rel_mlp is actually comparable in cost to Q*K^T! That's why fusing it is important.
If we precompute the first layer, we still have the second layer: 16*2 = 32 ops per edge, total ~18k ops for all queries, which is significant but less.
OK let me try a balanced approach:
A single fused kernel that processes everything. Since the total work per batch is small, using fewer blocks with more work per block is fine, as long as we minimize global memory traffic.
Design: one block per batch element (4 blocks total). Each block has T threads.
Tasks per block:
1. LayerNorm: for each of N=24 nodes, compute mean/var over D=96 elements, normalize.
2. QKV projections: for each node, compute Q, K, V (each involves D*D = 9216 FMAs per node per projection). Wait, no. Q is dim->dim, but per head we want head_dim. Actually for each head, Q = node_norm * W_q[:, h*head_dim:(h+1)*head_dim]. Since there are H heads, full Q is node_norm * W_q (D*D = 9216 FMAs). Same for K and V.
3. rel_mlp: for all N*N edges, compute bias for all heads
4. Attention per head: Q[b,h] * K[b,h]^T + bias -> softmax -> * V[b,h]
5. Out projection and residual
Let me think about thread assignment within a block.
Block size: 256 threads.
For LayerNorm: 256 threads handling 24 nodes * 96 dims = 2304 elements. Each thread handles ~9 elements. We need block-level reductions for mean and variance.
For QKV: 256 threads handling 24 nodes. We can assign each thread several (node, head_dim) elements. Actually, Q has shape (H, N, head_dim) = (4, 24, 24) = 2304 elements. With 256 threads, each thread computes 9 elements of Q. Each element is a dot product of node_norm (96 dims) with a column of W_q (96 dims). So each thread does 9 * 96 = 864 FMAs for Q, similarly for K and V.
For rel_mlp: 576 edges. With 256 threads, some threads handle 2 edges, some handle 3. Each edge computation is small.
For attention: for each head, 24 queries. Each query needs dot products with 24 keys + softmax + weighted sum with 24 values.
Actually, if I structure the kernel properly, the block-level parallelism should work.
Let me write a cleaner kernel structure:
```cuda
__global__ void fused_gfvla_kernel(
const float* nodes, // (B, N, D)
const float* pair_feats, // (B, N, N, P)
const float* edge_mask, // (B, N, N)
const float* norm_w, // (D,)
const float* norm_b, // (D,)
const float* q_w, // (D, D)
const float* k_w, // (D, D)
const float* v_w, // (D, D)
const float* rel_w1, // (P, P)
const float* rel_b1, // (P,)
const float* rel_w2, // (H, P)
const float* rel_b2, // (H,)
const float* out_w, // (D, D)
float* output, // (B, N, D)
int B, int N, int D, int P, int H
) {
int b = blockIdx.x;
int tid = threadIdx.x;
int num_threads = blockDim.x;
// Shared memory layout
// norm_nodes: N * D
// q: H * N * head_dim
// k: H * N * head_dim
// v: H * N * head_dim
// rel_bias: H * N * N
// attn_scores: N * N (per head, reused)
// softmax_sum: N (per head)
// out_acc: N * D (accumulator before residual)
extern __shared__ float smem[];
float* norm_nodes = smem;
float* q = norm_nodes + N * D;
float* k = q + H * N * head_dim;
float* v = k + H * N * head_dim;
float* rel_bias = v + H * N * head_dim;
float* attn_scores = rel_bias + H * N * N;
float* out_acc = attn_scores + N * N;
// Step 1: LayerNorm
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// Compute mean and variance per node - block-level reduction
for (int n = 0; n < N; n++) {
// Each thread computes partial sum
float sum = 0.0f;
float sq_sum = 0.0f;
for (int d = tid; d < D; d += num_threads) {
float val = norm_nodes[n * D + d];
sum += val;
sq_sum += val * val;
}
// Block reduction
// ... use warp shuffle then shared memory
float mean = ...;
float var = ...;
float inv_std = rsqrtf(var + 1e-5f);
for (int d = tid; d < D; d += num_threads) {
norm_nodes[n * D + d] = (norm_nodes[n * D + d] - mean) * inv_std * norm_w[d] + norm_b[d];
}
}
__syncthreads();
// Step 2: QKV projections
// For each head, compute Q, K, V
for (int h = 0; h < H; h++) {
for (int idx = tid; idx < N * head_dim; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int col = h * head_dim + hd;
float q_val = 0.0f, k_val = 0.0f, v_val = 0.0f;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
q_val += nd * q_w[d * D + col];
k_val += nd * k_w[d * D + col];
v_val += nd * v_w[d * D + col];
}
q[h * N * head_dim + idx] = q_val;
k[h * N * head_dim + idx] = k_val;
v[h * N * head_dim + idx] = v_val;
}
}
__syncthreads();
// Step 3: rel_mlp
for (int idx = tid; idx < H * N * N; idx += num_threads) {
int h = idx / (N * N);
int n = (idx / N) % N;
int m = idx % N;
// Read pair_feats[b, n, m, :]
float interm[16]; // P=16
for (int p = 0; p < P; p++) {
interm[p] = pair_feats[b * N * N * P + n * N * P + m * P + p];
}
// Linear1
float temp[16];
for (int p = 0; p < P; p++) {
temp[p] = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
temp[p] += interm[pp] * rel_w1[pp * P + p];
}
}
// SiLU
for (int p = 0; p < P; p++) {
temp[p] = temp[p] * sigmoid(temp[p]);
}
// Linear2
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += temp[p] * rel_w2[h * P + p];
}
rel_bias[idx] = bias;
}
__syncthreads();
// Step 4: Attention per head
for (int h = 0; h < H; h++) {
// Compute attention scores
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0.0f;
for (int hd = 0; hd < head_dim; hd++) {
score += q[h * N * head_dim + i * head_dim + hd] *
k[h * N * head_dim + j * head_dim + hd];
}
score *= scale;
score += rel_bias[h * N * N + idx];
// Apply mask
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f; // or -INFINITY but careful with softmax
}
attn_scores[idx] = score;
}
__syncthreads();
// Softmax per row
for (int i = 0; i < N; i++) {
// Find max for numerical stability
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn_scores[i * N + j]);
}
// Reduce max across block
// ...
// Compute exp and sum
float exp_sum = 0.0f;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn_scores[i * N + j] - max_val);
attn_scores[i * N + j] = e;
exp_sum += e;
}
// Reduce sum
// ...
// Normalize
for (int j = tid; j < N; j += num_threads) {
attn_scores[i * N + j] /= exp_sum;
}
}
__syncthreads();
// Compute context = attn @ v
for (int idx = tid; idx < N * head_dim; idx += num_threads) {
int i = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0.0f;
for (int j = 0; j < N; j++) {
ctx += attn_scores[i * N + j] * v[h * N * head_dim + j * head_dim + hd];
}
// Add to out_acc
int out_idx = i * D + h * head_dim + hd;
if (h == 0) out_acc[out_idx] = ctx;
else out_acc[out_idx] += ctx;
}
__syncthreads();
}
// Step 5: Output projection and residual
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float out_val = 0.0f;
for (int dd = 0; dd < D; dd++) {
out_val += out_acc[n * D + dd] * out_w[dd * D + d];
}
output[b * N * D + idx] = nodes[b * N * D + idx] + out_val;
}
}
```
Wait, this is getting very complicated with the block-level reductions. And I need to be very careful about shared memory layout and sizes.
Let me compute shared memory:
- norm_nodes: 2304 floats
- q: 2304 floats
- k: 2304 floats
- v: 2304 floats
- rel_bias: 2304 floats
- attn_scores: 576 floats
- out_acc: 2304 floats
Total: 13440 floats = 53.8 KB
That's well within limits. But the kernel is complex.
Actually, some optimizations:
1. Don't need to store Q, K, V for all heads simultaneously if we process heads one by one. But we need Q and K for the attention matmul. Actually, for head-by-head, we'd recompute Q and K each time or store norm_nodes and recompute. Storing norm_nodes (2304 floats) is cheaper than Q+K+V (6912 floats).
Wait, but we'd need to recompute Q, K, V for each head, which is redundant. Better to compute once and store, or compute and use immediately.
Alternative approach: compute Q, K, V on the fly per head from norm_nodes. Since norm_nodes is in shared memory, and Q for one head is stored in registers or shared memory, this is fine.
Let me redesign with head-by-head processing to reduce peak shared memory:
```cuda
__global__ void fused_kernel(...) {
int b = blockIdx.x;
int tid = threadIdx.x;
extern __shared__ float smem[];
float* norm_nodes = smem; // N * D
float* q = norm_nodes + N * D; // N * head_dim
float* k = q + N * head_dim; // N * head_dim
float* v = k + N * head_dim; // N * head_dim
float* attn = v + N * head_dim; // N * N
float* out_acc = attn + N * N; // N * D
// Load and normalize nodes
...
// Initialize out_acc to 0
...
for (int h = 0; h < H; h++) {
// Compute Q, K, V for this head
for (int idx = tid; idx < N * head_dim; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int col = h * head_dim + hd;
float qv = 0, kv = 0, vv = 0;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[d * D + col];
kv += nd * k_w[d * D + col];
vv += nd * v_w[d * D + col];
}
q[idx] = qv;
k[idx] = kv;
v[idx] = vv;
}
__syncthreads();
// Compute attention scores + rel_bias
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0;
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// Compute rel_bias on the fly
// Read pair_feats, do rel_mlp
float interm[16];
for (int p = 0; p < P; p++) {
interm[p] = pair_feats[b*N*N*P + i*N*P + j*P + p];
}
// linear1 + silu
float tmp[16];
for (int p = 0; p < P; p++) {
tmp[p] = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
tmp[p] += interm[pp] * rel_w1[pp*P + p];
}
float s = tmp[p];
tmp[p] = s / (1.0f + expf(-s)); // SiLU = x * sigmoid(x) -- wait
// Actually SiLU(x) = x * sigmoid(x)
tmp[p] = s * (1.0f / (1.0f + expf(-s)));
}
// linear2
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h*P + p];
}
score += bias;
if (edge_mask[b*N*N + i*N + j] < 0.5f) {
score = -INFINITY; // careful with exp later
}
attn[idx] = score;
}
__syncthreads();
// Softmax
for (int i = 0; i < N; i++) {
// Block reduction for max
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i*N + j]);
}
// warp reduce...
// shared reduce...
float row_max = ...;
float exp_sum = 0;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i*N + j] - row_max);
attn[i*N + j] = e;
exp_sum += e;
}
// reduce sum...
float sum = ...;
for (int j = tid; j < N; j += num_threads) {
attn[i*N + j] /= sum;
}
}
__syncthreads();
// Context
for (int idx = tid; idx < N * head_dim; idx += num_threads) {
int i = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0;
for (int j = 0; j < N; j++) {
ctx += attn[i*N + j] * v[j * head_dim + hd];
}
// Atomic add to out_acc or accumulate
atomicAdd(&out_acc[i * D + h * head_dim + hd], ctx);
}
__syncthreads();
}
// Output projection and residual
...
}
```
This gets really complex with the block reductions for softmax, and shared memory layout needs to be precise. Let me see if I can make it simpler by using a different thread organization.
Actually, for N=24 and head_dim=24, another approach is to have each warp process one query row. A warp has 32 threads, and there are 24 keys. For each query, we have 24 threads active (or 32 with some idle), each thread handles one key's dot product.
Wait, but the dot product of Q[i] with K[j] is 24 elements, which could be parallelized across threads too.
Actually, for query i, we need:
1. Compute scores[i,j] for j=0..N-1
2. Apply rel_bias and mask
3. Softmax across j
4. Accumulate V[j] weighted by attn[i,j]
If one warp handles one query:
- There are 24 queries. With 32 threads per warp, we could have 24 threads each computing one key's score.
- But computing the dot product of Q[i] with K[j] requires a reduction across head_dim=24 elements. We could have 24 threads each computing Q[i,d]*K[j,d] for one d, then warp reduce. But j varies, so this is tricky.
Alternative: one warp handles one query. Within the warp:
- Thread 0..23 compute dot product Q[i]·K[j] for j=0..23. Each thread computes the full dot product in a loop (head_dim=24, so 24 FMAs).
- Total for all scores: 24 queries * 24 dots = 576 dots. At 24 dots per warp, we need 24 warps. Or, a warp handles multiple queries.
With 8 warps per block (256 threads), each warp can handle 1 query, and each block handles 8 queries. For N=24 queries, we need 3 blocks per batch per head? No wait, if a block has 8 warps and handles 8 queries, and we have 4 batches * 4 heads = 16 (b,h) pairs, each needing 3 blocks = 48 blocks total.
Hmm, that's more blocks and might be easier to manage than one giant block per batch.
Actually, let me try a hybrid:
- Kernel 1: LayerNorm + QKV, outputting Q, K, V to global memory. Launch B blocks, each with enough threads.
- Kernel 2: For each (b, h), use one block with 256 threads. Within the block, use cooperative groups or simple warp-level operations.
For kernel 2, each block handles one (b, h). Block computes:
1. Load Q, K, V for this head to shared memory
2. Compute rel_bias on the fly (or precomputed)
3. Attention + softmax
4. Write context to a (B, N, D) accumulator
Then a third kernel does:
5. Output projection + residual
This 3-kernel approach reduces complexity. And we avoid atomic adds in kernel 2 by writing context directly to a (B, H, N, head_dim) buffer (which is 36 KB total - tiny), then a third kernel reads this and does out_proj + residual.
Hmm, but with 3 kernels we have more launch overhead. Still much better than 12+ kernels.
Actually, if I use 2 kernels:
- Kernel 1: LayerNorm + QKV + rel_mlp + attention. For each (b, h), one block.
- But then out_proj needs all heads...
OK, let me think about this differently. I can use a 2-kernel approach that stores intermediate context in global memory:
Kernel 1: LayerNorm + QKV
- Input: nodes
- Output: qkv (B, 3, H, N, head_dim)
Kernel 2: Attention + OutProj + Residual
- For each (b, h): load q, k, v for the head, compute attention and context
- Atomic-add context into a (B, N, D) buffer
- Wait, out_proj can't be done per-head easily with atomic add because it's a dense matmul.
Wait, actually it CAN! The output projection weight W_out has shape (D, D) or (D, D). If context per head has shape (N, head_dim), then for each head we can compute:
partial_out[n, d] = sum_{hd} context[n, hd] * W_out[h*head_dim + hd, d]
Then atomicAdd partial_out to output[n, d]. Since there are H=4 heads, each output element gets 4 atomic adds. Very manageable!
And then we just add the residual.
So kernel 2 can be:
```cuda
for each (b, h) block:
load q, k, v from global memory
compute attention with on-the-fly rel_bias
compute context
for each (n, d):
partial = sum_hd context[n, hd] * W_out[h*head_dim + hd, d]
atomicAdd(&out_buffer[b, n, d], partial)
```
Then after all kernel 2 blocks finish, a tiny kernel or even PyTorch op does `nodes + out_buffer`.
Actually, I can do the residual add inside kernel 2 as well, but then each block tries to add nodes to the final output. That would need 4 atomic adds per element. Or, one block per batch could do the final add after waiting for others... which is impossible without global synchronization.
So: kernel 2 produces the output projection result in a buffer (initialized to 0), and then a third tiny kernel adds the residual.
BUT, the residual add is trivial. We can even do it with PyTorch in ModelNew.forward: `return nodes + out_buffer`.
Wait, I can also use a different trick: have the first head's block in each batch do the residual add after ensuring all other heads have finished. But CUDA doesn't provide this without atomics.
Actually, here's another idea: each block in kernel 2 writes its partial output to a (B, H, N, D) tensor. Then a separate reduction kernel sums over H and adds residual. That's clean.
(B, H, N, D) = 4*4*24*96*4 = 147 KB. Very small.
And the reduction kernel has N*D=2304 work items per batch, with B=4 batches. Could be a simple kernel.
But the bandwidth of writing 147 KB and reading it back is tiny (less than 1% of memory bandwidth).
So the final design:
Kernel 1: `fused_norm_qkv`
- B blocks, each with T threads
- Loads nodes for batch b
- Computes LayerNorm
- Computes Q, K, V projections
- Writes to q, k, v buffers
Kernel 2: `fused_attention_per_head`
- B * H blocks
- For each (b, h):
- Load q, k, v
- Compute rel_bias on the fly from pair_feats
- Compute attention, softmax
- Compute context = attn @ v
- Compute out_partial = context @ W_out_slice
- Write out_partial[b, h, n, d] for all n, d
Kernel 3: `sum_heads_and_add_residual`
- B blocks, each handling N * D elements
- Sum over heads: out[b, n, d] = nodes[b, n, d] + sum_h out_partial[b, h, n, d]
This is clean and parallelizable.
Wait, I can merge kernel 2 and 3 by just computing the full output projection and adding residual in a smarter way... Hmm, not really since the partial sums need to be reduced.
Actually, can I avoid the reduction by doing the output projection differently? What if I compute the attention for all heads in one block?
For a block handling one batch:
- Process all 4 heads
- Accumulate context in a (N, D) shared memory buffer
- After all heads, do output projection and residual add
Then I only need B=4 blocks for the attention, and each block is bigger. With 256 threads per block:
- 4 blocks * 256 threads = 1024 threads. Still a bit low for an H100.
I could have each block handle multiple batches. Or each block handle one batch with more parallelism.
Wait, let me reconsider: 4 blocks of 256 threads each... Is that really too few? The work per block is substantial:
- For each of 4 heads: load Q, K (24*24*2 = 1152 floats); compute N*N attention scores (576 dots of 24 dims = 13,824 FMAs); softmax (24 rows); context (576 dots of 24 dims); output projection (24*96 dots of 24+24+24+24=96 dims? No - for each head, partial output is (N, head_dim) -> we need to store it).
Actually, if we accumulate the full (N, D) context in shared memory, then at the end we do one output projection.
Per batch, total attention work:
- Q*K^T: 4 * 24 * 24 * 24 = 55,296 FMAs
- Softmax: minimal
- Attn*V: 4 * 24 * 24 * 24 = 55,296 FMAs
- rel_mlp: 576 * 4 * (16*16*2 + 16*4*2) ≈ 576*4*576 = 1.3M? Let me recalculate.
- Linear1: 576 edges * 16 outputs * 16 inputs = 147,456 FMAs
- SiLU: small
- Linear2: 576 edges * 4 outputs * 16 inputs = 36,864 FMAs
Total: ~184k FMAs
- Output projection: 24 * 96 * 96 = 221,184 FMAs
Total: ~516k FMAs per batch. At 67 TFLOPS, this is ~8 microseconds. But memory bandwidth and kernel launch overhead will dominate.
With 4 blocks, each doing 516k FMAs, total is 2M FMAs. At 67 TFLOPS, minimum time is ~30 microseconds. But memory bandwidth might be the bottleneck.
Actually, I want to minimize global memory traffic. Let me trace memory accesses for the one-block-per-batch design:
Per batch block:
- Read nodes: N*D = 2,304 floats (9.2 KB) - once
- Write output: N*D = 2,304 floats (9.2 KB) - once
- Read pair_feats: N*N*P = 9,216 floats (36.8 KB) - once
- Read edge_mask: N*N = 576 floats (2.3 KB) - once
- Read weights (norm, qkv, rel_mlp, out): all cached or repeatedly read
Intermediate traffic in shared memory only. So total global memory traffic per batch is ~57 KB. For B=4, 228 KB.
At 3.35 TB/s, this is ~68 nanoseconds of bandwidth time. Completely negligible!
So the bottleneck is definitely compute. But 2M FMAs at 67 TFLOPS is ~30 microseconds. However, with launch overhead and actual efficiency, it might take 50-100 microseconds.
The PyTorch reference has ~12 kernels, each reading/writing global memory. Total traffic is much higher. So we should get a nice speedup.
Let me go with the one-block-per-batch approach. It's cleaner and more fused.
For the kernel, I'll use 256 threads per block.
Actually, one concern: with only 4 blocks, we can't utilize all SMs on the H100. But since each block has tons of work and the amount of work is small overall, having more blocks wouldn't help much because we quickly run out of work.
Wait, can I process multiple queries per block in parallel with better utilization?
Another idea: launch B*H blocks = 16 blocks for attention, plus use atomics for output projection, and then a separate kernel for residual. This gives more blocks and better SM utilization.
But then each block loads pair_feats (36 KB) and edge_mask (2.3 KB), which for 16 blocks is more total memory traffic: 16 * 39 KB = 624 KB, vs 4 * 39 KB = 156 KB. However, pair_feats is the same for all heads in a batch. So blocks for the same batch read the same pair_feats! That's wasteful.
If I preprocess rel_mlp first layer first, then pair_feats intermediate is H-independent. Each block would read a much smaller tensor.
Hmm, I think I'll stick with one block per batch for maximum fusion.
Actually, let me reconsider the block-level reduction for LayerNorm and softmax. Let me write cleaner reduction code.
For LayerNorm mean and variance:
```cuda
for each node n:
float sum = 0, sq_sum = 0;
for (int d = tid; d < D; d += blockDim.x) {
float val = nodes[b*N*D + n*D + d];
smem[n*D + d] = val;
sum += val;
sq_sum += val * val;
}
// block reduction
for (int offset = blockDim.x/2; offset > 0; offset /= 2) {
sum += __shfl_down_sync(0xffffffff, sum, offset);
sq_sum += __shfl_down_sync(0xffffffff, sq_sum, offset);
}
// broadcast from lane 0
if (tid % 32 == 0) {
s_sum[tid/32] = sum;
s_sq_sum[tid/32] = sq_sum;
}
__syncthreads();
if (tid < 32) {
sum = (tid < (blockDim.x+31)/32) ? s_sum[tid] : 0;
sq_sum = (tid < (blockDim.x+31)/32) ? s_sq_sum[tid] : 0;
for (int offset = 16; offset > 0; offset /= 2) {
sum += __shfl_down_sync(0xffffffff, sum, offset);
sq_sum += __shfl_down_sync(0xffffffff, sq_sum, offset);
}
}
__syncthreads();
// Now lane 0 in warp 0 has the full sum and sq_sum
if (tid == 0) {
s_mean[n] = sum / D;
s_var[n] = sq_sum / D - s_mean[n] * s_mean[n];
}
__syncthreads();
float mean = s_mean[n];
float inv_std = rsqrtf(s_var[n] + 1e-5f);
for (int d = tid; d < D; d += blockDim.x) {
float val = smem[n*D + d];
smem[n*D + d] = (val - mean) * inv_std * norm_w[d] + norm_b[d];
}
```
Wait, this has a bug - I need to be careful with warp-level vs block-level reductions.
Actually, I can simplify using cooperative groups or just use __syncthreads-based reduction in shared memory.
For simplicity and correctness, let me just use shared memory for reductions:
```cuda
// Block-level sum reduction
__device__ float block_sum(float val, float* shared) {
int tid = threadIdx.x;
int warp_id = tid / 32;
int lane_id = tid % 32;
// Warp reduce
for (int offset = 16; offset > 0; offset /= 2)
val += __shfl_down_sync(0xffffffff, val, offset);
if (lane_id == 0)
shared[warp_id] = val;
__syncthreads();
// Reduce across warps
if (tid < 32) {
val = (tid < (blockDim.x + 31) / 32) ? shared[tid] : 0;
for (int offset = 16; offset > 0; offset /= 2)
val += __shfl_down_sync(0xffffffff, val, offset);
}
return val;
}
```
Hmm, but only lane 0 of each warp in the second reduction participates. It's tricky. Let me use a simpler shared memory approach:
```cuda
// Even simpler: all threads write to shared, then tree reduction
// But with 256 threads, this takes 8 steps
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x/2; s > 0; s >>= 1) {
if (tid < s) shared[tid] += shared[tid + s];
__syncthreads();
}
float total = shared[0];
```
This is clean and works for any block size. It's slightly slower than warp shuffle but for D=96 and N=24, who cares?
Similarly for softmax, we need max and sum per row of 24 elements. With 256 threads and block-level syncs, we can do:
```cuda
for (int i = 0; i < N; i++) {
// Find max
float max_val = -1e9f;
for (int j = tid; j < N; j += blockDim.x)
max_val = fmaxf(max_val, attn[i*N + j]);
shared[tid] = max_val;
__syncthreads();
for (int s = blockDim.x/2; s > 0; s >>= 1) {
if (tid < s) shared[tid] = fmaxf(shared[tid], shared[tid + s]);
__syncthreads();
}
float row_max = shared[0];
// Compute exp and sum
float exp_sum = 0;
for (int j = tid; j < N; j += blockDim.x) {
float e = expf(attn[i*N + j] - row_max);
attn[i*N + j] = e;
exp_sum += e;
}
shared[tid] = exp_sum;
__syncthreads();
for (int s = blockDim.x/2; s > 0; s >>= 1) {
if (tid < s) shared[tid] += shared[tid + s];
__syncthreads();
}
float row_sum = shared[0];
for (int j = tid; j < N; j += blockDim.x)
attn[i*N + j] /= row_sum;
}
```
This is clean and correct.
Now let me think about the QKV projection. With norm_nodes in shared memory and weight matrices in global memory, each thread computes some elements of Q, K, V.
Actually, I realize I should also think about weight matrix access patterns. The weight matrices are:
- q_w: (D, D) = (96, 96) = 9216 floats
When computing Q[b, h, n, hd] = sum_d norm_nodes[b, n, d] * q_w[d, h*head_dim + hd]
To maximize cache efficiency, threads that compute elements with the same d should load the corresponding norm_nodes value together.
Actually, with D=96 and head_dim=24, and H=4:
Thread 0 computes Q[b,0,0,0] = dot(norm_nodes[0,:], q_w[:,0])
Thread 1 computes Q[b,0,0,1] = dot(norm_nodes[0,:], q_w[:,1])
...
Thread 23 computes Q[b,0,0,23] = dot(norm_nodes[0,:], q_w[:,23])
Thread 24 computes Q[b,0,1,0] = dot(norm_nodes[1,:], q_w[:,0])
etc.
Doing these as nested loops is fine for small sizes. The weights will be cached in L1 since threads access nearby columns.
Similarly for rel_mlp: pair_feats is (B, N, N, P), accessed as pair_feats[b, i, j, :]. For edge (i,j), all threads computing different heads need the same pair_feats values. Might as well compute it on the fly as needed.
Wait, but if I compute rel_bias on the fly inside the attention loop, each thread computes it once per edge. With 256 threads and 576 edges, edges are distributed across threads, so each thread computes a few edges. That's fine.
OK let me write the kernel more carefully. I'll focus on correctness first, then optimize.
One issue: I'm planning to put everything in one kernel with one block per batch. With 4 blocks total, and an H100 having ~132 SMs, many SMs will be idle. But since the total work is small, more blocks wouldn't help much - we're not compute bound, we're more like "kernel fusion bound" (reducing overhead).
But wait, to improve occupancy and potentially get better scheduling, maybe I should have more, smaller blocks. For example, one block per (batch, head) = 16 blocks for the attention, plus a separate block for layer norm and qkv.
Hmm, but the layer norm and qkv for all heads are highly correlated. If I split them, I increase memory traffic.
What about one block per batch for LayerNorm+QKV, and one block per (batch, head) for attention?
- Kernel 1: 4 blocks
- Kernel 2: 16 blocks
- Kernel 3: 4 blocks for output proj + residual
This gives better parallelism for the attention part.
Let me trace memory:
Kernel 1 reads nodes (36 KB) and writes q, k, v (110 KB total)
Kernel 2 reads q, k, v (110 KB), pair_feats (147 KB), edge_mask (9 KB), and writes context_per_head (147 KB)
Kernel 3 reads context_per_head (147 KB) and nodes (36 KB), writes output (36 KB)
Total global memory: ~600 KB. Still tiny.
And kernel 2 has 16 blocks, which is much better for SM utilization.
Let me design this 3-kernel approach:
Kernel 1: `fused_norm_qkv`
```cuda
__global__ void fused_norm_qkv(
const float* nodes, const float* norm_w, const float* norm_b,
const float* q_w, const float* k_w, const float* v_w,
float* q_out, float* k_out, float* v_out,
int B, int N, int D, int H, int head_dim
) {
int b = blockIdx.x;
int tid = threadIdx.x;
extern __shared__ float smem[];
float* norm_nodes = smem;
// Load nodes
for (int idx = tid; idx < N * D; idx += blockDim.x) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// LayerNorm
for (int n = 0; n < N; n++) {
float sum = 0, sq_sum = 0;
for (int d = tid; d < D; d += blockDim.x) {
float val = norm_nodes[n * D + d];
sum += val;
sq_sum += val * val;
}
// block reduction
sum = blockReduceSum(sum);
sq_sum = blockReduceSum(sq_sum);
float mean = sum / D;
float var = sq_sum / D - mean * mean;
float inv_std = rsqrtf(var + 1e-5f);
for (int d = tid; d < D; d += blockDim.x) {
float val = norm_nodes[n * D + d];
norm_nodes[n * D + d] = (val - mean) * inv_std * norm_w[d] + norm_b[d];
}
__syncthreads();
}
// QKV projections
int total_qkv = H * N * head_dim;
for (int idx = tid; idx < total_qkv; idx += blockDim.x) {
int h = idx / (N * head_dim);
int n = (idx / head_dim) % N;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float qv = 0, kv = 0, vv = 0;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[d * D + w_col];
kv += nd * k_w[d * D + w_col];
vv += nd * v_w[d * D + w_col];
}
q_out[b * H * N * head_dim + idx] = qv;
k_out[b * H * N * head_dim + idx] = kv;
v_out[b * H * N * head_dim + idx] = vv;
}
}
```
Kernel 2: `fused_attention_per_head`
```cuda
__global__ void fused_attention_per_head(
const float* q_in, const float* k_in, const float* v_in,
const float* pair_feats, const float* edge_mask,
const float* rel_w1, const float* rel_b1,
const float* rel_w2, const float* rel_b2,
float* context_out,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x / H;
int h = blockIdx.x % H;
int tid = threadIdx.x;
extern __shared__ float smem[];
float* q = smem;
float* k = q + N * head_dim;
float* v = k + N * head_dim;
float* attn = v + N * head_dim;
// Load Q, K, V
for (int idx = tid; idx < N * head_dim; idx += blockDim.x) {
q[idx] = q_in[b * H * N * head_dim + h * N * head_dim + idx];
k[idx] = k_in[b * H * N * head_dim + h * N * head_dim + idx];
v[idx] = v_in[b * H * N * head_dim + h * N * head_dim + idx];
}
__syncthreads();
// Compute attention scores + rel_bias
for (int idx = tid; idx < N * N; idx += blockDim.x) {
int i = idx / N;
int j = idx % N;
float score = 0;
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// rel_mlp
float interm[16];
for (int p = 0; p < P; p++) {
interm[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
for (int p = 0; p < P; p++) {
tmp[p] = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
tmp[p] += interm[pp] * rel_w1[pp * P + p];
}
// SiLU
tmp[p] = tmp[p] / (1.0f + expf(-tmp[p])); // sigmoid
tmp[p] = tmp[p] * (tmp[p] / (1.0f + expf(-tmp[p]))); // wait no
// SiLU(x) = x * sigmoid(x)
float s = tmp[p];
tmp[p] = s * (1.0f / (1.0f + expf(-s)));
}
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax
for (int i = 0; i < N; i++) {
float max_val = -1e9f;
for (int j = tid; j < N; j += blockDim.x) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = blockReduceMax(max_val);
float exp_sum = 0;
for (int j = tid; j < N; j += blockDim.x) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = blockReduceSum(exp_sum);
for (int j = tid; j < N; j += blockDim.x) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Context = attn @ v
for (int idx = tid; idx < N * head_dim; idx += blockDim.x) {
int i = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0;
for (int j = 0; j < N; j++) {
ctx += attn[i * N + j] * v[j * head_dim + hd];
}
context_out[b * H * N * head_dim + h * N * head_dim + idx] = ctx;
}
}
```
Wait, I can simplify. Instead of writing context_out and then doing output projection separately, I can do output projection directly in kernel 2!
Each block computes context for one head, which is (N, head_dim). Then it computes out_partial = context @ W_out[head_start:head_end, :] and atomically adds to a global accumulator.
Kernel 2 modified:
```cuda
__global__ void fused_attention_out(
...
const float* out_w,
float* out_acc, // (B, N, D), initialized to 0
...
) {
...
// After computing context
for (int idx = tid; idx < N * D; idx += blockDim.x) {
int n = idx / D;
int d = idx % D;
float partial = 0;
for (int hd = 0; hd < head_dim; hd++) {
partial += context[n * head_dim + hd] * out_w[(h * head_dim + hd) * D + d];
}
atomicAdd(&out_acc[b * N * D + idx], partial);
}
}
```
And then kernel 3 just does `output = nodes + out_acc`.
This saves writing/reading the context buffer!
Total global memory with this approach:
Kernel 1: read nodes (36KB), write qkv (110KB)
Kernel 2: read qkv (110KB), pair_feats (147KB), edge_mask (9KB), write out_acc (36KB, with atomic adds so read-modify-write)
Kernel 3: read out_acc (36KB), nodes (36KB), write output (36KB)
Total: ~556 KB. Very small.
Wait, atomic adds on global memory are more expensive than simple writes. But with 4 atomic adds per element and good cache behavior, should be fine.
Actually, another idea: instead of atomics, pre-allocate (B, H, N, D) and have each block write its partial output there. Then do a reduction kernel.
(B, H, N, D) = 4*4*24*96*4 = 147 KB. Read it once in reduction: 147 KB. Total traffic: ~700 KB. Similar. And atomics might be okay for 16 blocks.
Hmm, but atomics on the same memory locations from 4 blocks can have contention. The L2 cache should handle it since each batch's out_acc is accessed by only 4 blocks. And within a batch, the 4 blocks write to different (n, d) elements. So minimal contention.
Actually wait, the 4 blocks for the same batch all write to the same out_acc slice. If they happen to write the same (n, d) element at the same time, there is contention. But with 2304 elements and 4 blocks, each block writes 2304 elements. The write pattern is deterministic, but if threads from different blocks execute simultaneously, they might conflict.
However, with only 4 blocks contributing to each batch's out_acc, and H100 having many SMs, these 4 blocks might execute at different times. Or if they execute concurrently, the hardware handles atomic contention reasonably.
I'll try the atomic approach. If it's slow, I'll switch to the reduction approach.
Actually, I just realized: for the 3-kernel approach, kernel 2 has B*H = 16 blocks. That's good for occupancy. Each block has its own work. Let me also make sure the block size is reasonable. I'll use 256 threads.
Let me now write the full code. But wait, I should consider whether a simpler single-kernel approach per batch might actually be better despite fewer blocks. Because all the intermediate data stays in shared memory.
Let me think about single-kernel per batch again:
- Block size: maybe 512 threads to get more parallelism within the block
- 4 blocks total
- All data in shared memory except inputs/outputs
With 512 threads and the compute being ~500k FMAs per block, time should still be dominated by arithmetic. But 4 blocks is very few. A modern GPU likes to have thousands of threads in flight.
Actually, I think I'll go with a 2-kernel approach that's a compromise:
Kernel 1: One block per (batch, head) for attention. But to avoid atomicAdd, we write context to a (B, H, N, head_dim) buffer. Block count: 16.
Kernel 2: One block per batch for out_proj + residual. Block count: 4.
And for LayerNorm+QKV... we could either:
- Include it in Kernel 1 (recompute Q, K, V from nodes on the fly, or cache norm_nodes in shared memory and compute Q, K, V per head)
- Or have a separate kernel
If we include LayerNorm+QKV in Kernel 1:
- Each block loads nodes for its batch, computes norm_nodes, then computes Q, K, V for its head, then does attention and writes context.
- The problem: 4 blocks per batch do the same LayerNorm computation. Redundant but maybe acceptable since N=24, D=96 is small.
- Wait, actually only H=4 blocks per batch do it. Redundant work: 4x.
LayerNorm cost: N*D = 2304 ops + reductions. Negligible.
QKV projection cost per head: N * D * head_dim = 24*96*24 = 55,296 FMAs. For 4 heads: 221k FMAs. If each block does its own head, total is 4x the optimal but still only 221k per block.
Actually wait, the optimal is to compute all Q, K in one go and use them. If each block computes only its head's Q, K, V, each block does N * D * head_dim FMAs = 55k FMAs. Total across 4 blocks: 221k FMAs, same as optimal. No redundancy for QKV! Because Q, K, V per head is naturally parallel.
So yes, each block can load the raw nodes, compute LayerNorm (shared across 4 blocks in same batch... slightly redundant but tiny), then compute Q, K, V for its head only.
This eliminates the need for a separate QKV buffer!
So single kernel approach:
- One block per (batch, head) = 16 blocks
- Each block:
1. Load nodes for batch b into shared memory (or registers)
2. Compute LayerNorm for all nodes (redundant across 4 blocks, but tiny)
3. Compute Q, K, V for this head
4. Compute attention with rel_bias
5. Compute context
6. Atomically add output projection partial to global output buffer
Then a tiny kernel or PyTorch op adds residual.
Let me compute the cost of LayerNorm per block:
- Read N*D = 2304 floats
- Compute mean/var (N reductions over D elements each)
- Normalize
- Write back to shared memory
With 4 blocks per batch doing this, each reading the same nodes from global memory. Total reads: 16 * 2304 = 36,864 floats = 147 KB. That's 4x the optimal 36 KB.
For a bandwidth-bound operation, this is worse. But LayerNorm is small overall.
Alternatively, I can precompute LayerNorm in a separate kernel.
Actually, I think the cleanest approach is:
1. Precompute LayerNorm in one kernel, output normalized nodes (B, N, D)
2. Each attention block reads normalized nodes and computes Q, K, V for its head
3. Attention + softmax + context
4. Atomic-add output projection partial
5. Final kernel adds residual
Step 1: 4 blocks, reads 36KB, writes 36KB
Step 2-4: 16 blocks, each reads normalized nodes (36KB), pair_feats (36KB), edge_mask (2.3KB). Total reads: 16 * 74.3 KB = 1.19 MB. That's a lot! Because each of 16 blocks reads the same normalized nodes and pair_feats for its batch.
Wait, normalized nodes is per batch. For batch b, 4 blocks read it. So total normalized nodes read: 4 batches * 4 blocks * 36 KB = 576 KB. Pair_feats: same, 576 KB. Edge_mask: 36 KB.
Total reads for step 2-4: ~1.2 MB. Plus atomics writing to out_acc.
And step 1 reads/writes 72 KB. Total ~1.3 MB.
Compare with the one-kernel-per-batch approach:
- 4 blocks total
- Each reads nodes, pair_feats, edge_mask once
- Total reads: 4 * (36 + 147 + 9) = 768 KB
- Writes output: 144 KB
- Total: ~912 KB
So one-kernel-per-batch has less memory traffic! But fewer blocks.
Hmm, for a workload this small, the memory traffic difference might not matter much. The main win is reducing kernel launch overhead. With 4 kernel launches (including the residual add) vs 12+, we should still win.
But wait, the one-kernel-per-batch has 4 blocks, and step 2 of the multi-kernel has 16 blocks. Which is better for GPU utilization?
For an H100, having only 4 active blocks might be okay if each block has a lot of warps and the scheduler can overlap memory/compute within warps. But 4 blocks on 132 SMs leaves 128 SMs idle.
However, the work per block is ~500k FMAs, which at peak would take ~8 microseconds. In practice, with memory stalls and syncs, maybe 50-100 microseconds per block. If blocks execute serially on different SMs, total time is still just the time for 4 blocks. If they execute in parallel, it's even faster. With 132 SMs and only 4 blocks, we definitely have plenty of SMs to run them concurrently.
Actually, the limiting factor is more likely the block-level syncs and the fact that we're not saturating all the compute units within an SM. A single block with 256 threads only uses part of an SM's capacity.
If I use 512 threads or 1024 threads per block, I get more parallelism within a block. But the shared memory requirements grow. Let me check:
For one block per batch with head-by-head processing:
- norm_nodes: 2304 floats = 9.2 KB
- q: 576 floats = 2.3 KB
- k: 576 floats = 2.3 KB
- v: 576 floats = 2.3 KB
- attn: 576 floats = 2.3 KB
- out_acc: 2304 floats = 9.2 KB
- plus temporary arrays for reductions
Total: ~30 KB
30 KB * 4 blocks = 120 KB shared memory used simultaneously. H100 has 228 KB per SM, so 1 block per SM fits easily. We could even have 2 blocks per SM if we wanted.
With 4 blocks spread across 4 SMs (out of 132), each SM runs one block. The other SMs are idle. But that's fine - we only have 4 units of work.
Wait, CUDA runtime might not spread them to 4 different SMs. It depends on scheduling. But either way, having 4 blocks isn't great for latency hiding because if one block stalls on a memory operation, there's no other block on the same SM to switch to.
With 16 blocks, we have more opportunities for latency hiding. The scheduler can run a different block while one waits for memory.
So maybe having 16 blocks is better even with more memory traffic?
For 16 blocks spread across, say, 16 SMs, we still have many idle SMs. But within an SM, if multiple blocks can run concurrently, we get better latency hiding.
Actually, on H100, each SM can run multiple blocks if resources allow. With shared memory of 30 KB per block and 228 KB available, an SM can run up to 7 blocks. With 16 blocks total, they can pack into 3 SMs (each running ~5 blocks) or spread across 16 SMs (1 block each).
I think having 16 blocks is generally better for latency hiding.
Let me compromise: I'll precompute LayerNorm in one small kernel (or even use PyTorch for it, since it's just one kernel and efficient). Then I'll have 16 blocks for attention, each handling one (batch, head).
Actually, I can just implement it in PyTorch inside ModelNew.forward: `h = self.norm(nodes)` and then pass `h` to the CUDA kernel. The kernel then does QKV projection, attention, and output projection.
This gives me 1 PyTorch kernel for LayerNorm, 1 CUDA kernel for everything else. Total 2 kernels. Very minimal overhead.
Wait, but the PyTorch LayerNorm kernel is highly optimized. For a small tensor (4*24*96), it's fast. Using PyTorch for it is fine.
Then the custom kernel takes normalized nodes and does:
1. Q, K, V projections per head
2. Attention with rel_bias
3. Output projection and residual add
But with 16 blocks, I still need atomics or a buffer for the output projection. Let me use a (B, N, D) output buffer initialized to 0, and atomics.
Actually, I realize I can avoid atomics by having each block write to a (B, H, N, D) buffer, then doing a reduction over H in a separate kernel. But that adds another kernel launch.
Or, I can have the reduction and residual add in the same kernel that also handles pair_feats for a different batch... no, because all blocks need to finish first.
What if I do:
Kernel: one block per (batch, head)
- Computes Q, K, V from normalized nodes
- Computes attention + context
- Writes context (N, head_dim) to global memory buffer (B, H, N, head_dim)
Then: use PyTorch to reshape and do output projection + residual:
`out = nodes + F.linear(context.reshape(B, N, D), self.out_proj.weight)`
Or even do this in a small custom kernel.
Actually, context.reshape(B, N, D) is a view (free). Then F.linear does a matmul. This is 2 PyTorch ops (view is free, linear is one kernel). So total would be:
1. PyTorch LayerNorm
2. Custom attention kernel
3. PyTorch Linear + add
3 kernel launches. Not bad.
But wait, can I avoid the PyTorch Linear by doing it in custom code? The Linear is just a matmul, which PyTorch handles well. But it requires reading the context buffer (147 KB) and out_proj weight (36 KB), writing output (36 KB). Total ~200 KB.
Actually, I can do the output projection inside the custom kernel with atomics. Let me try that.
Here's my final design:
**ModelNew.forward:**
1. `h = self.norm(nodes)` [PyTorch kernel]
2. Call custom CUDA kernel `fused_attention` with:
- h (B, N, D)
- pair_feats (B, N, N, P)
- edge_mask (B, N, N)
- q_w, k_w, v_w, rel_w1, rel_b1, rel_w2, rel_b2, out_w
- output buffer initialized to 0
3. `return nodes + output_buffer` [PyTorch kernel, or fused into the output]
Wait, for step 3, if output_buffer is just zeros with the out_proj result, then `nodes + output_buffer` is an element-wise add, which is one PyTorch kernel. But if I initialize output_buffer to `nodes`, then the custom kernel can add directly... but then multiple blocks per batch would conflict.
Actually, if I initialize out_acc[b, n, d] = nodes[b, n, d] before the kernel, then each block for batch b atomically adds its partial output. Wait no, that would add nodes multiple times!
If I initialize out_acc to 0, run the custom kernel with atomics to add partials, then add nodes at the end: that's clean.
Or, I can have a final custom kernel that adds nodes to out_acc. That's one extra kernel for 4 blocks.
Let me just do `return nodes + out_acc` in PyTorch. One extra kernel but it's fused and fast.
Actually, wait. Instead of atomics, I can have each block for head h write its partial output to out_acc[b, h, n, d] = partial. Then do a sum over h with PyTorch: `out = nodes + out_acc.sum(dim=1)`. This is two PyTorch ops (sum + add).
Hmm, or I can write a tiny final kernel that sums over H and adds residual.
Let me use the atomic approach since it avoids extra memory and kernels.
OK, final plan:
```python
class ModelNew(nn.Module):
def __init__(self, dim, pair_dim, num_heads):
super().__init__()
# ... same weights as Model
def forward(self, nodes, pair_feats, edge_mask):
h = self.norm(nodes)
B, N, D = h.shape
P = self.pair_dim
H = self.num_heads
head_dim = D // H
out_acc = torch.zeros(B, N, D, device=nodes.device, dtype=nodes.dtype)
fused_attention_kernel(
h, pair_feats, edge_mask,
self.q_proj.weight, self.k_proj.weight, self.v_proj.weight,
self.rel_mlp[0].weight, self.rel_mlp[0].bias,
self.rel_mlp[2].weight, self.rel_mlp[2].bias,
self.out_proj.weight,
out_acc,
B, N, D, P, H, head_dim, 1.0 / math.sqrt(head_dim)
)
return nodes + out_acc
```
Wait, the weight matrices from nn.Linear have shape (out_features, in_features). So:
- `q_proj.weight`: (D, D)
- `rel_mlp[0].weight`: (pair_dim, pair_dim)
- `rel_mlp[2].weight`: (num_heads, pair_dim)
- `out_proj.weight`: (D, D)
In CUDA, when I access `q_w[d * D + col]`, this assumes column-major or row-major? Let me check PyTorch's memory layout.
PyTorch tensors are row-major (last dimension contiguous). For a weight matrix of shape (out, in):
- weight[out_idx, in_idx] is at weight[out_idx * in + in_idx]
When computing `x @ W^T` for a linear layer: output[i, j] = sum_k x[i, k] * W[j, k]
In CUDA: for output element j: `sum_k x[k] * W[j * in + k]`
Wait, let me verify. For `nn.Linear(in, out)`, the weight has shape (out, in).
The computation is: `y = x @ W^T`, where x is (..., in) and W is (out, in).
So `y[j] = sum_k x[k] * W[j, k]`.
In memory: W is stored as W[j, k] at `W[j * in + k]`.
So for Q projection: Q[n, j] = sum_d h[n, d] * q_w[j, d]
In my CUDA code, for head h and head_dim hd: `w_col = h * head_dim + hd`
`q_val = sum_d h[n, d] * q_w[w_col * D + d]`
OK that checks out.
For the output projection:
`out[n, d] = sum_{h, hd} context[n, h, hd] * out_w[d, h * head_dim + hd]`
Wait, no. For `nn.Linear(D, D)` applied to x of shape (..., D):
`y[j] = sum_k x[k] * W[j, k]`
So `out[n, d] = sum_{dd} context[n, dd] * out_w[d, dd]`
Where dd ranges over D elements. And context[n, dd] for dd in [h*head_dim, (h+1)*head_dim) is from head h.
So in my CUDA code, for block (b, h):
`partial[n, d] = sum_hd context[n, hd] * out_w[d, h * head_dim + hd]`
Wait, the formula is `y[j] = sum_k x[k] * W[j, k]`. Here j = d, k = h * head_dim + hd.
So `partial[n, d] = sum_hd context[n, hd] * out_w[d * D + (h * head_dim + hd)]`
Yes, that's correct.
OK, let me write the CUDA code. But first, let me think about using vectorized loads.
Actually, for D=96, it's divisible by 8 (float4), so we can use float4 for vectorized memory access.
But for simplicity, let me first write a scalar version and see if it compiles and passes correctness.
One more thing: the edge_mask. It's a float tensor with values in {0, 1}. In the reference, it uses `masked_fill(edge_mask.unsqueeze(1) < 0.5, float("-inf"))`. For fp32, comparing with 0.5 is safe.
For the softmax, when all values are -inf (which shouldn't happen if there are valid edges, but to be safe), softmax would produce NaN. In the reference, F.softmax handles this. But in our case, if a row has no valid edges (all masked), we might get NaN. I should handle this gracefully. However, the reference would also produce NaN in this case.
Actually, let me check: if edge_mask row is all zeros, then logits are all -inf, and F.softmax produces NaN. Our kernel would also produce NaN if we use exp(-inf - (-inf)) = NaN. To avoid this, I should check if max_val is still -inf (all masked), and if so, set all attention weights to 0.
But for simplicity, let me match the reference behavior.
Hmm, let me just make sure my -inf handling works. If I set masked values to -1e9f, and the non-masked values are finite, then softmax works correctly. But -1e9f isn't exactly -inf. The reference uses -inf. Let me use -1e9f since the numerical difference won't matter for valid rows.
Actually, in float, exp(-1e9) underflows to 0, which is fine. And for all-masked rows, if all are -1e9f, softmax gives uniform distribution. But the reference gives NaN for all-masked rows. Let me match this more closely by using -INFINITY and checking.
But actually, for numerical equivalence within tolerance, exp(-1e9) is effectively 0 for float32, and the softmax on an all -1e9 row would be uniform. The reference would give NaN. This might cause issues.
To be safe, let me use `-1e9f` and accept the difference for all-masked rows (which probably don't exist in the test data, since edge_mask is all ones in get_inputs).
Actually wait, get_inputs returns `torch.ones(...)`, but the evaluation might use different inputs. Let me just be safe and handle it properly.
In the original code: `logits = logits.masked_fill(edge_mask.unsqueeze(1) < 0.5, float("-inf"))`
Then `F.softmax(logits, dim=-1)`.
In PyTorch, softmax of all -inf is NaN.
I'll use `-1e9f` for numerical stability. For typical logits that are O(1), this is effectively -inf.
Let me write the code now.
Wait, for the SiLU gradient, I need to make sure I implement it correctly:
SiLU(x) = x * sigmoid(x) = x / (1 + exp(-x))
For negative x with large magnitude, exp(-x) can overflow. So:
sigmoid(x) = 1 / (1 + exp(-x))
For x < 0, exp(-x) = exp(|x|) which can overflow.
Better: sigmoid(x) = exp(x) / (1 + exp(x)) for x < 0.
Or just use the fast approximation:
```cpp
inline float sigmoid(float x) {
return 1.0f / (1.0f + expf(-x));
}
```
For x = -88, exp(88) overflows in float. For x = 88, exp(-88) underflows to 0. So sigmoid(88) = 1, sigmoid(-88) = 0. The intermediate values might overflow.
To be safe:
```cpp
inline float sigmoid(float x) {
if (x >= 0) {
return 1.0f / (1.0f + expf(-x));
} else {
float e = expf(x);
return e / (1.0f + e);
}
}
```
For SiLU: `x * sigmoid(x)`. When x = -88, sigmoid(x) = 0, so SiLU(x) = 0, no overflow.
Let me use this stable sigmoid.
Now, for the kernel code, I'll also need block reduction functions. Let me write them carefully.
Actually, to keep things simple and correct, for reductions with small N=24, I might not need warp shuffles. A simple loop-based reduction in shared memory with __syncthreads is fine.
```cpp
__device__ float blockReduceSum(float val, float* shared) {
int tid = threadIdx.x;
int blockSize = blockDim.x;
shared[tid] = val;
__syncthreads();
for (int s = blockSize / 2; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] += shared[tid + s];
}
__syncthreads();
}
return shared[0];
}
```
Wait, this doesn't handle block sizes that aren't powers of 2. But I'll use blockDim.x = 256 which is a power of 2. For the reduction to be correct, I need all threads in the block to participate in __syncthreads.
Actually, `blockReduceSum` with a tree reduction in shared memory is fine as long as the block size is a power of 2. If blockDim.x = 256:
- s = 128: threads 0..127 add shared[128..255] to shared[0..127]
- s = 64: threads 0..63 add shared[64..127] to shared[0..63]
- ... until s = 1: thread 0 adds shared[1] to shared[0]
Then all threads can read shared[0].
For blockReduceMax:
```cpp
__device__ float blockReduceMax(float val, float* shared) {
int tid = threadIdx.x;
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x / 2; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] = fmaxf(shared[tid], shared[tid + s]);
}
__syncthreads();
}
return shared[0];
}
```
These functions need a shared memory array of size blockDim.x. Since multiple reductions happen, I'll reuse the same shared memory region.
Now let me plan the shared memory layout for the attention kernel:
For a block handling (batch b, head h):
```cpp
extern __shared__ float smem[];
// First part: normalized nodes for the batch
// Size: N * D
float* norm_nodes = smem;
// Next: Q, K, V for this head
// Size: 3 * N * head_dim
float* q = norm_nodes + N * D;
float* k = q + N * head_dim;
float* v = k + N * head_dim;
// Next: attention scores
// Size: N * N
float* attn = v + N * head_dim;
// Next: reduction scratchpad
// Size: blockDim.x
float* scratch = attn + N * N;
```
Shared memory needed: N*D + 3*N*head_dim + N*N + blockDim.x
= 2304 + 1728 + 576 + 256 = 4864 floats = 19.5 KB
Very reasonable.
Wait, I also need to store the context before output projection. Actually, I can compute context and do output projection in the same loop without storing context explicitly.
For each (n, d) in the output projection:
partial[n, d] = sum_hd context[n, hd] * out_w[d, h*head_dim + hd]
Where context[n, hd] = sum_j attn[n*N + j] * v[j*head_dim + hd]
So:
partial[n, d] = sum_hd [sum_j attn[n*N + j] * v[j*head_dim + hd]] * out_w[d, h*head_dim + hd]
= sum_j attn[n*N + j] * [sum_hd v[j*head_dim + hd] * out_w[d, h*head_dim + hd]]
We can precompute `v_proj[j, d] = sum_hd v[j*head_dim + hd] * out_w[d, h*head_dim + hd]`
Then `partial[n, d] = sum_j attn[n*N + j] * v_proj[j, d]`
This changes the computation order. For fixed n and d, we need v_proj[j, d] for all j. We can compute v_proj[j, d] on the fly or precompute it.
Hmm, actually this might be more efficient:
1. Compute `v_proj[j, d] = sum_hd v[j, hd] * out_w[d, h*head_dim + hd]` for all j, d
This is a small matmul: (N, head_dim) @ (head_dim, D) = (N, D)
2. Then `partial[n, d] = sum_j attn[n, j] * v_proj[j, d]`
This is another small matmul: (N, N) @ (N, D) = (N, D)
Total work: N*head_dim*D + N*N*D = 24*24*96 + 576*96 = 55k + 55k = 110k FMAs
Direct computation:
For each n: context[n, hd] = sum_j attn[n,j] * v[j, hd] -> N*N*head_dim = 576*24 = 13.8k FMAs
Then partial[n, d] = sum_hd context[n, hd] * out_w[d, ...] -> N*D*head_dim = 24*96*24 = 55.3k FMAs
Total per head: 69k FMAs
So the direct approach is 69k vs 110k. Direct is better.
And we don't need extra shared memory for v_proj.
So I'll compute context per (n, hd), then compute partial output per (n, d) using the context.
Since context for one head is N*head_dim = 576 floats, I can store it in shared memory or compute on the fly.
Option A: Store context in shared memory
- context_smem: N * head_dim
- Then threads read context_smem to compute partial output
Option B: Each thread computes context for some (n, hd), and then we have a second phase where threads compute partial output.
Either way works. Let me use Option A for simplicity.
Shared memory addition:
`float* context = attn + N * N;` (reuse after softmax, or use separate)
Actually, attn is still needed for the context computation. After context is computed, attn can be reused. But context has the same size as q/k/v.
Let me use: after computing attention and context, I can reuse q's shared memory to store the context, since q is no longer needed.
Wait, q and k are needed for the attention scores. v is needed for context. After context is computed, none are needed.
So shared memory layout can be:
- norm_nodes: N*D
- q: N*head_dim
- k: N*head_dim
- v: N*head_dim
- attn: N*N
- context: N*head_dim (reuse q after computing attn scores)
Actually, when computing context, I still need v. But q and k are no longer needed. So I can reuse q's memory for context.
Layout:
- norm_nodes: N*D
- q: N*head_dim
- k: N*head_dim
- v: N*head_dim
- attn: N*N
- scratch: blockDim.x
After attention scores are computed and softmax is done:
- q and k are no longer needed
- Compute context into q (reuse q memory)
- Then compute output projection from context and atomically add
Shared memory total: 2304 + 576 + 576 + 576 + 576 + 256 = 4864 floats = 19.5 KB
Same as before.
Wait, I'm not reusing q for context. Let me be explicit:
```cpp
float* context = q; // Reuse q memory after softmax
```
For computing context:
```cpp
for (int idx = tid; idx < N * head_dim; idx += blockDim.x) {
int n = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0;
for (int j = 0; j < N; j++) {
ctx += attn[n * N + j] * v[j * head_dim + hd];
}
context[idx] = ctx;
}
__syncthreads();
```
Then for output projection:
```cpp
for (int idx = tid; idx < N * D; idx += blockDim.x) {
int n = idx / D;
int d = idx % D;
float partial = 0;
for (int hd = 0; hd < head_dim; hd++) {
partial += context[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
atomicAdd(&out_acc[b * N * D + idx], partial);
}
```
This looks correct.
Now, I need to be careful: when I write `float* context = q;`, the compiler might get confused. In CUDA, I can just use `q` directly as context after the softmax phase, since I know the memory is no longer used for Q. But to be clear, I'll just refer to `q` as the context buffer.
Actually, to be really safe and clear, let me allocate separate pointers:
```cpp
float* norm_nodes = smem;
float* q = norm_nodes + N * D;
float* k = q + N * head_dim;
float* v = k + N * head_dim;
float* attn = v + N * head_dim;
float* scratch = attn + N * N;
```
And I'll just reuse `q` for context in the code, with a comment.
OK, let me now draft the full kernel code. I'll also need to extract weight tensors from the PyTorch modules.
Actually, wait. I realize I need to handle the rel_mlp weights carefully. In PyTorch:
```python
self.rel_mlp = nn.Sequential(
nn.Linear(pair_dim, pair_dim),
nn.SiLU(),
nn.Linear(pair_dim, num_heads),
)
```
The weights are:
- `self.rel_mlp[0].weight`: shape (pair_dim, pair_dim)
- `self.rel_mlp[0].bias`: shape (pair_dim,)
- `self.rel_mlp[2].weight`: shape (num_heads, pair_dim)
- `self.rel_mlp[2].bias`: shape (num_heads,)
I need to check the exact shapes and access patterns.
For the first linear: `intermediate = pair_feats @ W1^T + b1`
Where W1 has shape (pair_dim, pair_dim). In CUDA: `intermediate[p] = b1[p] + sum_pp pair_feats[pp] * W1[p, pp]`
Wait, `nn.Linear(in, out)` computes `y = x @ W^T + b`, where W has shape (out, in).
So `y[p] = b[p] + sum_pp x[pp] * W[p, pp]`
In memory: `W[p, pp]` is at `W[p * in + pp]`.
For the second linear: `bias[h] = b2[h] + sum_p intermediate[p] * W2[h, p]`
Where W2 has shape (num_heads, pair_dim).
In the CUDA code, I compute:
```cpp
// Linear 1
temp[p] = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
temp[p] += pair_feats[pp] * rel_w1[p * P + pp];
}
```
Wait, should it be `rel_w1[p * P + pp]` or `rel_w1[pp * P + p]`?
Since W1 has shape (P, P), `W1[p, pp]` is at `W1[p * P + pp]`. Yes.
For the second linear:
```cpp
bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += temp[p] * rel_w2[h * P + p];
}
```
And W2 has shape (H, P), so `W2[h, p]` is at `W2[h * P + p]`. Yes.
OK, I think the indexing is correct.
Now let me write the full Python file with the kernel.
Actually, one more thing: I should use `const float*` for all inputs and `float*` for outputs. And make sure to pass the right tensor pointers.
Let me draft the kernel code:
```cpp
#include <torch/extension.h>
#include <cuda_runtime.h>
#include <math.h>
inline float sigmoid_fast(float x) {
if (x >= 0) {
return 1.0f / (1.0f + expf(-x));
} else {
float e = expf(x);
return e / (1.0f + e);
}
}
__device__ float block_reduce_sum(float val, float* shared) {
int tid = threadIdx.x;
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] += shared[tid + s];
}
__syncthreads();
}
return shared[0];
}
__device__ float block_reduce_max(float val, float* shared) {
int tid = threadIdx.x;
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] = fmaxf(shared[tid], shared[tid + s]);
}
__syncthreads();
}
return shared[0];
}
__global__ void fused_attention_kernel(
const float* __restrict__ nodes, // (B, N, D), pre-normalized
const float* __restrict__ pair_feats, // (B, N, N, P)
const float* __restrict__ edge_mask, // (B, N, N)
const float* __restrict__ q_w, // (D, D)
const float* __restrict__ k_w, // (D, D)
const float* __restrict__ v_w, // (D, D)
const float* __restrict__ rel_w1, // (P, P)
const float* __restrict__ rel_b1, // (P,)
const float* __restrict__ rel_w2, // (H, P)
const float* __restrict__ rel_b2, // (H,)
const float* __restrict__ out_w, // (D, D)
float* __restrict__ out_acc, // (B, N, D)
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x / H;
int h = blockIdx.x % H;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* norm_nodes = smem; // N * D
float* q = norm_nodes + N * D; // N * head_dim
float* k = q + N * head_dim; // N * head_dim
float* v = k + N * head_dim; // N * head_dim
float* attn = v + N * head_dim; // N * N
float* scratch = attn + N * N; // blockDim.x
// Load normalized nodes for this batch
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// Compute Q, K, V for this head
int qkv_size = N * head_dim;
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float qv = 0.0f, kv = 0.0f, vv = 0.0f;
#pragma unroll 4
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[w_col * D + d]; // Wait, is this right?
kv += nd * k_w[w_col * D + d];
vv += nd * v_w[w_col * D + d];
}
q[idx] = qv;
k[idx] = kv;
v[idx] = vv;
}
__syncthreads();
```
Wait! I made an indexing error. For `nn.Linear(in, out)` with weight shape (out, in), the element is `W[out_idx, in_idx]` at `W[out_idx * in + in_idx]`.
For `q_proj(h)`: `q[n, w_col] = sum_d h[n, d] * q_w[w_col, d]`
In memory: `q_w[w_col * D + d]`. Yes, this is correct!
Wait, but PyTorch uses row-major: `tensor[i, j]` is at `tensor[i * stride0 + j * stride1]`. For a contiguous 2D tensor, `stride1 = 1`, so `tensor[i, j]` is at `tensor[i * dim1 + j]`. So `q_w[w_col, d]` is at `q_w[w_col * D + d]`. Yes, that's correct.
OK, continuing:
```cpp
// Compute attention scores with rel_bias
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0.0f;
#pragma unroll 4
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// Compute rel_mlp for edge (i, j)
float pf[16];
#pragma unroll
for (int p = 0; p < P; p++) {
pf[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
#pragma unroll
for (int p = 0; p < P; p++) {
float val = rel_b1[p];
#pragma unroll
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
// SiLU
val = val * sigmoid_fast(val);
tmp[p] = val;
}
float bias = rel_b2[h];
#pragma unroll
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
// Apply mask
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax per row
for (int i = 0; i < N; i++) {
// Find max
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = block_reduce_max(max_val, scratch);
// Compute exp and sum
float exp_sum = 0.0f;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
// Normalize
for (int j = tid; j < N; j += num_threads) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Compute context = attn @ v, store in q (reuse)
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0.0f;
#pragma unroll 4
for (int j = 0; j < N; j++) {
ctx += attn[n * N + j] * v[j * head_dim + hd];
}
q[idx] = ctx; // Reuse q as context buffer
}
__syncthreads();
// Output projection and atomic add
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float partial = 0.0f;
#pragma unroll 4
for (int hd = 0; hd < head_dim; hd++) {
partial += q[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
atomicAdd(&out_acc[b * N * D + idx], partial);
}
}
```
Wait, there's a potential issue with the block_reduce functions. The `block_reduce_max` and `block_reduce_sum` functions use `__syncthreads()`. When called inside the softmax loop over N rows, for each row we call block_reduce_max and then block_reduce_sum. These functions internally call `__syncthreads()`, and the loop itself has a `__syncthreads()` at the end.
In CUDA, all threads in a block must execute `__syncthreads()`. The functions do this correctly. But the functions also have `return shared[0]` at the end, and `__syncthreads()` is called before reading `shared[0]`. Actually, `return shared[0]` after the final `__syncthreads()` is fine because all threads have synchronized and written to shared[0].
Wait, but in `block_reduce_sum`, after the final `__syncthreads()`, all threads return `shared[0]`. That's correct because shared[0] has the final result.
However, there's a subtle issue: in the loop `for (int i = 0; i < N; i++)`, all threads must participate in the loop. If some threads skip the loop body because of a condition near the beginning of the loop, they might not reach the `__syncthreads()` at the end. In my code, `for (int j = tid; j < N; j += num_threads)` - all threads participate. And both block_reduce_* functions have __syncthreads(). So all threads should reach all sync points.
But wait, the `block_reduce_max` and `block_reduce_sum` functions have a __syncthreads() inside them. When called inside the loop, the loop doesn't have a __syncthreads() before them. But since the previous iteration's __syncthreads() and the current iteration's __syncthreads() in the reduce functions handle it, it should be OK.
Actually, there's a problem. After computing the exp values:
```cpp
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
```
Here, `block_reduce_sum` uses `scratch[tid] = val` at the beginning. But `scratch` was just used by `block_reduce_max` in the same iteration. After `block_reduce_max` finishes, all threads have synchronized. Then threads compute their partial exp_sum, and `block_reduce_sum` starts. This should be fine.
Similarly for the normalize loop after block_reduce_sum - there's a __syncthreads() at the end of the iteration.
OK, let me make sure the shared memory for `scratch` is sized correctly. `scratch[threadIdx.x]`. Since threadIdx.x can be up to 255 (for blockDim.x=256), we need 256 floats. But actually blockDim.x might be 256. Let me compute based on actual launch.
Actually, to be safe and support variable block sizes, I can just use `blockDim.x` in the kernel. But I'll also hardcode to 256 for efficiency.
Now, one potential issue: the `scratch` array is used for both `block_reduce_max` and `block_reduce_sum` in the same iteration. After `block_reduce_max`, `scratch[0]` contains the max. Then `block_reduce_sum` overwrites `scratch[tid]` with the partial sums. This is fine.
Another issue: the Softmax loop has `__syncthreads()` at the very end, after the normalize loop. Then the next iteration of `i` starts. At the beginning of the next iteration, threads compute `max_val` from `attn[i*N + j]`, which was written in previous iterations. Since there was a `__syncthreads()` at the end of the previous iteration, all threads have finished writing. Good.
Now, what about `out_w` indexing? Let me double check.
`out_w` has shape (D, D). For `nn.Linear(D, D)`:
`out[n, d] = sum_{dd} context[n, dd] * out_w[d, dd]`
In my code:
```cpp
partial += q[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
```
Wait, `q` here is the context buffer, storing context for head h. So `q[n * head_dim + hd]` is `context[n, hd]`. And `out_w[d * D + h * head_dim + hd]` is `out_w[d, h * head_dim + hd]`. But the linear formula is `y[d] = sum_k x[k] * W[d, k]`.
So if x has length D, then `y[d] = sum_{k=0}^{D-1} x[k] * W[d, k]`.
For partial computation with only head h: `partial[d] = sum_{hd} context[hd] * W[d, h*head_dim + hd]`.
In my code: `out_w[d * D + h * head_dim + hd]`. Since `out_w[d, k]` is at `out_w[d * D + k]`, this is correct.
One thing: I need to make sure the indexing in the Q, K, V projections is right.
For Q projection: `q[n, w_col] = sum_d h[n, d] * q_w[w_col, d]`
In my code: `qv += nd * q_w[w_col * D + d]`
Wait, should it be `q_w[d * D + w_col]` or `q_w[w_col * D + d]`?
The weight `q_w` has shape (D, D). In PyTorch row-major: `q_w[out, in]` is at `q_w[out * in_dim + in]`.
So `q_w[w_col, d]` is at `q_w[w_col * D + d]`. Yes.
And `nn.Linear(in, out)` computes: `output[b, out] = sum_in input[b, in] * weight[out, in]`.
So `q[n, w_col] = sum_d h[n, d] * q_w[w_col, d]`. Correct.
OK, indexing is correct.
One more thing: how do I get `nodes` after LayerNorm? In ModelNew.forward:
```python
h = self.norm(nodes)
```
Then I pass `h` to the CUDA kernel. The kernel expects pre-normalized nodes.
So the kernel signature uses `nodes` but it should be the normalized nodes. I'll clarify in comments.
Let me now figure out what `dim`, `pair_dim`, `num_heads` correspond to:
- D = dim = 96
- P = pair_dim = 16
- H = num_heads = 4
- head_dim = dim / num_heads = 24
- N = num_objects = 24 (but dynamic at runtime from tensor shape)
- B = batch_size = 4 (but dynamic)
In the kernel, I'll extract B, N, D, P, H from the tensors.
Actually, I'll pass them as arguments since they're known. But also get N from nodes.shape[1].
Now, let me write the full Python code.
Wait, I need to be careful about SiLU. In PyTorch, `nn.SiLU()` is `x * sigmoid(x)`. My implementation uses `val = val * sigmoid_fast(val)`. Correct.
One more thing: for the `rel_mlp`, the first `nn.Linear(pair_dim, pair_dim)` has `bias=True` by default. The second `nn.Linear(pair_dim, num_heads)` also has `bias=True` by default. So I need to pass both biases.
In my kernel code, I do:
```cpp
float val = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
```
But wait, `pf` is loaded from `pair_feats`, which has shape (B, N, N, P). So `pf[pp]` is correct.
I also need to make sure the `pair_feats` indexing is right:
```cpp
pair_feats[b * N * N * P + i * N * P + j * P + p]
```
For a (B, N, N, P) tensor in row-major:
- stride for b: N * N * P
- stride for i: N * P
- stride for j: P
- stride for p: 1
So `element[b, i, j, p]` is at `b * N*N*P + i * N*P + j * P + p`. Yes, correct.
Similarly for edge_mask:
```cpp
edge_mask[b * N * N + i * N + j]
```
Correct.
Now, for the launch configuration:
```python
BLOCK_SIZE = 256
grid = (B * H,)
block = (BLOCK_SIZE,)
shared_mem = (N * D + 3 * N * head_dim + N * N + BLOCK_SIZE) * 4 # bytes
```
Wait, if I launch with `fused_attention_kernel<<<grid, block, shared_mem>>>`, I need to make sure the shared memory size is correct.
With B * H = 16 blocks, and block size 256:
- N * D = 2304
- N * head_dim = 576
- N * N = 576
- blockDim.x = 256
Total floats: 2304 + 3*576 + 576 + 256 = 4864
Bytes: 4864 * 4 = 19456
Actually let me recalculate:
norm_nodes: N*D = 2304
q: N*head_dim = 576
k: 576
v: 576
attn: N*N = 576
scratch: 256
Total: 2304 + 576 + 576 + 576 + 576 + 256 = 4864
Bytes: 19456
Wait, I counted q, k, v as 576 each = 1728. So 2304 + 1728 + 576 + 256 = 4864. Yes.
19456 bytes = ~19 KB. Well within limits.
Now, let me also consider using `__restrict__` and `__launch_bounds__` for better compiler optimization. But it's probably not necessary for this small kernel.
One more optimization: since P=16 is very small, the compiler can fully unroll the loops if I use `#pragma unroll`. I'll add `#pragma unroll` for the p and pp loops since P is not a template parameter but a runtime value... Actually, `#pragma unroll` with a number might help. But for variable P, the compiler won't know the exact number at compile time unless I pass it as a template parameter.
Actually, I could make P a template parameter, but it's overkill. For the test, P is always 16. The compiler might still unroll small loops heuristically.
For N and head_dim, they're also constants in practice (24), but I'll treat them as runtime values.
Actually, wait. I should make sure the kernel can handle arbitrary N, D, P, H, not just the test values. But for speed, hardcoding or templating the loop bounds would help. I'll keep it general but add `__restrict__` and hope the compiler optimizes well.
Alternatively, I can use if-statements with constants known at module load time, and pass them as template arguments via C++ function overloading. But that's getting complex.
Let me keep it simple first. I'll add some `#pragma unroll` directives manually where appropriate.
Hmm, `__syncthreads()` inside a loop can be expensive if the loop iterations have little work. For N=24 and block_size=256, the softmax loop has 24 iterations, each with a few syncs. That's fine.
Let me actually think about whether I can optimize the softmax further. With N=24, it's tiny. No real optimization needed.
Now, I'll write the full Python module.
I also realize I need to handle the residual add. I'll create an `out_acc` tensor initialized to zeros, run the kernel, then return `nodes + out_acc`.
One potential issue: numerical precision. The reference uses PyTorch's exact operations. My kernel uses float32 arithmetic, which should match well. But the order of operations differs, which can cause slight differences. The tolerance is atol=rtol=1e-4, which is relaxed.
However, for softmax, the reference uses `F.softmax` which has careful numerical stability. My implementation subtracts the max before exponentiating, which should match.
For the SiLU, `sigmoid_fast` might differ slightly from PyTorch's implementation. But within tolerance should be fine.
Let me also think about the output projection atomic adds. When multiple blocks write to the same `out_acc[b * N * D + idx]`, the order of atomic adds affects the final sum. But floating-point addition is commutative within rounding tolerance, so the final sum should be the same regardless of atomic add order.
Wait, floating-point addition is NOT associative. The order does matter for the exact result. With 4 atomic adds to the same element, the result depends on the order. Since the blocks run concurrently and the order is non-deterministic, there could be slight numerical differences.
However, the tolerance is 1e-4, so small rounding differences from atomic add order should be fine.
But to be safe, I could avoid atomics by having each block write its partial output to a separate buffer, and then do a reduction. Let me think...
With 16 blocks, 4 blocks per batch, writing to a (B, H, N, D) buffer. The buffer size is 4*4*24*96*4 = 147 KB. Then a reduction kernel sums over H.
Reduction kernel can be simple:
```cpp
__global__ void reduce_and_add_residual(
const float* partials, // (B, H, N, D)
const float* nodes,
float* output,
int B, int N, int D, int H
) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int total = B * N * D;
if (idx >= total) return;
float sum = 0.0f;
for (int h = 0; h < H; h++) {
sum += partials[idx * H + h]; // Wait, need correct indexing
}
output[idx] = nodes[idx] + sum;
}
```
Actually, the layout matters. If partials is (B, H, N, D), then:
`partials[b, h, n, d]` is at `((b * H + h) * N + n) * D + d`
The reduction over h for fixed (b, n, d):
```cpp
int b = idx / (N * D);
int rem = idx % (N * D);
int n = rem / D;
int d = rem % D;
float sum = 0;
for (int h = 0; h < H; h++) {
sum += partials[((b * H + h) * N + n) * D + d];
}
output[idx] = nodes[idx] + sum;
```
This reduction is very bandwidth-light (reads 147 KB, writes 36 KB) and easy to write.
But for numerical reproducibility, this deterministic reduction is better than atomics.
OK, I'll use the deterministic reduction approach to avoid atomic ordering issues.
Actually, with only 147 KB for the partials buffer, the memory cost is negligible.
Let me update the design:
Kernel 1: `fused_attention`
- One block per (batch, head)
- Computes attention + context + output projection partial
- Writes to `partials[b, h, n, d]`
Kernel 2: `reduce_residual`
- One thread per (b, n, d) = B*N*D = 9216 threads
- Sums over H and adds residual
This is deterministic and clean.
Actually, I can even make the reduction kernel have one block per batch, with threads handling multiple (n, d) elements.
With 9216 elements and block size 256: 9216 / 256 = 36 blocks.
Or I can launch B*N blocks of D threads each. Many options.
Let me just use a simple 1D grid:
```cpp
int total_threads = B * N * D;
int blocks = (total_threads + 255) / 256;
reduce_residual<<<blocks, 256>>>(partials, nodes, nodes_out, B, N, D, H);
```
Actually wait, the output should be `nodes + out_proj(context)`. The reference returns `nodes + self.out_proj(ctx)` where `ctx` is the attended context. So yes, residual add.
In my kernel, `partial[b, h, n, d]` is the contribution of head h to output element (b, n, d). The reduction sums over h.
OK, let me write the complete code.
Oh wait, I should double check: `ctx = torch.matmul(attn, v).transpose(1, 2).reshape(B, N, D)` has shape (B, N, D). Then `self.out_proj(ctx)` maps (B, N, D) -> (B, N, D).
In my kernel, for head h, I compute context_h of shape (N, head_dim). Then I compute partial_out_h = context_h @ W_out[head_start:head_end, :]^T?
Wait, let me re-check the linear algebra.
`out_proj` is `nn.Linear(D, D)` with weight shape (D, D).
`out = ctx @ out_proj.weight.T`
For a single output channel d:
`out[b, n, d] = sum_{k=0}^{D-1} ctx[b, n, k] * out_w[d, k]`
Now ctx[b, n, k] for k in [h*head_dim, (h+1)*head_dim) is from head h.
So for the partial from head h:
`partial[b, h, n, d] = sum_{hd=0}^{head_dim-1} context_h[n, hd] * out_w[d, h*head_dim + hd]`
In my CUDA code, I'm computing:
```cpp
partial += q[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
```
Where `q` stores `context_h`. This is `context_h[n, hd] * out_w[d, h*head_dim + hd]`. Correct.
After summing over hd, we get `partial[b, h, n, d]`. Correct.
Then the reduction kernel sums over h:
`sum_h partial[b, h, n, d] = sum_h sum_hd context_h[n, hd] * out_w[d, h*head_dim + hd]`
`= sum_{k=0}^{D-1} ctx[b, n, k] * out_w[d, k]`
`= out[b, n, d]`
Then `output[b, n, d] = nodes[b, n, d] + out[b, n, d]`. Correct.
OK, everything checks out.
One optimization: for the reduction kernel, instead of storing to partials and reading back, could I just use the custom kernel to compute all heads and then do out_proj? The problem is that a block only handles one head, and out_proj needs all heads. I could have a block handle all heads for one batch, but then I'd have fewer blocks (4 vs 16).
I'll stick with 16 blocks for attention and a separate reduction kernel.
Now let me code this up.
Wait, one more thing: the LayerNorm. In the reference, it happens before QKV. I can just use PyTorch's LayerNorm: `h = self.norm(nodes)`. This is one kernel launch. Then pass `h` to my custom kernel.
But actually, I can also fuse LayerNorm into the custom kernel. Each of the 16 blocks computes LayerNorm for its batch. It's redundant (4 times per batch), but the cost is small.
If I don't fuse it, I have one extra PyTorch kernel. But PyTorch's LayerNorm kernel is already very efficient. Let me just use PyTorch for LayerNorm to keep things simple.
So `ModelNew.forward`:
1. `h = self.norm(nodes)` (PyTorch kernel)
2. `partials = torch.zeros(B, H, N, D, device=nodes.device, dtype=nodes.dtype)`
3. Launch `fused_attention_kernel<<<B*H, 256, shared_mem>>>` with h, pair_feats, edge_mask, weights
4. Launch `reduce_residual_kernel<<<...>>>` with partials, nodes, output
5. Return output
Step 2 and 4 can be fused into one kernel that initializes the output and reduces, but for simplicity I'll keep them separate.
Actually, I don't even need a separate reduce kernel if I do the reduction in PyTorch:
```python
out = nodes + partials.sum(dim=1) # Sum over heads
```
`partials.sum(dim=1)` is a reduction kernel, and `nodes + ...` is an element-wise add kernel. Two PyTorch kernels. Total: LayerNorm + custom attention + sum + add = 4 kernels. Still very few.
Actually, I can even do: `out = nodes + partials.sum(dim=1)` where partials is computed by the custom kernel.
But PyTorch's sum might not be fully optimized for this small size. A custom reduction kernel should be better. But let's see - for a tensor of shape (4, 4, 24, 96), summing dim=1 gives (4, 24, 96). The reduction is over 4 elements. PyTorch will launch a kernel. For 9216 output elements, each summing 4 inputs, it's a trivial kernel.
I'll write a custom tiny kernel to be safe and fast.
Let me draft the complete Python file:
```python
import math
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.cpp_extension import load_inline
CUDA_CODE = r'''
#include <torch/extension.h>
#include <cuda_runtime.h>
inline __device__ float sigmoid_fast(float x) {
if (x >= 0.0f) {
return 1.0f / (1.0f + expf(-x));
} else {
float e = expf(x);
return e / (1.0f + e);
}
}
__device__ float block_reduce_sum(float val, float* shared) {
int tid = threadIdx.x;
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] += shared[tid + s];
}
__syncthreads();
}
return shared[0];
}
__device__ float block_reduce_max(float val, float* shared) {
int tid = threadIdx.x;
shared[tid] = val;
__syncthreads();
for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
if (tid < s) {
shared[tid] = fmaxf(shared[tid], shared[tid + s]);
}
__syncthreads();
}
return shared[0];
}
__global__ void fused_attention_kernel(
const float* __restrict__ nodes,
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
const float* __restrict__ q_w,
const float* __restrict__ k_w,
const float* __restrict__ v_w,
const float* __restrict__ rel_w1,
const float* __restrict__ rel_b1,
const float* __restrict__ rel_w2,
const float* __restrict__ rel_b2,
const float* __restrict__ out_w,
float* __restrict__ partials,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x / H;
int h = blockIdx.x % H;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* norm_nodes = smem;
float* q = norm_nodes + N * D;
float* k = q + N * head_dim;
float* v = k + N * head_dim;
float* attn = v + N * head_dim;
float* scratch = attn + N * N;
// Load normalized nodes
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// Compute Q, K, V for this head
int qkv_size = N * head_dim;
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float qv = 0.0f, kv = 0.0f, vv = 0.0f;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[w_col * D + d];
kv += nd * k_w[w_col * D + d];
vv += nd * v_w[w_col * D + d];
}
q[idx] = qv;
k[idx] = kv;
v[idx] = vv;
}
__syncthreads();
// Compute attention scores + rel_bias
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0.0f;
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// rel_mlp on pair_feats[b, i, j, :]
float pf[16];
for (int p = 0; p < P; p++) {
pf[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
for (int p = 0; p < P; p++) {
float val = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
val = val * sigmoid_fast(val);
tmp[p] = val;
}
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax per row
for (int i = 0; i < N; i++) {
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = block_reduce_max(max_val, scratch);
float exp_sum = 0.0f;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
for (int j = tid; j < N; j += num_threads) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Context = attn @ v, reuse q as context buffer
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
float ctx = 0.0f;
for (int j = 0; j < N; j++) {
ctx += attn[n * N + j] * v[j * head_dim + hd];
}
q[idx] = ctx;
}
__syncthreads();
// Output projection partial for this head
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float partial = 0.0f;
for (int hd = 0; hd < head_dim; hd++) {
partial += q[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
int partial_idx = ((b * H + h) * N + n) * D + d;
partials[partial_idx] = partial;
}
}
__global__ void reduce_and_residual_kernel(
const float* __restrict__ partials,
const float* __restrict__ nodes,
float* __restrict__ output,
int B, int N, int D, int H
) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int total = B * N * D;
if (idx >= total) return;
int b = idx / (N * D);
int rem = idx % (N * D);
int n = rem / D;
int d = rem % D;
float sum = 0.0f;
for (int h = 0; h < H; h++) {
sum += partials[((b * H + h) * N + n) * D + d];
}
output[idx] = nodes[idx] + sum;
}
'''
class ModelNew(nn.Module):
def __init__(self, dim: int, pair_dim: int, num_heads: int):
super().__init__()
assert dim % num_heads == 0
self.dim = dim
self.pair_dim = pair_dim
self.num_heads = num_heads
self.head_dim = dim // num_heads
self.scale = 1.0 / math.sqrt(self.head_dim)
self.norm = nn.LayerNorm(dim)
self.q_proj = nn.Linear(dim, dim, bias=False)
self.k_proj = nn.Linear(dim, dim, bias=False)
self.v_proj = nn.Linear(dim, dim, bias=False)
self.rel_mlp = nn.Sequential(
nn.Linear(pair_dim, pair_dim),
nn.SiLU(),
nn.Linear(pair_dim, num_heads),
)
self.out_proj = nn.Linear(dim, dim, bias=False)
# Compile CUDA extension
self._cuda_ext = load_inline(
name="gfvla_fused",
cpp_sources="",
cuda_sources=CUDA_CODE,
functions=["fused_attention_kernel", "reduce_and_residual_kernel"],
extra_cuda_cflags=["-O3", "--use_fast_math"],
verbose=False,
)
def forward(self, nodes, pair_feats, edge_mask):
B, N, D = nodes.shape
P = self.pair_dim
H = self.num_heads
head_dim = self.head_dim
h = self.norm(nodes)
# Allocate partial output buffer
partials = torch.zeros(B, H, N, D, device=nodes.device, dtype=nodes.dtype)
BLOCK_SIZE = 256
grid1 = (B * H,)
block1 = (BLOCK_SIZE,)
shared_mem = (N * D + 3 * N * head_dim + N * N + BLOCK_SIZE) * 4
self._cuda_ext.fused_attention_kernel(
h, pair_feats, edge_mask,
self.q_proj.weight, self.k_proj.weight, self.v_proj.weight,
self.rel_mlp[0].weight, self.rel_mlp[0].bias,
self.rel_mlp[2].weight, self.rel_mlp[2].bias,
self.out_proj.weight,
partials,
B, N, D, P, H, head_dim, self.scale,
grid=grid1, block=block1, shared_mem=shared_mem
)
output = torch.empty_like(nodes)
total_elems = B * N * D
grid2 = ((total_elems + BLOCK_SIZE - 1) // BLOCK_SIZE,)
block2 = (BLOCK_SIZE,)
self._cuda_ext.reduce_and_residual_kernel(
partials, nodes, output,
B, N, D, H,
grid=grid2, block=block2
)
return output
```
Wait, I need to check the weight tensor shapes and ensure they're contiguous.
PyTorch nn.Linear weights are stored as (out_features, in_features). When I pass `self.q_proj.weight` to CUDA, it needs to be contiguous. PyTorch weights are contiguous by default.
But let me verify the exact shapes:
- `self.q_proj.weight`: (D, D) = (96, 96)
- `self.rel_mlp[0].weight`: (pair_dim, pair_dim) = (16, 16)
- `self.rel_mlp[0].bias`: (16,)
- `self.rel_mlp[2].weight`: (num_heads, pair_dim) = (4, 16)
- `self.rel_mlp[2].bias`: (4,)
This looks correct.
One issue: `load_inline` with empty `cpp_sources` might have issues. The documentation says `cpp_sources` can be empty. Let me check.
Actually, the standard usage is:
```python
load_inline(
name="...",
cpp_sources=cpp_wrapper,
cuda_sources=cuda_code,
functions=[...],
)
```
But since we don't need a C++ wrapper, we can set `cpp_sources=""`. I think this works.
Also, I need to make sure the `load_inline` call happens once, ideally in `__init__`. Since `ModelNew` is instantiated once, this is fine.
One more potential issue: `load_inline` might not compile correctly with `__restrict__` in kernel arguments. Actually, in CUDA C++, `__restrict__` is valid on pointer arguments.
Let me also think about whether I need `extern "C"` for the kernel functions... No, the pybind11 wrapper generated by `load_inline` will handle name mangling. The functions just need to be callable.
Actually, wait. If I pass `cpp_sources=""`, will `load_inline` still generate bindings? The documentation says it auto-generates a pybind11 module from `functions=[...]`. So it should work.
Let me test this by compiling it first.
But before compiling, I want to double-check the indexing for `rel_mlp` weights.
In the kernel:
```cpp
val += pf[pp] * rel_w1[p * P + pp];
```
`rel_w1` comes from `self.rel_mlp[0].weight`, which has shape (16, 16).
In PyTorch: `weight[p, pp]` is at `weight[p, pp]`.
In CUDA: `rel_w1[p * P + pp]` corresponds to `weight[p, pp]`. Correct.
Similarly:
```cpp
bias += tmp[p] * rel_w2[h * P + p];
```
`rel_w2` comes from `self.rel_mlp[2].weight`, shape (4, 16).
`rel_w2[h * P + p]` is `weight[h, p]`. Correct.
OK, I think the code is correct. Let me also consider whether I should vectorize the pair_feats reads. With P=16, I could use float4 to read 4 floats at a time. But the gains might be minimal for such a small kernel.
Let me now write the code and compile it.
Wait, one more thing. In the kernel, when I write:
```cpp
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
```
With N=24, D=96, N*D=2304. With 256 threads, this does 9 iterations per thread for the first 2304 % 256 =... 2304 / 256 = 9 exactly! So all 256 threads execute exactly 9 iterations. Perfect load balancing.
Similarly for qkv_size = 576, with 256 threads: 256 threads do 2 iterations, 64 threads do 3. Not perfect but fine.
For N*N = 576: same.
For N*D = 2304 in the output projection: same as norm_nodes, 9 iterations each. Perfect.
The softmax loop over N=24 rows does 24 iterations, each with small work. Might have some overhead, but 24 is tiny.
Actually, in the softmax loop, for each row I call block_reduce_max and block_reduce_sum. With a block size of 256 and N=24, most threads have only a partial role in the reductions. This is somewhat inefficient. Let me think about whether I can optimize this.
For a block handling one (b, h), the attention matrix is N x N = 24 x 24 = 576 elements.
If I use fewer threads per block, say 64 or 128, the reductions would be more utilized.
With 64 threads:
- N*D load: 2304 / 64 = 36 iterations per thread
- qkv: 576 / 64 = 9 iterations
- attn: 576 / 64 = 9 iterations
- softmax: 64 threads for 24 elements per row, most threads participate
- context: 576 / 64 = 9 iterations
- output: 2304 / 64 = 36 iterations
With 64 threads, reductions are better utilized, and shared memory scratch is smaller.
Total threads: 16 blocks * 64 threads = 1024 threads. Still not fully utilizing the GPU, but better.
With 128 threads:
- Load: 18 iterations
- Reductions: 128 threads for 24 elements, still somewhat wasteful but better than 256
- Total threads: 2048
I think 128 or 64 threads per block might be better. Let me try 128.
But actually, for the H100, warp size is 32. With 64 threads, we have 2 warps. Block-level sync between 2 warps is fast. With 128 threads, 4 warps.
Actually, for this kernel, most of the time is spent in the QKV projection and output projection loops, which are compute-heavy. The softmax is a small part. So having more threads to parallelize the compute-heavy parts is probably better.
With 256 threads:
- QKV: each thread does ~2.25 dot products (576 elements / 256 threads)
- Wait, each thread doesn't do a dot product; each thread computes one element of Q/K/V.
- 576 Q elements / 256 threads = ~2.25 elements per thread. Each element is a dot product of 96.
- Total FMAs per thread for QKV: ~2.25 * 3 * 96 = ~648 FMAs
That's actually not a lot! The kernel might be execution-latency dominated rather than throughput dominated.
Hmm, with such small work per thread, having fewer threads with more work might be better. Or having each thread compute multiple elements.
Actually, maybe I should restructure the kernel to have each thread compute an entire row or multiple elements. But that would make the softmax harder since it requires inter-thread communication.
Alternatively, I can process the attention in "tiles" where each thread handles one query and loops over keys and values.
Actually, let me reconsider: with 16 blocks of 256 threads, each thread has:
- Load: 9 loads
- QKV: ~2.25 dot products of length 96
- Attention: ~2.25 dots of length 24
- rel_mlp: ~2.25 MLPs of size 16->16->1
- Softmax: participates in 24 row reductions (mostly idle in the reduction tree)
- Context: ~2.25 dots of length 24
- Output: ~9 dots of length 24
Total work per thread is actually moderate. The main inefficiency is in softmax where many threads are idle or underutilized.
To optimize softmax, I could use a smaller block size. Let me try 64 threads.
With 64 threads and 16 blocks = 1024 total threads:
- Load: 36 loads
- QKV: 9 dots of length 96 = 864 FMAs
- Attention: 9 dots of length 24 = 216 FMAs
- rel_mlp: 9 MLPs
- Softmax: 64 threads for 24 elements. In block_reduce, 32 threads do useful work for the first step (s=32). 24 threads participate in the actual loop.
- Context: 9 dots of length 24 = 216 FMAs
- Output: 36 dots of length 24 = 864 FMAs
With 64 threads, work per thread is higher, and reductions are better utilized. Total threads is 1024, which is a decent number.
Let me try block_dim = 64. But I want to make sure occupancy is good. With shared memory = 19.5 KB per block, and 228 KB per SM, each SM can hold up to 11 blocks. With 1024 threads total and 64 threads per block, we have 16 blocks. On 132 SMs, each SM runs at most 1 block (since 16 < 132), but could run up to 11. Actually, the runtime will schedule 16 blocks across available SMs.
Wait, 16 blocks on 132 SMs means about 1 block per 8 SMs. Many SMs are idle.
If I had more blocks, they'd be distributed better. But we only have 16 (B*H) independent work units.
Unless I tile further... For example, each block could handle a TILE of queries instead of all queries. But since N=24, tilting doesn't help much.
Alternatively, I could process multiple (b, h) per block. But then shared memory needs to be larger.
What if I process one batch per block (4 blocks total), handling all heads sequentially? That's what I considered earlier. With 4 blocks, even fewer.
Or what if I process one batch per block but with multiple threadblocks per batch? For example, 4 blocks per batch, each handling H/4 = 1 head? That's back to 16 blocks.
OK, 16 blocks is a natural partition based on the problem structure.
One trick: I can use 4 warps per block but launch with 128 threads? Or 32 threads per warp, 2 warps = 64 threads.
Actually, let me think about another angle. What if I fuse everything into one kernel with 4 blocks (one per batch), and each block processes all heads with 256 threads?
With 4 blocks and 256 threads, total threads = 1024. Same as 16 blocks * 64 threads.
But the 4-block approach has less memory traffic. And within each block, threads can work on multiple things.
For 4 blocks, each block handles one batch:
- 256 threads work on N*D=2304 elements for loading and output
- For QKV of all heads: 2304 elements per projection, 3 projections = 6912 elements. With 256 threads, ~27 elements per thread. Each element is a dot product of 96. 27*96 = 2592 FMAs per thread just for QKV.
- For attention: 4 heads * 576 = 2304 attention scores. 2304 / 256 = 9 per thread.
- For softmax: 4 heads * 24 rows = 96 row softmaxs. With 256 threads, some threads help multiple rows or are idle.
- For context: 4 heads * 576 = 2304 elements. ~9 per thread.
- For output: 2304 elements. ~9 per thread, each dot product of 96. 864 FMAs.
With 4 blocks * 256 threads, the work per thread is higher, especially for QKV and output projection. This seems better for utilization.
But the shared memory needs to be larger:
norm_nodes: 2304
q: 2304 (for all heads)
k: 2304
v: 2304
attn: 576 (one head at a time)
context_all_heads: 2304 (or reuse q)
out_acc: 2304
scratch: 256
Total: 2304*5 + 576 + 256 = 12496 floats = ~50 KB
50 KB shared memory per block. With 228 KB per SM, each SM can hold up to 4 blocks. But we only have 4 blocks total. They'll be distributed to 4 different SMs.
Actually, what if I DO use one block per batch? Let me rewrite the kernel for this approach.
Fused kernel per batch:
1. Load nodes
2. LayerNorm (already done by PyTorch, so load normalized nodes)
3. Compute Q, K, V for all heads
4. For each head:
a. Compute attention scores with rel_bias
b. Softmax
c. Compute context
d. Accumulate into out_acc
5. Output projection (from out_acc) + residual
Wait, but if I compute Q, K, V for all heads at once, I need more shared memory. Let me use an approach where I compute Q, K, V for each head from norm_nodes on the fly.
Actually, a better approach:
1. norm_nodes in shared memory
2. For each head h:
a. Compute Q_h, K_h, V_h in shared memory (each N*head_dim)
b. Compute attention + softmax
c. Compute context in Q_h's space (reuse)
d. Add context to out_acc
3. After all heads: output projection from out_acc + residual
Shared memory:
- norm_nodes: N*D = 2304
- q, k, v: 3*N*head_dim = 1728
- attn: N*N = 576
- out_acc: N*D = 2304
- scratch: 256
Total: 7168 floats = 28.7 KB
With 4 blocks, this is fine.
Processing per head sequentially means we have a loop of H=4 iterations. Each iteration involves computing Q, K, V for one head. If we compute Q, K, V for all heads at once, we'd use more SMEM but avoid recomputation.
Option A: compute Q, K, V for all heads at once, store in SMEM (6912 floats = 27.6 KB)
Then for each head, use the precomputed Q, K, V.
Total SMEM: 2304 + 6912 + 576 + 2304 + 256 = 12352 floats = 49.4 KB
Option B: compute Q, K, V per head on demand from norm_nodes. No extra storage for all QKVs.
Total SMEM: 2304 + 1728 + 576 + 2304 + 256 = 7168 floats = 28.7 KB
Option A has larger SMEM but avoids recomputing QKV projections 4 times. Option B saves SMEM but recomputes QKV.
Since QKV for one head is N*D*head_dim = 24*96*24 = 55,296 FMAs. For 4 heads: 221k FMAs. Recomputing 4 times means 4x55k = 221k total, same as computing once. Wait no - computing once for all heads is 221k FMAs. Computing per head 4 times: each head is 55k, total 4*55k = 221k. Same amount of arithmetic!
But computing per head means each head block computes its own. Computing all at once in one block means the block does 221k FMAs for QKV. There's no redundancy since the work is partitioned differently.
Actually wait, I was confused. In the one-block-per-batch approach, if I compute all QKVs at once, the block does 221k FMAs for QKV. If I compute per head in a loop, the block still does 221k FMAs for QKV (4 * 55k). Same arithmetic.
So Option B is just Option A with a different loop structure. The total work is the same.
Let me go with Option B (one block per batch, compute QKV per head) since it uses less shared memory.
Actually, there's another consideration: with 16 blocks (one per head), the 221k QKV FMAs are distributed across 16 blocks. With 4 blocks, they're concentrated in 4 blocks. Since the total arithmetic is the same and we're likely latency-limited, having more blocks with less work each might be better.
Hmm, but the 16-block approach also requires reading `nodes` 16 times (once per block, 4 times per batch). With 4 blocks, `nodes` is read 4 times (once per batch). So 4 blocks have less memory traffic for the input.
Let me actually compute memory traffic for both:
16-block approach (one per head):
- norm_nodes read: 16 * 2304 floats = 921.6 KB (since each of 16 blocks reads its batch's nodes)
- pair_feats read: 16 * 9216 floats = 3.69 MB
- edge_mask read: 16 * 576 floats = 230.4 KB
- partials written: 16 * 2304 floats = 921.6 KB
Total traffic: ~5.76 MB
4-block approach (one per batch):
- norm_nodes read: 4 * 2304 = 36.8 KB
- pair_feats read: 4 * 9216 = 147 KB
- edge_mask read: 4 * 576 = 9.2 KB
- output written: 4 * 2304 = 36.8 KB
Total traffic: ~230 KB
The 4-block approach has 25x less memory traffic! This is huge.
But wait, in the 4-block approach, we compute everything in shared memory and only write the final output. In the 16-block approach, we write partials and then read them back.
So the 4-block approach is much more bandwidth-efficient.
However, with only 4 blocks, we might not saturate all SMs. But the work per block is substantial:
- QKV: 221k FMAs
- Attention: 110k FMAs (for all heads)
- rel_mlp: ~184k FMAs (for all heads, since each block does it once)
- Out projection: 221k FMAs
Total: ~736k FMAs per block.
At 67 TFLOPS = 67e12 FLOPS/s, and each FMA is 2 FLOPs, the peak time for 736k FMAs is 736k / (67e12 / 2) = ~22 nanoseconds. That's impossibly fast - meaning we can't even measure it at peak. In reality, with memory stalls and instruction overhead, it will take microseconds.
But the key point: 4 blocks means only 4 SMs can run this kernel at a time. With an H100 having ~132 SMs, most are idle.
But wait, each SM can run multiple thread blocks concurrently if resources permit. With shared memory = 29 KB per block and 228 KB per SM, each SM can run up to 7 blocks. With 4 blocks total and block size 256 (8192 threads per SM max, but 256 is small), an SM could run all 4 blocks if they were assigned to it. But CUDA runtime distributes blocks to SMs. With 4 blocks, they might be placed on 4 different SMs.
However, the kernel is very compute-heavy relative to its thread count, so even with 4 SMs active, if each SM is fully utilized, the time might still be dominated by the compute.
Let me think about whether launching more blocks would help. Since we only have 4 batches, we can't launch more blocks for the main computation unless we tile further.
One tiling strategy: process half the queries in one block and half in another, for the same (batch, head). But then we'd need atomic adds for the output projection. Or we could have multiple blocks cooperate within a batch.
Actually, for the 4-block approach, we already have all queries in one block. There's no natural tiling point unless we tile over D or head_dim.
I think the 4-block approach is the best for minimizing memory traffic and kernel complexity.
But let me think about one more thing: if I use 4 blocks of 256 threads, each thread has around 736k FMAs / 1024 threads = ~718 FMAs per thread. That's not a lot. The main kernel launch overhead and synchronization costs will dominate.
Actually, maybe I should use MORE threads per block to reduce iterations in the loops. With 512 threads:
- QKV all heads: 2304 * 3 = 6912 elements. Each is a dot of 96. 6912 / 512 = ~13.5 per thread.
- Attention: 2304 elements. 2304 / 512 = 4.5 per thread.
- Output: 2304 elements. 2304 / 512 = 4.5 per thread.
With 1024 threads:
- QKV: ~6.75 per thread
- Attention: ~2.25 per thread
- Output: ~2.25 per thread
But CUDA has a limit of 1024 threads per block. And shared memory of 29 KB is fine.
However, for N=24 and head_dim=24, having 1024 threads means many threads are idle during the softmax (which operates on 24 elements per row). With 1024 threads and 24 rows * 4 heads = 96 softmax rows, some rows get multiple threads helping.
Actually, in the softmax, I have a loop `for (int i = 0; i < N; i++)` inside which each thread loops over j. With 1024 threads and N=24, each thread only participates in at most 1 j per row (since 1024 > 24). Actually no, the inner loop is `for (int j = tid; j < N; j += num_threads)`, so with 1024 threads, each thread handles at most 0 or 1 j. 24 threads handle 1 j each, 1000 threads handle 0. That's very wasteful!
For the softmax, the block size should ideally be around N or slightly more, so most threads participate. With N=24, a block size of 32 or 64 is ideal for softmax. But we want larger blocks for the compute-heavy parts.
This is the classic tension. A block size of 256 might be a good compromise:
- For compute-heavy loops: good parallelism
- For softmax: 256 threads for 24 elements, so ~10% utilization. A bit wasteful but the loop runs only 24 * 4 = 96 times total.
Actually, in the 4-block approach, the softmax loop runs N=24 times per head, and H=4 heads, so 96 iterations. Each iteration has a block-level sync. That's 96 syncs per block. The sync overhead might be noticeable.
In the 16-block approach, each block does N=24 softmax rows, so 24 syncs per block. Total syncs across all blocks: 16 * 24 = 384. More total syncs but distributed.
Hmm, this is getting complicated. Let me just go with the 16-block approach, with block size 64 or 128.
Actually, I'll go with the 4-block approach and use 256 threads. The kernel is cleaner with fewer blocks, and memory traffic is minimal. If the performance isn't optimal, I can optimize later.
Wait, I've been overthinking this. Let me just write the one-big-kernel-per-batch code and test it. If it's correct and faster, great. If not, I'll iterate.
Actually, I just realized there's an even simpler 2-kernel approach:
- Kernel 1: Fused QKV + attention for all heads, one block per batch. Output is (B, N, D).
- Kernel 2: Output projection + residual, one block per batch.
For kernel 1:
1. Load normalized nodes
2. Compute Q, K, V for all heads
3. For each head:
a. Compute attention (reuse attn SMEM)
b. Compute context and accumulate into out_acc (N, D)
4. Write out_acc
For kernel 2:
1. Load out_acc
2. out = nodes + out_acc @ out_w^T
Wait, but then kernel 2 is just another matmul + add, which is what PyTorch does efficiently. Actually, I can include the output projection in kernel 1.
For kernel 1 with output projection included:
1. Load normalized nodes
2. Compute Q, K, V for all heads
3. Initialize out_acc to 0 in shared memory
4. For each head:
a. Compute attention
b. Compute context
c. partial_out = context @ out_w[head_start:head_end, :]
d. Add partial_out to out_acc
5. Write out_acc to global memory
Then kernel 2: `output = nodes + out_acc`
This is still one fused kernel per batch (4 blocks) plus a simple add kernel.
The main issue with this design is that for step 4c, `partial_out` is (N, D), computed from context (N, head_dim) and out_w_slice (head_dim, D). This is a matmul of size (24, 24) @ (24, 96). For each head.
I can compute this efficiently with the 256 threads.
OK let me write this version. It's cleaner and avoids partials buffer.
```cpp
__global__ void fused_kernel(
const float* __restrict__ nodes, // normalized nodes (B, N, D)
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
const float* __restrict__ q_w,
const float* __restrict__ k_w,
const float* __restrict__ v_w,
const float* __restrict__ rel_w1,
const float* __restrict__ rel_b1,
const float* __restrict__ rel_w2,
const float* __restrict__ rel_b2,
const float* __restrict__ out_w,
float* __restrict__ out_acc,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* norm_nodes = smem;
float* qkv = norm_nodes + N * D; // 3 * H * N * head_dim
float* attn = qkv + 3 * H * N * head_dim; // N * N
float* out_smem = attn + N * N; // N * D
float* scratch = out_smem + N * D; // blockDim.x
// Load normalized nodes
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// Compute Q, K, V for all heads
int total_qkv = H * N * head_dim;
for (int idx = tid; idx < total_qkv; idx += num_threads) {
int h = idx / (N * head_dim);
int n = (idx / head_dim) % N;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float* q_ptr = qkv + 0 * total_qkv;
float* k_ptr = qkv + 1 * total_qkv;
float* v_ptr = qkv + 2 * total_qkv;
float qv = 0, kv = 0, vv = 0;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[w_col * D + d];
kv += nd * k_w[w_col * D + d];
vv += nd * v_w[w_col * D + d];
}
q_ptr[idx] = qv;
k_ptr[idx] = kv;
v_ptr[idx] = vv;
}
__syncthreads();
// Initialize out_smem to 0
for (int idx = tid; idx < N * D; idx += num_threads) {
out_smem[idx] = 0.0f;
}
__syncthreads();
// Attention for each head
float* q_ptr = qkv + 0 * total_qkv;
float* k_ptr = qkv + 1 * total_qkv;
float* v_ptr = qkv + 2 * total_qkv;
for (int h = 0; h < H; h++) {
// Compute attention scores
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0;
for (int hd = 0; hd < head_dim; hd++) {
int q_idx = h * N * head_dim + i * head_dim + hd;
int k_idx = h * N * head_dim + j * head_dim + hd;
score += q_ptr[q_idx] * k_ptr[k_idx];
}
score *= scale;
// rel_mlp
float pf[16];
for (int p = 0; p < P; p++) {
pf[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
for (int p = 0; p < P; p++) {
float val = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
val = val * sigmoid_fast(val);
tmp[p] = val;
}
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax
for (int i = 0; i < N; i++) {
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = block_reduce_max(max_val, scratch);
float exp_sum = 0;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
for (int j = tid; j < N; j += num_threads) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Context = attn @ v, then add to out_smem via output projection
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float ctx = 0;
for (int j = 0; j < N; j++) {
float attn_val = attn[n * N + j];
for (int hd = 0; hd < head_dim; hd++) {
ctx += attn_val * v_ptr[h * N * head_dim + j * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
}
out_smem[idx] += ctx;
}
__syncthreads();
}
// Write out_smem to global memory
for (int idx = tid; idx < N * D; idx += num_threads) {
out_acc[b * N * D + idx] = out_smem[idx];
}
}
```
Wait! The context + output projection can be combined:
`partial_out[n, d] = sum_j attn[n, j] * sum_hd v[j, hd] * out_w[d, h*head_dim + hd]`
`= sum_j attn[n, j] * v_proj[j, d]`
Where `v_proj[j, d] = sum_hd v[j, hd] * out_w[d, h*head_dim + hd]`
I can precompute `v_proj` for all j, d, then do `out_smem += attn @ v_proj`.
For each head, `v_proj` is (N, D) = (24, 96). Computing it:
```cpp
for (int jd = 0; jd < N * D; jd++) {
int j = jd / D;
int d = jd % D;
float vp = 0;
for (int hd = 0; hd < head_dim; hd++) {
vp += v_ptr[...] * out_w[...];
}
v_proj[jd] = vp;
}
```
Then:
```cpp
for (int nd = 0; nd < N * D; nd++) {
int n = nd / D;
int d = nd % D;
float out_val = 0;
for (int j = 0; j < N; j++) {
out_val += attn[n * N + j] * v_proj[j * D + d];
}
out_smem[nd] += out_val;
}
```
This is two separate loops. The first computes v_proj (N*D*head_dim FMAs = 55k), the second computes the sum (N*N*D FMAs = 55k). Total 110k FMAs per head.
My earlier combined approach:
```cpp
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float ctx = 0;
for (int j = 0; j < N; j++) {
float attn_val = attn[n * N + j];
for (int hd = 0; hd < head_dim; hd++) {
ctx += attn_val * v_ptr[...] * out_w[...];
}
}
out_smem[idx] += ctx;
}
```
This computes the same thing but in one triple loop. The number of FMAs is the same (N * N * D * head_dim = 24*24*96*24 = 1.3M per head). Wait, that's 1.3M FMAs, not 110k!
My decomposition was wrong. Let me recalculate:
`partial_out[n, d] = sum_j attn[n,j] * sum_hd v[j,hd] * out_w[d, h*head_dim+hd]`
`= sum_j sum_hd attn[n,j] * v[j,hd] * out_w[d, h*head_dim+hd]`
This is indeed N * N * head_dim * D operations if done naively.
With N=24, head_dim=24, D=96:
24 * 24 * 24 * 96 = 1,327,104 FMAs per head!
For H=4 heads: 5.3M FMAs.
But wait, if I factorize:
`v_proj[j, d] = sum_hd v[j,hd] * out_w[d, h*head_dim+hd]` - this is N*D*head_dim = 24*96*24 = 55,296 FMAs
Then `out[n, d] = sum_j attn[n,j] * v_proj[j, d]` - this is N*N*D = 24*24*96 = 55,296 FMAs
Total: 110,592 FMAs per head!
My earlier combined code was wrong - it had a triple nested loop that did redundant work. I should use the factorized approach.
In the factorized approach, I need to store v_proj somewhere. v_proj for one head is N*D = 2304 floats = 9.2 KB. I can store it in shared memory.
Shared memory layout:
- norm_nodes: 2304
- qkv: 3*H*N*head_dim = 6912
- attn: 576
- v_proj: 2304
- out_smem: 2304
- scratch: 256
Total: 14656 floats = 58.6 KB
Still fine.
But wait, I can reuse qkv memory for v_proj, since v for the current head is still needed but q and k for the current head can be overwritten... Actually, no. v_proj needs v for the current head. But q and k for the current head are no longer needed after the attention scores are computed.
So I can use q's space or k's space for v_proj. q for one head is N*head_dim = 576 floats. v_proj is 2304 floats. So q's space is too small.
Actually, I can reuse the entire qkv space! After computing attention, q and k are not needed. But v is needed for v_proj. Once v_proj is computed, v is also not needed directly.
Let me think:
1. After softmax, q and k are not needed.
2. We need v to compute v_proj. v_proj = v @ W_out_slice.
3. Once v_proj is computed, v is not needed.
4. Then out += attn @ v_proj.
I can compute v_proj and store it in the qkv space, overwriting q and k.
qkv space is 3*H*N*head_dim = 6912 floats.
v_proj for one head is N*D = 2304 floats.
So I can store v_proj in, say, the first 2304 floats of qkv.
After computing v_proj for head h, I overwrite q space.
Then compute out += attn @ v_proj.
For the next head, I need q, k, v again. But I've overwritten q and k! So I can't reuse qkv across heads unless I compute q, k, v for all heads upfront and then process heads sequentially.
Actually, I can compute Q, K, V for all heads once, store them in qkv. Then for each head h, I use q[h], k[h], v[h] to compute attention, then compute v_proj from v[h], store v_proj somewhere, and accumulate out_smem.
If I store v_proj in the q[0] and k[0] slots for the current head, that would only give me 2*N*head_dim = 1152 floats. But v_proj needs N*D = 2304 floats.
Alternatively, I can compute v_proj and accumulate out_smem on the fly without storing full v_proj:
For each n, d:
`out[n,d] += sum_j attn[n,j] * sum_hd v[j,hd] * out_w[d, ...]`
We can compute this as:
```cpp
for (int j = 0; j < N; j++) {
float a = attn[n * N + j];
for (int hd = 0; hd < head_dim; hd++) {
float v_val = v_ptr[h * N * head_dim + j * head_dim + hd];
float w_val = out_w[d * D + h * head_dim + hd];
out_smem[n * D + d] += a * v_val * w_val;
}
}
```
But this is N*N*D*head_dim FMAs, which is 1.3M per head. That's inefficient.
Wait, I miscalculated. Let me re-examine.
The efficient factorization:
`partial_out[n,d] = sum_j attn[n,j] * (sum_hd v[j,hd] * out_w[d, h*head_dim+hd])`
`= sum_j attn[n,j] * v_proj[j,d]`
Where v_proj[j,d] = sum_hd v[j,hd] * out_w[d, h*head_dim+hd]
Computing v_proj: for each (j,d), dot product of head_dim elements.
- FMAs: N * D * head_dim = 24 * 96 * 24 = 55,296
Computing partial_out: for each (n,d), dot product of N elements.
- FMAs: N * D * N = 24 * 96 * 24 = 55,296
Total: 110,592 FMAs per head. For 4 heads: 442k FMAs.
The naive approach:
For each (n,d), triple sum over j and hd:
- FMAs: N * D * N * head_dim = 24 * 96 * 24 * 24 = 1,327,104 per head.
For 4 heads: 5.3M FMAs.
So the factorized approach is 12x faster! I must use it.
Now, where to store v_proj?
Each head needs v_proj of size N*D = 2304 floats.
I can store it in shared memory alongside the other arrays.
Shared memory total with v_proj explicitly stored:
- norm_nodes: 2304
- qkv: 6912
- attn: 576
- v_proj: 2304
- out_smem: 2304
- scratch: 256
Total: 14656 floats = 58.6 KB
Still fine. But qkv is 6912 floats and v_proj is 2304. Total with norm_nodes: 2304 + 6912 + 2304 + 2304 + 576 + 256 = 14656.
If I want to reduce shared memory, I can compute v_proj on the fly per (j,d) without storing the full array, and then accumulate for all n. But then I'd compute v_proj[j,d] multiple times (once per n).
Or I can loop over j first, compute v_proj for a single j across all d, then accumulate `attn[n,j] * v_proj[j,d]` for all n, d. This way I only need v_proj for one j at a time.
For j from 0 to N-1:
1. Compute v_proj_j[d] = sum_hd v[j,hd] * out_w[d, h*head_dim+hd] for all d. This is a vector of D elements.
2. For each n: out_smem[n,d] += attn[n,j] * v_proj_j[d]
If I store v_proj_j in shared memory (D elements = 96 floats), the loop is:
```cpp
for (int j = 0; j < N; j++) {
// Compute v_proj_j[d] for all d - D * head_dim FMAs
for (int d = tid; d < D; d += num_threads) {
float vp = 0;
for (int hd = 0; hd < head_dim; hd++) {
vp += v_ptr[...] * out_w[...];
}
v_proj_j[d] = vp;
}
__syncthreads();
// Accumulate for all n, d
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
out_smem[idx] += attn[n * N + j] * v_proj_j[d];
}
__syncthreads();
}
```
This requires v_proj_j to be size D = 96 floats. Very small!
And the inner loops are well-parallelized.
For j=0..N-1 (24 iterations):
- Iteration j: compute v_proj_j (96 elements) - 96 dot products of 24. With 256 threads, each thread handles a few d values.
- Then accumulate into out_smem for all n, d.
Actually, the accumulation can be done with 256 threads parallelizing over N*D=2304 elements. Each thread handles 2304/256 = 9 elements. For each element, it does one multiply and add. Very efficient.
And for computing v_proj_j, each d requires a dot product of head_dim=24. With 256 threads doing 96 dot products, each thread does about 1/3 of a dot product. A bit underutilized, but fine.
Wait, with 256 threads and D=96, `for (int d = tid; d < D; d += num_threads)` means threads 0..95 each handle 1 d value, and threads 96..255 are idle. That's 62.5% utilization for this loop.
But this loop runs 24 times (for each j). Overall, it's not too bad.
Alternatively, I can compute v_proj_j with fewer threads by having each thread handle multiple d values. But the outer loop over j already parallelizes the work.
Actually, we can restructure: have each thread precompute a slice of v_proj for all j and d, storing nothing, and then loop over j inside the n,d loop. No, that would be the inefficient triple loop.
The j-loop approach with v_proj_j buffer is good. Let me estimate the thread utilization:
For the v_proj_j computation:
- 96 useful threads out of 256. ~37.5% utilization.
For the accumulation:
- 2304 elements / 256 = 9 iterations. All threads busy.
Since the accumulation is a larger portion, overall utilization is decent.
Another optimization: compute v_proj for all j at once and store in shared memory (N*D = 2304 floats). Then the accumulation is just one loop without the __syncthreads between j iterations.
Memory: 2304 floats = 9.2 KB. With the current SMEM budget, adding this is fine.
Let me check total SMEM with v_proj array:
- norm_nodes: 2304
- qkv: 6912
- attn: 576
- v_proj: 2304
- out_smem: 2304
- scratch: 256
Total: 14656 floats = 58.6 KB
This is well within the 227 KB limit. And for 4 blocks, total SMEM used by all blocks simultaneously is 4 * 58.6 = 234.4 KB. Very close to the 228 KB per SM limit! But the 4 blocks are on different SMs, so each SM only uses 58.6 KB. Fine.
Wait, blocks are distributed to SMs, and each SM has 228 KB. If a block uses 58.6 KB, then each SM can hold at most floor(228 / 58.6) = 3 blocks. With 4 blocks total, most SMs run 1 block, and a few might run 2 or 3. Fine.
OK, I'll use the v_proj array.
But wait, I just realized: I don't need to store QKV for all heads simultaneously if I process one head at a time! This would save a lot of shared memory.
Alternative: process one head at a time, compute Q_h, K_h, V_h from norm_nodes, compute attention, compute v_proj, accumulate to out_smem, then move to next head.
Shared memory for head-by-head:
- norm_nodes: 2304
- q: 576
- k: 576
- v: 576
- attn: 576
- v_proj: 2304
- out_smem: 2304
- scratch: 256
Total: 9472 floats = 37.9 KB
Much less! And Q, K, V per head are computed on the fly from norm_nodes.
QKV computation per head: each thread handles some (n, hd) elements. Total elements per head: 576. With 256 threads, 576/256 = ~2.25 per thread. Each is a dot product of D=96. Total FMAs per thread for QKV: ~2.25 * 3 * 96 = ~648.
For all heads: 4 * 648 = ~2592 FMAs per thread. Still not a lot, but reasonable.
And we save shared memory. This seems like the better approach.
But wait, earlier I thought computing QKV per head on the fly means the block does 55k FMAs per head * 4 = 221k FMAs for QKV. With 256 threads, that's 221k / 256 = 864 FMAs per thread. Right, 864 FMAs per thread. That's a decent amount.
Actually, it's 864 FMAs for QKV. Plus attention and output, total maybe 2k-3k FMAs per thread. With 67 TFLOPS peak, minimum time is in the hundreds of nanoseconds. But with memory stalls and syncs, probably microseconds.
OK, I'll use the head-by-head approach with less shared memory. Here's the final kernel design:
```cpp
__global__ void fused_kernel(
const float* __restrict__ norm_nodes_in,
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
const float* __restrict__ q_w,
const float* __restrict__ k_w,
const float* __restrict__ v_w,
const float* __restrict__ rel_w1,
const float* __restrict__ rel_b1,
const float* __restrict__ rel_w2,
const float* __restrict__ rel_b2,
const float* __restrict__ out_w,
float* __restrict__ out_acc,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* norm_nodes = smem; // N * D
float* q = norm_nodes + N * D; // N * head_dim
float* k = q + N * head_dim; // N * head_dim
float* v = k + N * head_dim; // N * head_dim
float* attn = v + N * head_dim; // N * N
float* v_proj = attn + N * N; // N * D
float* out_smem = v_proj + N * D; // N * D
float* scratch = out_smem + N * D; // blockDim.x
// Load normalized nodes for this batch
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = norm_nodes_in[b * N * D + idx];
}
__syncthreads();
// Initialize output accumulator in shared memory
for (int idx = tid; idx < N * D; idx += num_threads) {
out_smem[idx] = 0.0f;
}
__syncthreads();
// Process each head
for (int h = 0; h < H; h++) {
// Compute Q, K, V for this head
int qkv_size = N * head_dim;
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float qv = 0, kv = 0, vv = 0;
for (int d = 0; d < D; d++) {
float nd = norm_nodes[n * D + d];
qv += nd * q_w[w_col * D + d];
kv += nd * k_w[w_col * D + d];
vv += nd * v_w[w_col * D + d];
}
q[idx] = qv;
k[idx] = kv;
v[idx] = vv;
}
__syncthreads();
// Compute attention scores + rel_bias
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0;
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// rel_mlp
float pf[16];
for (int p = 0; p < P; p++) {
pf[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
for (int p = 0; p < P; p++) {
float val = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
val = val * sigmoid_fast(val);
tmp[p] = val;
}
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax per row
for (int i = 0; i < N; i++) {
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = block_reduce_max(max_val, scratch);
float exp_sum = 0;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
for (int j = tid; j < N; j += num_threads) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Compute v_proj[j, d] = sum_hd v[j, hd] * out_w[d, h*head_dim + hd]
for (int idx = tid; idx < N * D; idx += num_threads) {
int j = idx / D;
int d = idx % D;
float vp = 0;
for (int hd = 0; hd < head_dim; hd++) {
vp += v[j * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
v_proj[idx] = vp;
}
__syncthreads();
// Accumulate: out_smem[n, d] += sum_j attn[n, j] * v_proj[j, d]
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float accum = 0;
for (int j = 0; j < N; j++) {
accum += attn[n * N + j] * v_proj[j * D + d];
}
out_smem[idx] += accum;
}
__syncthreads();
}
// Write to global memory
for (int idx = tid; idx < N * D; idx += num_threads) {
out_acc[b * N * D + idx] = out_smem[idx];
}
}
```
Shared memory:
- norm_nodes: 2304
- q: 576
- k: 576
- v: 576
- attn: 576
- v_proj: 2304
- out_smem: 2304
- scratch: 256
Total: 9472 floats = 37.9 KB
Good. And for the residual add, I'll write a tiny kernel:
```cpp
__global__ void add_residual_kernel(
const float* __restrict__ nodes,
const float* __restrict__ out_acc,
float* __restrict__ output,
int total_elements
) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= total_elements) return;
output[idx] = nodes[idx] + out_acc[idx];
}
```
And in Python:
```python
out_acc = torch.empty(B, N, D, device=nodes.device, dtype=nodes.dtype)
# ... launch fused_kernel ...
output = torch.empty_like(nodes)
# ... launch add_residual_kernel ...
return output
```
Actually, `out_acc` can be initialized and then used directly. But the kernel writes to all elements, so no initialization needed.
Wait, I need to make sure the `nodes` passed to `add_residual_kernel` is the original `nodes`, and `out_acc` is what the fused kernel wrote.
Let me finalize the code. But first, I should consider whether there's a way to avoid the second kernel entirely by passing `nodes` to the fused kernel and doing the add there.
If I load `nodes` in the fused kernel and add it to `out_smem` before writing, I can write the final output directly.
In the fused kernel:
```cpp
// Before writing, add residual
for (int idx = tid; idx < N * D; idx += num_threads) {
out_acc[b * N * D + idx] = nodes[b * N * D + idx] + out_smem[idx];
}
```
This saves a second kernel launch! And since `nodes` is already in cache (we just read `norm_nodes` from it), it should be efficient.
Wait, in my design, `norm_nodes_in` is already the result of LayerNorm, not the original nodes. To add the residual, I need the original `nodes`. So I should pass both `norm_nodes_in` and `nodes` to the kernel.
Actually, yes: in `ModelNew.forward`:
```python
h = self.norm(nodes)
# launch fused_kernel with h and nodes
```
And in the kernel, at the end:
```cpp
// Write final output = original_nodes + out_proj(context)
for (int idx = tid; idx < N * D; idx += num_threads) {
output[b * N * D + idx] = nodes[b * N * D + idx] + out_smem[idx];
}
```
This eliminates the second kernel entirely!
So the final design is:
1. PyTorch LayerNorm
2. Custom fused kernel (attention + output proj + residual)
2 kernels total.
Let me rewrite the kernel signature:
```cpp
__global__ void fused_kernel(
const float* __restrict__ norm_nodes_in, // LayerNorm output (B, N, D)
const float* __restrict__ orig_nodes, // Original nodes for residual (B, N, D)
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
... weights ...,
float* __restrict__ output, // Final output (B, N, D)
...
)
```
This is great. One custom kernel, plus PyTorch LayerNorm.
Actually, I can also integrate LayerNorm into the custom kernel to reduce to ONE total kernel. But PyTorch's LayerNorm is very efficient, and integrating it adds complexity. Plus we only save one kernel launch.
Hmm, but one kernel launch for LayerNorm plus one custom kernel = 2 kernels. With PyTorch reference having ~12, that's still a 6x reduction.
Actually wait, the PyTorch LayerNorm kernel is the first of ~12. I already cut it down to 2 (LayerNorm + custom). That's good enough.
But if I want to be even faster, I can integrate LayerNorm. Each block would load `nodes`, compute LayerNorm, store in shared memory, and then proceed. This saves writing/reading the normalized nodes.
Memory traffic with PyTorch LayerNorm:
- LayerNorm reads nodes (36 KB), writes h (36 KB)
- Custom kernel reads h (36 KB), writes output (36 KB)
- Plus pair_feats, edge_mask, etc.
Total for nodes/h: 144 KB (3 reads, 1 write of 36KB each)
With integrated LayerNorm:
- Custom kernel reads nodes (36 KB), writes output (36 KB)
- pair_feats, edge_mask, weights
Total for nodes: 72 KB (1 read, 1 write)
So integrated LayerNorm saves 72 KB of memory traffic. For a workload this small, maybe it doesn't matter much. But it also saves a kernel launch, which might matter.
The cost: each of the 4 blocks computes LayerNorm for its batch. With N=24 and D=96, LayerNorm is not expensive. And since a block handles one batch, there's no redundancy.
Let me integrate LayerNorm to make it fully fused!
LayerNorm computation per batch:
For each node n:
1. Compute mean: sum_d nodes[n, d] / D
2. Compute variance: sum_d (nodes[n, d] - mean)^2 / D
3. Normalize: (nodes[n, d] - mean) / sqrt(var + 1e-5) * weight[d] + bias[d]
This is done in shared memory.
Here's the updated kernel:
```cpp
__global__ void fused_kernel(
const float* __restrict__ nodes, // (B, N, D)
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
const float* __restrict__ norm_w,
const float* __restrict__ norm_b,
const float* __restrict__ q_w,
const float* __restrict__ k_w,
const float* __restrict__ v_w,
const float* __restrict__ rel_w1,
const float* __restrict__ rel_b1,
const float* __restrict__ rel_w2,
const float* __restrict__ rel_b2,
const float* __restrict__ out_w,
float* __restrict__ output,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* norm_nodes = smem; // N * D
float* q = norm_nodes + N * D; // N * head_dim
float* k = q + N * head_dim; // N * head_dim
float* v = k + N * head_dim; // N * head_dim
float* attn = v + N * head_dim; // N * N
float* v_proj = attn + N * N; // N * D
float* out_smem = v_proj + N * D; // N * D
float* scratch = out_smem + N * D; // blockDim.x
float* reduction_smem = scratch + num_threads; // N * 2 (for mean/var reduction)
// Load nodes and compute LayerNorm
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_nodes[idx] = nodes[b * N * D + idx];
}
__syncthreads();
// LayerNorm per node
for (int n = 0; n < N; n++) {
float sum = 0.0f, sq_sum = 0.0f;
for (int d = tid; d < D; d += num_threads) {
float val = norm_nodes[n * D + d];
sum += val;
sq_sum += val * val;
}
sum = block_reduce_sum(sum, scratch);
sq_sum = block_reduce_sum(sq_sum, scratch);
float mean = sum / D;
float var = sq_sum / D - mean * mean;
float inv_std = rsqrtf(var + 1e-5f);
for (int d = tid; d < D; d += num_threads) {
float val = norm_nodes[n * D + d];
norm_nodes[n * D + d] = (val - mean) * inv_std * norm_w[d] + norm_b[d];
}
__syncthreads();
}
// Initialize output accumulator
for (int idx = tid; idx < N * D; idx += num_threads) {
out_smem[idx] = 0.0f;
}
__syncthreads();
// Then proceed with attention as before...
}
```
Wait, the LayerNorm requires `block_reduce_sum` twice per node (for sum and sq_sum). For N=24 nodes, that's 48 block reductions. With `__syncthreads()` in each reduction, that's 48 * ~9 = ~432 syncs per block. That's a lot of syncs, but each sync is cheap.
Actually, I can optimize LayerNorm by having threads compute both sum and sq_sum in one pass, then do two reductions. Still 24 passes and 48 reductions.
Since D=96 and block_size=256, only 96 threads participate in the sum loop. The rest are idle but participate in the reduction. This is somewhat wasteful.
Alternatively, I can have each thread handle one node, computing the full mean/var for that node. Then no syncs are needed for LayerNorm!
With 256 threads and N=24 nodes, I can assign roughly 10 threads per node. Or simpler: have 24 threads handle one node each, and the remaining 232 threads are idle.
But 24 < 32 (warp size), so one warp handles 24 nodes, each thread doing one node. The other 7 warps are completely idle for LayerNorm. That's very wasteful.
Or: each thread handles multiple nodes. 256 threads / 24 nodes = ~10.7 threads per node. Each thread handles a subset of D=96 elements for its assigned node, then a small reduction within the node group.
This gets complicated. Maybe using PyTorch LayerNorm is cleaner and not much slower for this small size. Let me stick with PyTorch LayerNorm for simplicity.
Actually, I just realized: if I use 4 blocks of 64 threads each, and assign one warp (32 threads) per node for LayerNorm... no, still complicated.
Let me just use PyTorch for LayerNorm. It's good enough.
Actually, wait. Let me consider yet another design: use 4 blocks, each with D=96 threads. Each block handles one batch.
With 96 threads:
- LayerNorm: each node gets 96/24 = 4 threads. Actually, D=96 for each of N=24 nodes.
- For node n, threads handle the D elements. But if I only have 96 threads, I could assign 4 threads per node (96/24 = 4). Each of 4 threads handles 24 elements of a node. Then I need a reduction of 4 threads.
- Or assign 1 thread per dimension across all nodes. Not straightforward.
Alternatively, block size = N = 24 threads. Each thread handles one node, computing the full mean, var, and normalized values over D=96 elements. No sync needed for LayerNorm!
With 24 threads per block, N=24 nodes, each thread handles exactly one node.
For QKV: each thread computes all (h, hd) elements for its node? That would be H*head_dim = 96 elements per node, each a dot of D=96. That's 96*96 = 9216 FMAs per thread for Q. A lot of work per thread. Total 24 threads * 9216 FMAs = 221k FMAs per block. Good.
For attention: each thread handles one query row. For query i, thread i computes all N=24 scores (dots of head_dim=24), does softmax over 24 elements, then computes context for all D=96 elements.
Wait, if each thread handles one query, it computes:
1. Attn scores for all j: 24 dots of 24 = 576 FMAs
2. rel_mlp for 24 edges: 24 * small MLP
3. Softmax over 24 values: sequential on one thread
4. Context = attn @ v: 24 dots of 24 dims for each hd? No, context is (head_dim,), so 24 dots of 24 = 576 FMAs.
5. Output projection partial: context (head_dim) @ out_w[:, :] = 24 dots of D=96 dims each = 2304 FMAs.
Wait, the output projection should be done after accumulating all heads. Or each thread can accumulate its partial output across heads.
With block size = 24 threads (one per query), and 4 heads:
For each head h:
- Thread i computes score[i, j] for j=0..N-1 - 576 FMAs
- Computes 24 rel_mlps
- Softmax - sequential
- Computes context[i, :] - 576 FMAs
- Accumulates into out_smem[i, :] += context[i, :] @ out_w[h*head_dim:(h+1)*head_dim, :] - 2304 FMAs
Total per head per thread: ~3500 FMAs. For 4 heads: ~14k FMAs. Total per block: 24 * 14k = 336k FMAs.
This is actually very efficient in terms of thread utilization! Each thread has a lot of sequential work, which is great for avoiding syncs. And the block has only 24 threads, so no wasted threads in the softmax (each thread handles exactly one row).
But with only 24 threads per block and 4 blocks, total threads = 96. This is very few for an H100. The GPU would be severely underutilized. Each SM can run up to ~2048 threads, but we only have 96 total. That's about 5% occupancy at best.
Hmm, but if each thread has 14k FMAs, and the H100 can do about 2000 FMAs per cycle per SM (?), the runtime per thread is ~7 cycles. But with instruction latency and memory stalls, maybe thousands of cycles.
Actually, the H100 can issue many instructions per cycle per warp, but with only 3 warps (24 threads = 1 warp on some arch, actually 24 threads is not a full warp). Wait, a warp is 32 threads. 24 threads per block means the block has 1 warp (with 8 threads masked/tailing). Actually, CUDA still executes it as a warp of 32, but only 24 threads have valid work. The other 8 are "inactive" but still consume resources.
So effectively, 4 blocks * 1 warp = 4 warps total. With ~132 SMs, each SM runs at most 1 warp. That's extremely poor utilization.
So block size of 24 is not good for GPU utilization.
Let me try a block size of 128. For LayerNorm, it's wasteful. For softmax, it's somewhat wasteful. But for compute-heavy parts, it's good.
Actually, let me just go with 256 threads and accept the softmax inefficiency. The key metric is end-to-end time, and fusing everything will help a lot.
OK, final decision:
- PyTorch LayerNorm
- Custom fused kernel with 4 blocks (one per batch), 256 threads each
- The custom kernel handles QKV per head, attention, output projection, and residual add
Wait, I realize I should also consider using 4 blocks with 128 threads. The compute-heavy parts are more balanced, and softmax reductions have better utilization (128 threads for 24 elements - first reduction step only 64 threads useful, then 32, then 16, 8, 4, 2, 1 - about 7 steps).
With 64 threads: 32 useful in first step, then 16, 8, 4, 2, 1 - 6 steps.
With 256 threads: 128, 64, 32, 16, 8, 4, 2, 1 - 8 steps. But 232 threads are idle for the first step.
Actually, for a tree reduction, the number of steps is log2(block_size). The total work is proportional to block_size (each thread writes to shared mem). So smaller block sizes have fewer steps but also fewer total threads for the compute-heavy parts.
I think 128 threads is a good middle ground. Let me use 128.
Actually, let me just test with 256 first and see. I can change it later.
One more concern: `load_inline` compilation time. I should make the CUDA code as small and clean as possible to compile quickly.
Let me write the complete Python module and test compilation.
Actually, I realize that for the `reduce_and_residual_kernel` (if I use it), or for the atomic approach, the code is simpler. Let me reconsider.
Atomic approach:
- 16 blocks, one per (batch, head)
- Each block computes attention + context, then does output projection and atomic-add to output buffer
- Then a simple element-wise kernel adds residual
16 blocks is much better for GPU utilization than 4 blocks! And the atomic-add contention is minimal (4 threads write to each output element on average).
Let me estimate performance:
- Kernel 1: 16 blocks * 256 threads = 4096 total threads
- Kernel 2: 4 * 96 / 256 = ~2 blocks * 256 threads = 512 threads
Total threads: 4608. Still not great for 132 SMs, but better.
Memory traffic:
- LayerNorm: read nodes (36KB), write h (36KB)
- Kernel 1: 16 blocks each read h (36KB * 16? No - only 4 batches, 4 blocks per batch read the same h). So total h reads: 4*36KB = 144KB. Pair_feats: 4*147KB = 588KB. Edge_mask: 4*9KB = 36KB. Output buffer atomics: 4*36KB = 144KB (read-modify-write).
- Kernel 2: read output buffer (36KB), read nodes (36KB), write output (36KB)
Total: ~1 MB.
Wait, each of the 4 blocks per batch reads the same h. But when they run on different SMs, the L2 cache might serve the second, third, and fourth reads if the first block's read is still in cache. Since h is only 36 KB, it easily fits in L2. So the actual memory bandwidth might not be 4x.
Similarly for pair_feats and edge_mask.
So the 16-block approach might not be that bad for memory traffic in practice.
And it has much better parallelism. With 16 blocks, we can keep up to 16 SMs busy simultaneously.
I'll go with the 16-block atomic approach. The code is also simpler (no v_proj array, no accumulation across heads in shared memory).
Final kernel for 16-block atomic approach:
```cpp
__global__ void fused_attention_head_kernel(
const float* __restrict__ norm_nodes,
const float* __restrict__ pair_feats,
const float* __restrict__ edge_mask,
const float* __restrict__ q_w,
const float* __restrict__ k_w,
const float* __restrict__ v_w,
const float* __restrict__ rel_w1,
const float* __restrict__ rel_b1,
const float* __restrict__ rel_w2,
const float* __restrict__ rel_b2,
const float* __restrict__ out_w,
float* __restrict__ out_acc,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int b = blockIdx.x / H;
int h = blockIdx.x % H;
int tid = threadIdx.x;
int num_threads = blockDim.x;
extern __shared__ float smem[];
float* q = smem;
float* k = q + N * head_dim;
float* v = k + N * head_dim;
float* attn = v + N * head_dim;
float* scratch = attn + N * N;
float* norm_buf = scratch + num_threads;
// Load norm_nodes for this batch
for (int idx = tid; idx < N * D; idx += num_threads) {
norm_buf[idx] = norm_nodes[b * N * D + idx];
}
__syncthreads();
// Compute Q, K, V for this head
int qkv_size = N * head_dim;
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
int w_col = h * head_dim + hd;
float qv = 0, kv = 0, vv = 0;
for (int d = 0; d < D; d++) {
float nd = norm_buf[n * D + d];
qv += nd * q_w[w_col * D + d];
kv += nd * k_w[w_col * D + d];
vv += nd * v_w[w_col * D + d];
}
q[idx] = qv;
k[idx] = kv;
v[idx] = vv;
}
__syncthreads();
// Attention scores
for (int idx = tid; idx < N * N; idx += num_threads) {
int i = idx / N;
int j = idx % N;
float score = 0;
for (int hd = 0; hd < head_dim; hd++) {
score += q[i * head_dim + hd] * k[j * head_dim + hd];
}
score *= scale;
// rel_mlp
float pf[16];
for (int p = 0; p < P; p++) {
pf[p] = pair_feats[b * N * N * P + i * N * P + j * P + p];
}
float tmp[16];
for (int p = 0; p < P; p++) {
float val = rel_b1[p];
for (int pp = 0; pp < P; pp++) {
val += pf[pp] * rel_w1[p * P + pp];
}
// SiLU
float s = val;
val = s / (1.0f + expf(-s));
val = s * val;
tmp[p] = val;
}
float bias = rel_b2[h];
for (int p = 0; p < P; p++) {
bias += tmp[p] * rel_w2[h * P + p];
}
score += bias;
if (edge_mask[b * N * N + i * N + j] < 0.5f) {
score = -1e9f;
}
attn[idx] = score;
}
__syncthreads();
// Softmax
for (int i = 0; i < N; i++) {
float max_val = -1e9f;
for (int j = tid; j < N; j += num_threads) {
max_val = fmaxf(max_val, attn[i * N + j]);
}
max_val = block_reduce_max(max_val, scratch);
float exp_sum = 0;
for (int j = tid; j < N; j += num_threads) {
float e = expf(attn[i * N + j] - max_val);
attn[i * N + j] = e;
exp_sum += e;
}
exp_sum = block_reduce_sum(exp_sum, scratch);
for (int j = tid; j < N; j += num_threads) {
attn[i * N + j] /= exp_sum;
}
__syncthreads();
}
// Context + output projection partial, atomic add
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float ctx = 0;
for (int hd = 0; hd < head_dim; hd++) {
float c = 0;
for (int j = 0; j < N; j++) {
c += attn[n * N + j] * v[j * head_dim + hd];
}
ctx += c * out_w[d * D + h * head_dim + hd];
}
atomicAdd(&out_acc[b * N * D + idx], ctx);
}
}
```
Wait, the context computation can be optimized:
`partial[n, d] = sum_hd [sum_j attn[n,j] * v[j,hd]] * out_w[d, h*head_dim+hd]`
I can compute `context[n, hd] = sum_j attn[n,j] * v[j,hd]` first, then dot with out_w.
```cpp
// Compute context into q (reuse)
for (int idx = tid; idx < qkv_size; idx += num_threads) {
int n = idx / head_dim;
int hd = idx % head_dim;
float c = 0;
for (int j = 0; j < N; j++) {
c += attn[n * N + j] * v[j * head_dim + hd];
}
q[idx] = c; // reuse q as context buffer
}
__syncthreads();
// Output projection partial
for (int idx = tid; idx < N * D; idx += num_threads) {
int n = idx / D;
int d = idx % D;
float partial = 0;
for (int hd = 0; hd < head_dim; hd++) {
partial += q[n * head_dim + hd] * out_w[d * D + h * head_dim + hd];
}
atomicAdd(&out_acc[b * N * D + idx], partial);
}
```
This is cleaner: first compute context (N*head_dim elements), then output projection (N*D elements). The context is stored in shared memory (reusing q).
Total shared memory:
- q (context): N*head_dim = 576
- k: 576 (unused after attention, but we keep it)
- v: 576
- attn: 576
- scratch: 256
- norm_buf: N*D = 2304
Total: 4864 floats = 19.5 KB
Wait, k is unused after the attention scores are computed. But since we're reusing q for context, and k is just sitting there, that's fine. We don't need to reclaim it.
Actually, if we stored norm_buf at the beginning of smem, then q, k, v, attn, scratch after it:
- norm_buf: 2304
- q: 576
- k: 576
- v: 576
- attn: 576
- scratch: 256
Total: 4864 floats = 19.5 KB
This is very compact.
OK, I'm going with this 16-block atomic approach. Then a simple element-wise add kernel for residual.
Combined, the code is simpler and should have good performance.
Let me write the final complete Python code.
Wait, I want to double check the rel_mlp SiLU implementation:
```cpp
float s = val;
val = s / (1.0f + expf(-s));
val = s * val;
```
This is `val = s * sigmoid(s)` which is SiLU(x) = x * sigmoid(x). Correct.
And for the edge mask comparison: the reference does `edge_mask.unsqueeze(1) < 0.5`. Since edge_mask is (B, N, N), unsqueezing gives (B, 1, N, N), and comparing gives a bool tensor. `masked_fill` fills where True, i.e., where edge_mask < 0.5, meaning edge_mask is 0.
So: `if (edge_mask[b*N*N + i*N + j] < 0.5f) score = -1e9f;`
This matches.
One final thing: the `load_inline` call will compile the CUDA code. I should make sure the host-side wrapper functions are listed in `functions=[...]`.
Let me write the complete code now.
Actually, I just want to verify: is it `functions=["fused_attention_head_kernel", "add_residual_kernel"]` or do I need to specify the exact names? Yes, the names should match the CUDA kernel function names.
And I need to know how to pass grid/block parameters. In `load_inline`, the generated Python bindings accept the arguments plus `grid` and `block` kwargs. Let me verify the standard pattern.
Actually, when `load_inline` generates bindings, it creates Python functions that accept the function arguments and then `grid`, `block`, and optionally `stream`, `shared_mem`. So I can call:
```python
cuda_ext.fused_attention_head_kernel(
arg1, arg2, ...,
grid=(grid_size,), block=(block_size,), shared_mem=shared_mem_size
)
```
This should work.
Now, I also need to consider whether I should use `torch::Tensor` types in the CUDA code rather than raw pointers. But `load_inline` generates pybind11 wrappers from the C++ signatures. If I use raw `float*` pointers, I need to pass the `.data_ptr<float>()` from Python.
Actually, I think `load_inline` can handle `torch::Tensor` arguments directly. Let me check common usage patterns.
Looking at common `load_inline` usage:
```cpp
void my_kernel(torch::Tensor input, torch::Tensor output) {
float* in_ptr = input.data_ptr<float>();
float* out_ptr = output.data_ptr<float>();
kernel<<<...>>>(in_ptr, out_ptr, ...);
}
```
And functions=["my_kernel"].
This is the standard pattern. The C++ wrapper handles tensor arguments.
But wait, in my case, I wrote the kernels directly without a C++ wrapper that takes tensors. The `load_inline` with empty `cpp_sources` might not generate the wrappers properly.
Hmm, let me check. `load_inline` auto-generates a pybind11 module. For functions listed in `functions=[...]`, it binds the CUDA functions as-is. But CUDA kernel functions can't be called directly from Python! They need to be launched with `<<< >>>` syntax, which requires a C++ wrapper function.
Actually, looking at the KernelBench examples and the standard usage, we DO need a C++ wrapper function that takes tensors, extracts pointers, and launches the kernel. The wrapper is part of what we write in `cpp_sources` or `cuda_sources`. Since `load_inline` with `cpp_sources=""` might not work, or the wrappers might need to be in `cpp_sources`.
Wait, let me re-read the problem statement: "Do not write your own `PYBIND11_MODULE(...)` block in `cpp_sources`. `load_inline` auto-generates one from the `functions=[...]` argument."
This implies that `functions=[...]` lists C++ functions (not CUDA kernels) that should be exposed to Python. These C++ functions can launch CUDA kernels.
So I need to write C++ wrapper functions that launch the kernels, and list those wrapper functions in `functions=[...]`.
Let me write:
```cpp
// C++ wrappers
void launch_fused_attention(
torch::Tensor norm_nodes,
torch::Tensor pair_feats,
torch::Tensor edge_mask,
torch::Tensor q_w,
torch::Tensor k_w,
torch::Tensor v_w,
torch::Tensor rel_w1,
torch::Tensor rel_b1,
torch::Tensor rel_w2,
torch::Tensor rel_b2,
torch::Tensor out_w,
torch::Tensor out_acc,
int B, int N, int D, int P, int H, int head_dim, float scale
) {
int block_size = 256;
int grid_size = B * H;
int shared_mem = (N * D + 3 * N * head_dim + N * N + block_size) * sizeof(float);
fused_attention_head_kernel<<<grid_size, block_size, shared_mem>>>(
norm_nodes.data_ptr<float>(),
pair_feats.data_ptr<float>(),
edge_mask.data_ptr<float>(),
q_w.data_ptr<float>(),
k_w.data_ptr<float>(),
v_w.data_ptr<float>(),
rel_w1.data_ptr<float>(),
rel_b1.data_ptr<float>(),
rel_w2.data_ptr<float>(),
rel_b2.data_ptr<float>(),
out_w.data_ptr<float>(),
out_acc.data_ptr<float>(),
B, N, D, P, H, head_dim, scale
);
}
void launch_add_residual(
torch::Tensor nodes,
torch::Tensor out_acc,
torch::Tensor output,
int total_elements
) {
int block_size = 256;
int grid_size = (total_elements + block_size - 1) / block_size;
add_residual_kernel<<<grid_size, block_size>>>(
nodes.data_ptr<float>(),
out_acc.data_ptr<float>(),
output.data_ptr<float>(),
total_elements
);
}
```
And list `functions=["launch_fused_attention", "launch_add_residual"]`.
The CUDA kernels go in `cuda_sources`, and the C++ wrappers go in `cpp_sources`.
Wait, but the problem says `cpp_sources` should not contain `PYBIND11_MODULE`. And `load_inline` auto-generates it. So I can put the C++ wrappers in `cpp_sources` and the CUDA kernels in `cuda_sources`.
Let me double check: does `load_inline` parse `cuda_sources` for device code and `cpp_sources` for host code? I think yes.
So my setup would be:
- `cpp_sources`: the C++ wrapper functions
- `cuda_sources`: the CUDA kernels
- `functions`: ["launch_fused_attention", "launch_add_residual"]
Actually, looking at PyTorch docs, `load_inline` signature is:
```python
load_inline(name, cpp_sources, cuda_sources=None, functions=None, ...)
```
The `functions` list tells it which functions from the combined source to bind. So I should have the wrapper functions available.
Let me write the complete code.
Wait, one more thing. The tensors might not be contiguous or might not be on CUDA. I should ensure they are by using `.contiguous()` in Python before passing.
Also, weight tensors from PyTorch modules might have specific strides but are typically contiguous.
OK, here's the final code I'll compile and test: