KDA源码剖析之FLA

鱿鱼圈 Lv4

前置知识:KDA算法公式、triton、cuda、c++

0 公式回顾

1 整体梳理

代码链接:flash-linear-attention/fla/ops/kda/chunk_fwd.py at main · fla-org/flash-linear-attention

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def chunk_kda_fwd(
q: torch.Tensor,
k: torch.Tensor,
v: torch.Tensor,
g: torch.Tensor,
beta: torch.Tensor,
scale: float,
initial_state: torch.Tensor,
output_final_state: bool,
state_v_first: bool = False,
cu_seqlens: torch.LongTensor | None = None,
cu_seqlens_cpu: torch.LongTensor | None = None,
chunk_indices: torch.LongTensor | None = None,
chunk_size: int = 64,
safe_gate: bool = False,
lower_bound: float | None = None,
use_gate_in_kernel: bool = False,
A_log: torch.Tensor | None = None,
dt_bias: torch.Tensor | None = None,
disable_recompute: bool = False,
return_intermediate_states: bool = False,
cp_context: FLACPContext | None = None,
):
# Apply gate activation
g_org = None
if use_gate_in_kernel:
g_org = g
g = kda_gate_chunk_cumsum(
g=g_org,
A_log=A_log,
dt_bias=dt_bias,
scale=RCP_LN2,
chunk_size=chunk_size,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
lower_bound=lower_bound,
)
else:
g = chunk_local_cumsum(
g=g,
scale=RCP_LN2,
chunk_size=chunk_size,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices
)

# qg = None if disable_recompute is False
w, u, qg, kg, Aqk, Akk = chunk_kda_fwd_intra(
q=q,
k=k,
v=v,
gk=g,
beta=beta,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_size=chunk_size,
chunk_indices=chunk_indices,
safe_gate=safe_gate,
disable_recompute=disable_recompute
)

if cp_context is not None:
initial_state = chunk_gated_delta_rule_fwd_h_pre_process(
k=kg,
w=w,
u=u,
gk=g,
cu_seqlens=cu_seqlens,
initial_state=initial_state,
context=cp_context,
chunk_size=chunk_size,
state_v_first=state_v_first,
)

h, v_new, final_state = chunk_gated_delta_rule_fwd_h(
k=kg,
w=w,
u=u,
gk=g,
initial_state=initial_state,
output_final_state=output_final_state,
cu_seqlens=cu_seqlens,
cu_seqlens_cpu=cu_seqlens_cpu,
chunk_indices=chunk_indices,
chunk_size=chunk_size,
state_v_first=state_v_first,
)

if cp_context is not None:
# In Context Parallel (CP) mode, global initial states are not supported at the entry point.
# The `initial_state` here is computed internally via inter-rank communication.
# Since only the first sequence in the local batch can be a continuation of a cross-rank sequence,
# only the first state in the tensor is relevant. We compress it to optimize memory for `save_for_backward`.
initial_state = compress_h0(initial_state, context=cp_context)

o = chunk_gla_fwd_o_gk(
q=q,
v=v_new,
g=g,
A=Aqk,
h=h,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_size=chunk_size,
chunk_indices=chunk_indices,
state_v_first=state_v_first,
)
if disable_recompute is False:
# Delete to save memory
w, u, qg, kg, v_new = None, None, None, None, None
if not return_intermediate_states:
h = None
if use_gate_in_kernel:
g = None
return o, final_state, g, Aqk, Akk, w, u, qg, kg, v_new, h, initial_state
  • kda_gate_chunk_cumsum / chunk_local_cumsum : chunk 间独立,各 chunk 内做局部前缀和 (cumsum)
  • chunk_kda_fwd_intra:chunk 间独立,各 chunk 内计算 w, u, qg, kg, Aqk, Akk
    • chunk_kda_fwd_kernel_intra_sub_chunk / chunk_kda_fwd_intra_token_parallel :计算各chunk的 Aqk, Akkd 的对角 token
    • chunk_kda_fwd_kernel_inter_solve_fused :计算各chunk的 Aqk, Akkd 的非对角 token + 计算 (I + Akk)的逆矩阵
    • recompute_w_u_fwd :计算w, u, qg, kg
  • chunk_gated_delta_rule_fwd_h : chunk 间有依赖,串行执行,计算 h(状态矩阵), v_new, final_state
  • chunk_gla_fwd_o_gk :计算输出 o

2 kda_gate_chunk_cumsum

代码链接:flash-linear-attention/fla/ops/kda/gate.py at main · fla-org/flash-linear-attention

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
@input_guard
def kda_gate_chunk_cumsum(
g: torch.Tensor,
A_log: torch.Tensor,
chunk_size: int,
scale: float = None,
dt_bias: torch.Tensor | None = None,
cu_seqlens: torch.Tensor | None = None,
output_dtype: torch.dtype | None = torch.float,
chunk_indices: torch.LongTensor | None = None,
lower_bound: float | None = None,
**kwargs,
) -> torch.Tensor:
if cu_seqlens is not None:
assert g.shape[0] == 1, "Only batch size 1 is supported when cu_seqlens are provided"
assert len(g.shape) == 4
B, T, H, S = g.shape
BT = chunk_size
if chunk_indices is None and cu_seqlens is not None:
chunk_indices = prepare_chunk_indices(cu_seqlens, BT)
NT = triton.cdiv(T, BT) if cu_seqlens is None else len(chunk_indices)
assert chunk_size == 2**(chunk_size.bit_length()-1), "chunk_size must be a power of 2"

g_org, g = g, torch.empty_like(g, dtype=output_dtype or g.dtype)
def grid(meta): return (triton.cdiv(meta['S'], meta['BS']), NT, B * H)
kda_gate_chunk_cumsum_vector_kernel[grid](
s=g_org,
A_log=A_log,
dt_bias=dt_bias,
o=g,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
lower_bound=lower_bound,
T=T,
H=H,
S=S,
BT=BT,
REVERSE=False,
)
return g

@triton.heuristics({
"HAS_BIAS": lambda args: args["dt_bias"] is not None,
'HAS_SCALE': lambda args: args['scale'] is not None,
'IS_VARLEN': lambda args: args['cu_seqlens'] is not None,
'USE_LOWER_BOUND': lambda args: args['lower_bound'] is not None,
})
@fla_cache_autotune(
configs=[
triton.Config({'BS': BS}, num_warps=num_warps)
for BS in BS_LIST
for num_warps in [2, 4, 8]
],
key=['H', 'S', 'BT', 'IS_VARLEN', 'REVERSE'],
**autotune_cache_kwargs,
)
@triton.jit(do_not_specialize=['T'])
def kda_gate_chunk_cumsum_vector_kernel(
s,
A_log,
dt_bias,
o,
scale,
cu_seqlens,
chunk_indices,
lower_bound,
T,
H: tl.constexpr,
S: tl.constexpr,
BT: tl.constexpr,
BS: tl.constexpr,
REVERSE: tl.constexpr,
HAS_BIAS: tl.constexpr,
HAS_SCALE: tl.constexpr,
IS_VARLEN: tl.constexpr,
USE_LOWER_BOUND: tl.constexpr,
):
i_s, i_t, i_bh = tl.program_id(0), tl.program_id(1), tl.program_id(2)
i_b, i_h = i_bh // H, i_bh % H
if IS_VARLEN:
i_n, i_t = tl.load(chunk_indices + i_t * 2).to(tl.int32), tl.load(chunk_indices + i_t * 2 + 1).to(tl.int32)
bos, eos = tl.load(cu_seqlens + i_n).to(tl.int32), tl.load(cu_seqlens + i_n + 1).to(tl.int32)
T = eos - bos
else:
bos, eos = i_b * T, i_b * T + T

p_s = tl.make_block_ptr(s + (bos * H + i_h) * S, (T, S), (H*S, 1), (i_t * BT, i_s * BS), (BT, BS), (1, 0))
p_o = tl.make_block_ptr(o + (bos * H + i_h) * S, (T, S), (H*S, 1), (i_t * BT, i_s * BS), (BT, BS), (1, 0))
# [BT, BS]
b_s = tl.load(p_s, boundary_check=(0, 1)).to(tl.float32)

# Apply dt_bias if exists
if HAS_BIAS:
p_b = tl.make_block_ptr(dt_bias + i_h * S, (S,), (1,), (i_s * BS,), (BS,), (0,))
b_bias = tl.load(p_b, boundary_check=(0,)).to(tl.float32)
b_s = b_s + b_bias[None, :]

b_A = tl.load(A_log + i_h).to(tl.float32)
if not USE_LOWER_BOUND:
# Apply gate: -exp(A_log) * softplus(g + bias)
b_gate = -exp(b_A) * softplus(b_s)
else:
b_gate = lower_bound * tl.sigmoid(exp(b_A) * b_s)

# Apply chunk local cumsum
if REVERSE:
b_o = tl.cumsum(b_gate, axis=0, reverse=True)
else:
b_o = tl.cumsum(b_gate, axis=0)

if HAS_SCALE:
b_o *= scale
tl.store(p_o, b_o.to(p_o.dtype.element_ty), boundary_check=(0, 1))

线程分配:

1
2
3
grid = (triton.cdiv(meta['S'], meta['BS']), NT, B * H)
# ↑ ↑ ↑
# X轴 Y轴 Z轴

每个 block 负责计算一个(batch, head, chunk, s/BS),对于一个block来说,负责 [BT, K] 的大小(这里的 BS 就是 BK)

image

沿 BT 维度做前缀和即可

3 chunk_kda_fwd_intra

代码链接:flash-linear-attention/fla/ops/kda/chunk_intra.py at main · fla-org/flash-linear-attention

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
def chunk_kda_fwd_intra(
q: torch.Tensor,
k: torch.Tensor,
v: torch.Tensor,
gk: torch.Tensor | None = None,
beta: torch.Tensor | None = None,
scale: float | None = None,
cu_seqlens: torch.LongTensor | None = None,
chunk_size: int = 64,
chunk_indices: torch.LongTensor | None = None,
safe_gate: bool = False,
disable_recompute: bool = False,
):
B, T, H, K, HV = *k.shape, gk.shape[2]
BT = chunk_size
if BT not in (32, 64):
raise ValueError(f"KDA intra chunk kernel only supports chunk_size 32 or 64, got {BT}.")
BC = 16
if chunk_indices is None and cu_seqlens is not None:
chunk_indices = prepare_chunk_indices(cu_seqlens, BT)
NT = triton.cdiv(T, BT) if cu_seqlens is None else len(chunk_indices)
NC = triton.cdiv(BT, BC)

Aqk = torch.empty(B, T, HV, BT, device=k.device, dtype=k.dtype)
# Akk must be zero-initialized - kernel only writes lower triangular
Akk = torch.zeros(B, T, HV, BT, device=k.device, dtype=k.dtype)
# Separate fp32 buffer for diagonal 16x16 blocks (for precision in solve_tril)
Akkd = torch.empty(B, T, HV, BC, device=k.device, dtype=torch.float32)

# Step 1: Run token_parallel first to compute diagonal blocks into Akkd (fp32)
# Step 1: compute diagonal blocks into Akk_diag (fp32)
if safe_gate:
grid = (NT, NC, B * HV)
BK = triton.next_power_of_2(K)
chunk_kda_fwd_kernel_intra_sub_chunk[grid](
q=q,
k=k,
g=gk,
beta=beta,
Aqk=Aqk,
Akk=Akkd,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
T=T,
H=H,
HV=HV,
K=K,
BT=BT,
BC=BC,
BK=BK,
USE_GATHER=IS_GATHER_SUPPORTED,
)
else:
Aqk, Akkd = chunk_kda_fwd_intra_token_parallel(
q=q,
k=k,
gk=gk,
beta=beta,
Aqk=Aqk,
Akk=Akkd,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_size=BT,
sub_chunk_size=BC,
)

# Step 2: Fused inter + solve_tril (works for both fixed-len and varlen)
grid = (NT, B * HV)
chunk_kda_fwd_kernel_inter_solve_fused[grid](
q=q,
k=k,
g=gk,
beta=beta,
Aqk=Aqk,
Akkd=Akkd,
Akk=Akk,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
T=T,
H=H,
HV=HV,
K=K,
BT=BT,
BC=BC,
NC=NC,
USE_SAFE_GATE=safe_gate,
)
w, u, qg, kg = recompute_w_u_fwd(
k=k,
v=v,
beta=beta,
A=Akk,
q=q if disable_recompute else None,
gk=gk,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
)
return w, u, qg, kg, Aqk, Akk

3.1 chunk_kda_fwd_intra_token_parallel

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
@triton.heuristics({
'IS_VARLEN': lambda args: args['cu_seqlens'] is not None,
})
@fla_cache_autotune(
configs=[
triton.Config({'BH': BH}, num_warps=num_warps)
for BH in [1, 2, 4, 8]
for num_warps in [1, 2, 4, 8]
],
key=["K", "H", "HV"],
**autotune_cache_kwargs,
)
@triton.jit(do_not_specialize=['T', 'N'])
def chunk_kda_fwd_kernel_intra_token_parallel(
q,
k,
g,
beta,
Aqk,
Akk,
scale,
cu_seqlens,
N,
T,
H: tl.constexpr,
HV: tl.constexpr,
K: tl.constexpr,
BT: tl.constexpr,
BC: tl.constexpr,
BH: tl.constexpr,
IS_VARLEN: tl.constexpr,
):
i_tg, i_hg = tl.program_id(0), tl.program_id(1)

if IS_VARLEN:
i_n = 0
left, right = 0, N

# Unrolled binary search (max B=2^32)
# We can limit iterations based on expected max batch size if needed
# 20 iterations covers B=1M, usually enough
for _ in range(20):
if left < right:
mid = (left + right) // 2
if i_tg < tl.load(cu_seqlens + mid + 1).to(tl.int32):
right = mid
else:
left = mid + 1
i_n = left

bos, eos = tl.load(cu_seqlens + i_n).to(tl.int32), tl.load(cu_seqlens + i_n + 1).to(tl.int32)
T = eos - bos
i_t = i_tg - bos
else:
bos = (i_tg // T) * T
i_t = i_tg % T # 当前 batch 内第几个 token

if i_t >= T:
return

i_c = i_t // BT # T 中第几个 chunk
i_s = (i_t % BT) // BC # chunk 内第几个 sun-chunk
i_tc = i_c * BT # chunk的起始位置
i_ts = i_tc + i_s * BC # sub-chunk的起始位置

G: tl.constexpr = HV // H

q += bos * H*K
k += bos * H*K
g += bos * HV*K
Aqk += bos * HV*BT
Akk += bos * HV*BC
beta += bos * HV

# BK = next_power_of_2(3) = 4。Triton 要求向量长度是 2 的幂,所以 K=3 要补到 4。
BK: tl.constexpr = triton.next_power_of_2(K)
o_hv = i_hg * BH + tl.arange(0, BH)
o_h = o_hv // G
o_k = tl.arange(0, BK)
# 例:o_hv=[0,1], HV=4 --> m_hv = [True, True]
# 例:o_k=[0,1,2,3], K=3 --> m_k = [True, True, True, False]
# m_hk 广播后形状 [BH=2, BK=4]
# [[True, True, True, False],
# [True, True, True, False]]
m_hv = o_hv < HV
m_k = o_k < K
m_hk = m_hv[:, None] & m_k[None, :]

# q/k: [B, T, H, K], manual load via mapped qk head index
# q 在当前 token 位置的布局是 [H=2, K=3],展平后:
# 位置: 0 1 2 3 4 5
# 对应: h0k0 h0k1 h0k2 h1k0 h1k1 h1k2
p_qk = o_h[:, None] * K + o_k[None, :]
b_q = tl.load(q + i_t * H * K + p_qk, mask=m_hk, other=0).to(tl.float32)
b_k = tl.load(k + i_t * H * K + p_qk, mask=m_hk, other=0).to(tl.float32)
# 基地址 = q + 5*2*3 = q + 30
# 加上 p_qk = [[0,1,2,3],[0,1,2,3]]
# 实际读取位置:
# [[q[30], q[31], q[32], q[33]],
# [q[30], q[31], q[32], q[33]]]
# mask 把第4列屏蔽(填0),最终:
# b_q = [[q[30], q[31], q[32], 0],
# [q[30], q[31], q[32], 0]]

# g: [B, T, HV, K], beta: [B, T, HV]
p_g = tl.make_block_ptr(g + i_t * HV * K, (HV, K), (K, 1), (i_hg * BH, 0), (BH, BK), (1, 0))
p_beta = tl.make_block_ptr(beta + i_t * HV, (HV,), (1,), (i_hg * BH,), (BH,), (0,))
b_g = tl.load(p_g, boundary_check=(0, 1)).to(tl.float32)
b_k = b_k * tl.load(p_beta, boundary_check=(0,)).to(tl.float32)[:, None]

for j in range(i_ts, min(i_t + 1, min(T, i_ts + BC))):
# 跟加载 token i 的方式完全一样,只是换成 token j 的位置。
b_kj = tl.load(k + j * H * K + p_qk, mask=m_hk, other=0).to(tl.float32)
p_gj = tl.make_block_ptr(g + j * HV * K, (HV, K), (K, 1), (i_hg * BH, 0), (BH, BK), (1, 0))
b_gj = tl.load(p_gj, boundary_check=(0, 1)).to(tl.float32)

b_kgj = tl.where(m_k[None, :], b_kj * exp2(b_g - b_gj), 0.0)
b_Aqk = tl.sum(b_q * b_kgj, axis=1) * scale
b_Akk = tl.sum(b_k * b_kgj, axis=1) * tl.where(j < i_t, 1.0, 0.0)

tl.store(Aqk + i_t * HV * BT + o_hv * BT + j % BT, b_Aqk.to(Aqk.dtype.element_ty), mask=m_hv)
tl.store(Akk + i_t * HV * BC + o_hv * BC + j - i_ts, b_Akk.to(Akk.dtype.element_ty), mask=m_hv)


def chunk_kda_fwd_intra_token_parallel(
q: torch.Tensor,
k: torch.Tensor,
gk: torch.Tensor,
beta: torch.Tensor,
Aqk: torch.Tensor,
Akk: torch.Tensor,
scale: float,
cu_seqlens: torch.LongTensor | None = None,
chunk_size: int = 64,
sub_chunk_size: int = 16,
) -> None:
"""
Token-parallel implementation: each token gets its own thread block.
Supports both fixed-length and variable-length sequences.
Reduces wasted computation on padding.

Writes directly to Aqk and Akk tensors (in-place).

Args:
q: [B, T, H, K]
k: [B, T, H, K]
gk: [B, T, HV, K] cumsum of gates (HV >= H for GVA)
beta: [B, T, HV]
Aqk: [B, T, HV, BT] output tensor to write to
Akk: [B, T, HV, BC] output tensor for diagonal blocks (fp32)
scale: attention scale
chunk_size: BT (default 64)
sub_chunk_size: BC (default 16)
"""
B, T, H, K, HV = *q.shape, gk.shape[2]
N = len(cu_seqlens) - 1 if cu_seqlens is not None else B
BT = chunk_size
BC = sub_chunk_size

def grid(meta): return (B * T, triton.cdiv(HV, meta['BH']))
chunk_kda_fwd_kernel_intra_token_parallel[grid](
q=q,
k=k,
g=gk,
beta=beta,
Aqk=Aqk,
Akk=Akk,
scale=scale,
cu_seqlens=cu_seqlens,
N=N,
T=T,
H=H,
HV=HV,
K=K,
BT=BT,
BC=BC,
)
return Aqk, Akk
1
2
3
4
5
6
7
8
9
10
11
12
axis1 (i_hg):     0         1
[hv0,hv1] [hv2,hv3]
axis0 (i_tg):
0 batch0-token0
1 batch0-token1
2 batch0-token2
...
7 batch0-token7
8 batch1-token0
9 batch1-token1
...
15 batch1-token7

block 视角

img

img

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Aqk 的非零模式(batch0, hv0 这一片)
chunk 内位置 j:
0 1 2 3 4 5 6 7
token 0: [ a00 . . . . . . . ]
token 1: [ a10 a11 . . . . . . ]
token 2: [ a20 a21 a22 . . . . . ]
token 3: [ a30 a31 a32 a33 . . . . ]
token 4: [ . . . . a44 . . . ]
token 5: [ . . . . a54 a55 . . ]
token 6: [ . . . . a64 a65 a66 . ]
token 7: [ . . . . a74 a75 a76 a77]

. = 零(未填充)
Akk 的非零模式(batch0, hv0 这一片,BC=4列)
token 0: [ 0 . . . ] 对角线=0
token 1: [ x 0 . . ]
token 2: [ x x 0 . ]
token 3: [ x x x 0 ]
token 4: [ 0 . . . ] 新sub-chunk,重新开始
token 5: [ x 0 . . ]
token 6: [ x x 0 . ]
token 7: [ x x x 0 ]

3.2 chunk_kda_fwd_kernel_inter_solve_fused

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
################################################################################
# Fused inter + solve_tril kernel: compute off-diagonal Akk and solve in one pass
################################################################################


@triton.heuristics({
'IS_VARLEN': lambda args: args['cu_seqlens'] is not None,
})
@fla_cache_autotune(
configs=[
triton.Config({'BK': BK}, num_warps=num_warps)
for BK in [32, 64]
for num_warps in [1, 2, 4]
],
key=["H", "HV", "K", "BT", "BC", "NC"],
**autotune_cache_kwargs,
)
@triton.jit(do_not_specialize=['T'])
def chunk_kda_fwd_kernel_inter_solve_fused(
q,
k,
g,
beta,
Aqk,# 输出: 注意力矩阵 [B, T, HV, BT]
Akkd,# 输入: Step1 算好的对角块逆 [B, T, HV, BC]
Akk,# 输出: 完整逆矩阵 [B, T, HV, BT]
scale,
cu_seqlens,
chunk_indices,
T,
H: tl.constexpr,
HV: tl.constexpr,
K: tl.constexpr,
BT: tl.constexpr,
BC: tl.constexpr,
NC: tl.constexpr,
BK: tl.constexpr,
IS_VARLEN: tl.constexpr,
USE_SAFE_GATE: tl.constexpr,
):
"""
Fused kernel: compute inter-subchunk Akk + solve_tril in one pass.
Prerequisite: token_parallel has already computed diagonal Akk blocks in Akkd.

This kernel:
1. Computes off-diagonal Aqk blocks -> writes to global
2. Computes off-diagonal Akk blocks -> keeps in registers
3. Loads diagonal Akk blocks from Akkd (fp32)
4. Does forward substitution on diagonals
5. Computes merged Akk_inv
6. Writes Akk_inv to Akk
"""
i_t, i_bh = tl.program_id(0), tl.program_id(1)
i_b, i_hv = i_bh // HV, i_bh % HV
i_h = i_hv // (HV // H)

if IS_VARLEN:
i_n, i_t = tl.load(chunk_indices + i_t * 2).to(tl.int32), tl.load(chunk_indices + i_t * 2 + 1).to(tl.int32)
bos, eos = tl.load(cu_seqlens + i_n).to(tl.int32), tl.load(cu_seqlens + i_n + 1).to(tl.int32)
T = eos - bos
else:
bos, eos = i_b * T, i_b * T + T

if i_t * BT >= T:
return

i_tc0 = i_t * BT
i_tc1 = i_t * BT + BC
i_tc2 = i_t * BT + 2 * BC
i_tc3 = i_t * BT + 3 * BC

q += (bos * H + i_h) * K
k += (bos * H + i_h) * K
g += (bos * HV + i_hv) * K
Aqk += (bos * HV + i_hv) * BT
Akk += (bos * HV + i_hv) * BT
Akkd += (bos * HV + i_hv) * BC

o_i = tl.arange(0, BC)
m_tc1 = (i_tc1 + o_i) < T
m_tc2 = (i_tc2 + o_i) < T
m_tc3 = (i_tc3 + o_i) < T

b_Aqk10 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk10 = tl.zeros([BC, BC], dtype=tl.float32)

b_Aqk20 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk20 = tl.zeros([BC, BC], dtype=tl.float32)
b_Aqk21 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk21 = tl.zeros([BC, BC], dtype=tl.float32)

b_Aqk30 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk30 = tl.zeros([BC, BC], dtype=tl.float32)
b_Aqk31 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk31 = tl.zeros([BC, BC], dtype=tl.float32)
b_Aqk32 = tl.zeros([BC, BC], dtype=tl.float32)
b_Akk32 = tl.zeros([BC, BC], dtype=tl.float32)

################################################################################
# off-diagonal blocks
################################################################################
# 计算非对角块 — 核心循环
# 外层循环:沿特征维 K 分块累加
for i_k in range(tl.cdiv(K, BK)):
o_k = i_k * BK + tl.arange(0, BK)
m_k = o_k < K

p_k0 = tl.make_block_ptr(k, (T, K), (H*K, 1), (i_tc0, i_k * BK), (BC, BK), (1, 0))
p_g0 = tl.make_block_ptr(g, (T, K), (HV*K, 1), (i_tc0, i_k * BK), (BC, BK), (1, 0))
b_k0 = tl.load(p_k0, boundary_check=(0, 1)).to(tl.float32)
b_g0 = tl.load(p_g0, boundary_check=(0, 1)).to(tl.float32)

if i_tc1 < T:
p_q1 = tl.make_block_ptr(q, (T, K), (H*K, 1), (i_tc1, i_k * BK), (BC, BK), (1, 0))
p_k1 = tl.make_block_ptr(k, (T, K), (H*K, 1), (i_tc1, i_k * BK), (BC, BK), (1, 0))
p_g1 = tl.make_block_ptr(g, (T, K), (HV*K, 1), (i_tc1, i_k * BK), (BC, BK), (1, 0))
# [BC, BK]
b_q1 = tl.load(p_q1, boundary_check=(0, 1)).to(tl.float32)
b_k1 = tl.load(p_k1, boundary_check=(0, 1)).to(tl.float32)
b_g1 = tl.load(p_g1, boundary_check=(0, 1)).to(tl.float32)
# [BK]
b_gn1 = tl.load(g + i_tc1 * HV*K + o_k, mask=m_k, other=0).to(tl.float32)
# [BC, BK]
b_gqn = tl.where(m_tc1[:, None], exp2(b_g1 - b_gn1[None, :]), 0)
# [BK, BC]
b_kgt = tl.trans(b_k0 * exp2(b_gn1[None, :] - b_g0))
# 2^(g[SC1_i] - g[SC0_j])
# = 2^(g[SC1_i] - gn1) × 2^(gn1 - g[SC0_j])
# ├──── b_gqn ────┤ ├──── 在b_kgt里 ────┤
# 防止溢出,保证数值稳定性
# [BC, BC]
b_Aqk10 += tl.dot(b_q1 * b_gqn, b_kgt)
b_Akk10 += tl.dot(b_k1 * b_gqn, b_kgt)

if NC >= 3 and i_tc2 < T:
p_q2 = tl.make_block_ptr(q, (T, K), (H*K, 1), (i_tc2, i_k * BK), (BC, BK), (1, 0))
p_k2 = tl.make_block_ptr(k, (T, K), (H*K, 1), (i_tc2, i_k * BK), (BC, BK), (1, 0))
p_g2 = tl.make_block_ptr(g, (T, K), (HV*K, 1), (i_tc2, i_k * BK), (BC, BK), (1, 0))
# [BC, BK]
b_q2 = tl.load(p_q2, boundary_check=(0, 1)).to(tl.float32)
b_k2 = tl.load(p_k2, boundary_check=(0, 1)).to(tl.float32)
b_g2 = tl.load(p_g2, boundary_check=(0, 1)).to(tl.float32)
# [BK]
b_gn2 = tl.load(g + i_tc2 * HV*K + o_k, mask=m_k, other=0).to(tl.float32)
# [BC, BK]
b_gqn2 = tl.where(m_tc2[:, None], exp2(b_g2 - b_gn2[None, :]), 0)
b_qg2 = b_q2 * b_gqn2
b_kg2 = b_k2 * b_gqn2
# [BK, BC]
b_kgt = tl.trans(b_k0 * exp2(b_gn2[None, :] - b_g0))
b_Aqk20 += tl.dot(b_qg2, b_kgt)
b_Akk20 += tl.dot(b_kg2, b_kgt)
# [BC, BC]
b_kgt = tl.trans(b_k1 * exp2(b_gn2[None, :] - b_g1))
# [BC, BC]
b_Aqk21 += tl.dot(b_qg2, b_kgt)
b_Akk21 += tl.dot(b_kg2, b_kgt)

if NC >= 4 and i_tc3 < T:
p_q3 = tl.make_block_ptr(q, (T, K), (H*K, 1), (i_tc3, i_k * BK), (BC, BK), (1, 0))
p_k3 = tl.make_block_ptr(k, (T, K), (H*K, 1), (i_tc3, i_k * BK), (BC, BK), (1, 0))
p_g3 = tl.make_block_ptr(g, (T, K), (HV*K, 1), (i_tc3, i_k * BK), (BC, BK), (1, 0))
# [BC, BK]
b_q3 = tl.load(p_q3, boundary_check=(0, 1)).to(tl.float32)
b_k3 = tl.load(p_k3, boundary_check=(0, 1)).to(tl.float32)
b_g3 = tl.load(p_g3, boundary_check=(0, 1)).to(tl.float32)
# [BK]
b_gn3 = tl.load(g + i_tc3 * HV*K + o_k, mask=m_k, other=0).to(tl.float32)
# [BC, BK]
b_gqn3 = tl.where(m_tc3[:, None], exp2(b_g3 - b_gn3[None, :]), 0)
b_qg3 = b_q3 * b_gqn3
b_kg3 = b_k3 * b_gqn3
# [BK, BC]
b_kgt = tl.trans(b_k0 * exp2(b_gn3[None, :] - b_g0))
# [BC, BC]
b_Aqk30 += tl.dot(b_qg3, b_kgt)
b_Akk30 += tl.dot(b_kg3, b_kgt)
# [BK, BC]
b_kgt = tl.trans(b_k1 * exp2(b_gn3[None, :] - b_g1))
# [BC, BC]
b_Aqk31 += tl.dot(b_qg3, b_kgt)
b_Akk31 += tl.dot(b_kg3, b_kgt)
# [BK, BC]
b_kgt = tl.trans(b_k2 * exp2(b_gn3[None, :] - b_g2))
# [BC, BC]
b_Aqk32 += tl.dot(b_qg3, b_kgt)
b_Akk32 += tl.dot(b_kg3, b_kgt)

################################################################################
# save off-diagonal Aqk blocks and prepare Akk
################################################################################
if i_tc1 < T:
p_Aqk10 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc1, 0), (BC, BC), (1, 0))
tl.store(p_Aqk10, (b_Aqk10 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))

p_b1 = tl.make_block_ptr(beta + bos * HV + i_hv, (T,), (HV,), (i_tc1,), (BC,), (0,))
b_b1 = tl.load(p_b1, boundary_check=(0,)).to(tl.float32)
b_Akk10 = b_Akk10 * b_b1[:, None]
if NC >= 3 and i_tc2 < T:
p_Aqk20 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc2, 0), (BC, BC), (1, 0))
p_Aqk21 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc2, BC), (BC, BC), (1, 0))
tl.store(p_Aqk20, (b_Aqk20 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Aqk21, (b_Aqk21 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))

p_b2 = tl.make_block_ptr(beta + bos * HV + i_hv, (T,), (HV,), (i_tc2,), (BC,), (0,))
b_b2 = tl.load(p_b2, boundary_check=(0,)).to(tl.float32)
b_Akk20 = b_Akk20 * b_b2[:, None]
b_Akk21 = b_Akk21 * b_b2[:, None]
if NC >= 4 and i_tc3 < T:
p_Aqk30 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc3, 0), (BC, BC), (1, 0))
p_Aqk31 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc3, BC), (BC, BC), (1, 0))
p_Aqk32 = tl.make_block_ptr(Aqk, (T, BT), (HV*BT, 1), (i_tc3, 2*BC), (BC, BC), (1, 0))
tl.store(p_Aqk30, (b_Aqk30 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Aqk31, (b_Aqk31 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Aqk32, (b_Aqk32 * scale).to(Aqk.dtype.element_ty), boundary_check=(0, 1))

p_b3 = tl.make_block_ptr(beta + bos * HV + i_hv, (T,), (HV,), (i_tc3,), (BC,), (0,))
b_b3 = tl.load(p_b3, boundary_check=(0,)).to(tl.float32)
b_Akk30 = b_Akk30 * b_b3[:, None]
b_Akk31 = b_Akk31 * b_b3[:, None]
b_Akk32 = b_Akk32 * b_b3[:, None]

p_Akk00 = tl.make_block_ptr(Akkd, (T, BC), (HV*BC, 1), (i_tc0, 0), (BC, BC), (1, 0))
p_Akk11 = tl.make_block_ptr(Akkd, (T, BC), (HV*BC, 1), (i_tc1, 0), (BC, BC), (1, 0))
b_Ai00 = tl.load(p_Akk00, boundary_check=(0, 1)).to(tl.float32)
b_Ai11 = tl.load(p_Akk11, boundary_check=(0, 1)).to(tl.float32)
if NC >= 3:
p_Akk22 = tl.make_block_ptr(Akkd, (T, BC), (HV*BC, 1), (i_tc2, 0), (BC, BC), (1, 0))
b_Ai22 = tl.load(p_Akk22, boundary_check=(0, 1)).to(tl.float32)
if NC >= 4:
p_Akk33 = tl.make_block_ptr(Akkd, (T, BC), (HV*BC, 1), (i_tc3, 0), (BC, BC), (1, 0))
b_Ai33 = tl.load(p_Akk33, boundary_check=(0, 1)).to(tl.float32)

################################################################################
# forward substitution on diagonals
################################################################################

if not USE_SAFE_GATE:
m_A = o_i[:, None] > o_i[None, :]
m_I = o_i[:, None] == o_i[None, :]

b_Ai00 = -tl.where(m_A, b_Ai00, 0)
b_Ai11 = -tl.where(m_A, b_Ai11, 0)
if NC >= 3:
b_Ai22 = -tl.where(m_A, b_Ai22, 0)
if NC >= 4:
b_Ai33 = -tl.where(m_A, b_Ai33, 0)

for i in range(2, min(BC, T - i_tc0)):
b_a00 = -tl.load(Akkd + (i_tc0 + i) * HV*BC + o_i)
b_a00 = tl.where(o_i < i, b_a00, 0.)
b_a00 += tl.sum(b_a00[:, None] * b_Ai00, 0)
b_Ai00 = tl.where((o_i == i)[:, None], b_a00, b_Ai00)
for i in range(BC + 2, min(2*BC, T - i_tc0)):
b_a11 = -tl.load(Akkd + (i_tc0 + i) * HV*BC + o_i)
b_a11 = tl.where(o_i < i - BC, b_a11, 0.)
b_a11 += tl.sum(b_a11[:, None] * b_Ai11, 0)
b_Ai11 = tl.where((o_i == i - BC)[:, None], b_a11, b_Ai11)
if NC >= 3:
for i in range(2*BC + 2, min(3*BC, T - i_tc0)):
b_a22 = -tl.load(Akkd + (i_tc0 + i) * HV*BC + o_i)
b_a22 = tl.where(o_i < i - 2*BC, b_a22, 0.)
b_a22 += tl.sum(b_a22[:, None] * b_Ai22, 0)
b_Ai22 = tl.where((o_i == i - 2*BC)[:, None], b_a22, b_Ai22)
if NC >= 4:
for i in range(3*BC + 2, min(4*BC, T - i_tc0)):
b_a33 = -tl.load(Akkd + (i_tc0 + i) * HV*BC + o_i)
b_a33 = tl.where(o_i < i - 3*BC, b_a33, 0.)
b_a33 += tl.sum(b_a33[:, None] * b_Ai33, 0)
b_Ai33 = tl.where((o_i == i - 3*BC)[:, None], b_a33, b_Ai33)

b_Ai00 += m_I
b_Ai11 += m_I
if NC >= 3:
b_Ai22 += m_I
if NC >= 4:
b_Ai33 += m_I

################################################################################
# compute merged inverse using off-diagonals
################################################################################

# we used tf32 to maintain matrix inverse's precision whenever possible.
b_Ai10 = -tl.dot(
tl.dot(b_Ai11, b_Akk10, input_precision=SOLVE_TRIL_DOT_PRECISION),
b_Ai00,
input_precision=SOLVE_TRIL_DOT_PRECISION
)

if NC >= 3:
b_Ai21 = -tl.dot(
tl.dot(b_Ai22, b_Akk21, input_precision=SOLVE_TRIL_DOT_PRECISION),
b_Ai11,
input_precision=SOLVE_TRIL_DOT_PRECISION
)
b_Ai20 = -tl.dot(
b_Ai22,
tl.dot(b_Akk20, b_Ai00, input_precision=SOLVE_TRIL_DOT_PRECISION) +
tl.dot(b_Akk21, b_Ai10, input_precision=SOLVE_TRIL_DOT_PRECISION),
input_precision=SOLVE_TRIL_DOT_PRECISION
)
if NC >= 4:
b_Ai32 = -tl.dot(
tl.dot(b_Ai33, b_Akk32, input_precision=SOLVE_TRIL_DOT_PRECISION),
b_Ai22,
input_precision=SOLVE_TRIL_DOT_PRECISION
)
b_Ai31 = -tl.dot(
b_Ai33,
tl.dot(b_Akk31, b_Ai11, input_precision=SOLVE_TRIL_DOT_PRECISION) +
tl.dot(b_Akk32, b_Ai21, input_precision=SOLVE_TRIL_DOT_PRECISION),
input_precision=SOLVE_TRIL_DOT_PRECISION
)
b_Ai30 = -tl.dot(
b_Ai33,
tl.dot(b_Akk30, b_Ai00, input_precision=SOLVE_TRIL_DOT_PRECISION) +
tl.dot(b_Akk31, b_Ai10, input_precision=SOLVE_TRIL_DOT_PRECISION) +
tl.dot(b_Akk32, b_Ai20, input_precision=SOLVE_TRIL_DOT_PRECISION),
input_precision=SOLVE_TRIL_DOT_PRECISION
)

################################################################################
# store full Akk_inv to Akk
################################################################################

p_Akk00 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc0, 0), (BC, BC), (1, 0))
p_Akk10 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc1, 0), (BC, BC), (1, 0))
p_Akk11 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc1, BC), (BC, BC), (1, 0))

tl.store(p_Akk00, b_Ai00.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk10, b_Ai10.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk11, b_Ai11.to(Akk.dtype.element_ty), boundary_check=(0, 1))
if NC >= 3:
p_Akk20 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc2, 0), (BC, BC), (1, 0))
p_Akk21 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc2, BC), (BC, BC), (1, 0))
p_Akk22 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc2, 2*BC), (BC, BC), (1, 0))
tl.store(p_Akk20, b_Ai20.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk21, b_Ai21.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk22, b_Ai22.to(Akk.dtype.element_ty), boundary_check=(0, 1))
if NC >= 4:
p_Akk30 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc3, 0), (BC, BC), (1, 0))
p_Akk31 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc3, BC), (BC, BC), (1, 0))
p_Akk32 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc3, 2*BC), (BC, BC), (1, 0))
p_Akk33 = tl.make_block_ptr(Akk, (T, BT), (HV*BT, 1), (i_tc3, 3*BC), (BC, BC), (1, 0))
tl.store(p_Akk30, b_Ai30.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk31, b_Ai31.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk32, b_Ai32.to(Akk.dtype.element_ty), boundary_check=(0, 1))
tl.store(p_Akk33, b_Ai33.to(Akk.dtype.element_ty), boundary_check=(0, 1))


grid = (NT, B * HV)
chunk_kda_fwd_kernel_inter_solve_fused[grid](
q=q,
k=k,
g=gk,
beta=beta,
Aqk=Aqk,
Akkd=Akkd,
Akk=Akk,
scale=scale,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
T=T,
H=H,
HV=HV,
K=K,
BT=BT,
BC=BC,
NC=NC,
USE_SAFE_GATE=safe_gate,
)

每个program处理一个(chunk, batch, hv_head)组合

把一个 chunk 切成了 四个 sub-chunk

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
64×64 矩阵的分块结构 (行=query的SC, 列=key的SC):

SC0列 SC1列 SC2列 SC3列
┌───────────┬───────────┬───────────┬───────────┐
SC0行 │ (对角块) │ 0 │ 0 │ 0 │
├───────────┼───────────┼───────────┼───────────┤
SC1行 │ Aqk10 │ (对角块) │ 0 │ 0 │
│ Akk10 │ │ │ │
├───────────┼───────────┼───────────┼───────────┤
SC2行 │ Aqk20 │ Aqk21 │ (对角块) │ 0 │
│ Akk20 │ Akk21 │ │ │
├───────────┼───────────┼───────────┼───────────┤
SC3行 │ Aqk30 │ Aqk31 │ Aqk32 │ (对角块) │
│ Akk30 │ Akk31 │ Akk32 │ │
└───────────┴───────────┴───────────┴───────────┘

这些非对角块就是要计算的目标

block视角

img

整体流程如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
═══════════════════════════════════════════════════════════════
一个线程块的完整工作流
═══════════════════════════════════════════════════════════════

输入: 一个 chunk (64 token), 一个 batch, 一个 head
q[64,K], k[64,K], g[64,K], beta[64], Akkd[64,16](来自Step1)

Step A: 定位 (Line 76-100)
─────────────────────────────────
"我是第 i_t 个 chunk, batch i_b, head i_hv"
计算指针偏移

Step B: 沿 K 维度循环累加非对角块 (Line 125-204)
─────────────────────────────────
for 每 BK 维特征:
加载 SC0 的 k, g
├─→ 计算 SC1→SC0: Aqk10, Akk10
├─→ 计算 SC2→SC0: Aqk20, Akk20
│ 计算 SC2→SC1: Aqk21, Akk21
└─→ 计算 SC3→SC0: Aqk30, Akk30
计算 SC3→SC1: Aqk31, Akk31
计算 SC3→SC2: Aqk32, Akk32

Step C: 后处理 (Line 209-249)
─────────────────────────────────
Aqk 非对角块 × scale → 写入显存
Akk 非对角块 × beta → 留在寄存器
加载 Step1 算好的对角块 Akkd → b_Ai00..33

Step D: 对角块前向代入求逆 (Line 255-294)
─────────────────────────────────
(仅 USE_SAFE_GATE=False 时)
对 b_Ai00..33 逐行做前向代入
得到 (I - L_diag)^{-1}

Step E: 分块逆矩阵合并 (Line 301-337)
─────────────────────────────────
Ai10 = -(Ai11 · L10 · Ai00)
Ai21 = -(Ai22 · L21 · Ai11)
Ai20 = -(Ai22 · (L20·Ai00 + L21·Ai10))
Ai32 = -(Ai33 · L32 · Ai22)
Ai31 = -(Ai33 · (L31·Ai11 + L32·Ai21))
Ai30 = -(Ai33 · (L30·Ai00 + L31·Ai10 + L32·Ai20))

Step F: 写入最终结果 (Line 343-365)
─────────────────────────────────
把 10 个 16×16 块拼成 64×64 矩阵写入 Akk

输出: Aqk[64,64] (注意力分数), Akk[64,64] (逆矩阵)
═══════════════════════════════════════════════════════════════

3.2.1 为什么要分对角和非对角块计算

假设 BT=16, BC=4, 所以 NC=BT/BC=4 个 sub-chunk。

一个 chunk 内有 16 个 token(编号 0~15),因果注意力矩阵是 16x16 的下三角。

现在把这个 16x16 矩阵按 sub-chunk 分成 4x4 的块矩阵,每块 4x4:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
token:  0  1  2  3 | 4  5  6  7 | 8  9 10 11 |12 13 14 15
sc0 | sc1 | sc2 | sc3

0 [ x . . . | . . . . | . . . . | . . . . ]
1 [ x x . . | . . . . | . . . . | . . . . ]
2 [ x x x . | . . . . | . . . . | . . . . ]
3 [ x x x x | . . . . | . . . . | . . . . ]
+------------+------------+------------+-------------
4 [ x x x x | x . . . | . . . . | . . . . ]
5 [ x x x x | x x . . | . . . . | . . . . ]
6 [ x x x x | x x x . | . . . . | . . . . ]
7 [ x x x x | x x x x | . . . . | . . . . ]
+------------+------------+------------+-------------
8 [ x x x x | x x x x | x . . . | . . . . ]
9 [ x x x x | x x x x | x x . . | . . . . ]
10 [ x x x x | x x x x | x x x . | . . . . ]
11 [ x x x x | x x x x | x x x x | . . . . ]
+------------+------------+------------+-------------
12 [ x x x x | x x x x | x x x x | x . . . ]
13 [ x x x x | x x x x | x x x x | x x . . ]
14 [ x x x x | x x x x | x x x x | x x x . ]
15 [ x x x x | x x x x | x x x x | x x x x ]

把它看成 4x4 的块矩阵:

1
2
3
4
5
      sc0    sc1    sc2    sc3
sc0 [ D00 . . . ]
sc1 [ O10 D11 . . ]
sc2 [ O20 O21 D22 . ]
sc3 [ O30 O31 O32 D33 ]

D00, D11, D22, D33 是"对角块":它们在主对角线上,每个都是 4x4 下三角矩阵。 O10, O20, O21, O30, O31, O32 是"非对角块":它们在主对角线下方,每个都是 4x4 全满矩阵。

用颜色标记(D=对角,O=非对角):

1
2
3
4
5
      sc0    sc1    sc2    sc3
sc0 [ DDDD .... .... .... ]
sc1 [ OOOO DDDD .... .... ]
sc2 [ OOOO OOOO DDDD .... ]
sc3 [ OOOO OOOO OOOO DDDD ]

对角块内部啥样?

拿 D11(sc1 对应的对角块)来看,它对应 token 4~7 对 token 4~7:

1
2
3
4
5
      j=4  j=5  j=6  j=7
i=4 [ x . . . ]
i=5 [ x x . . ]
i=6 [ x x x . ]
i=7 [ x x x x ]

这是下三角!因为因果性:token i 只能看到 j <= i。

每一行的有效元素个数不同:

  • i=4: 1个
  • i=5: 2个
  • i=6: 3个
  • i=7: 4个

这就是"不规则"的。

非对角块内部长什么样?

拿 O20(sc2 对应行,sc0 对应列)来看,它对应 token 8~11 对 token 0~3:

1
2
3
4
5
       j=0  j=1  j=2  j=3
i=8 [ x x x x ]
i=9 [ x x x x ]
i=10 [ x x x x ]
i=11 [ x x x x ]

全满!因为 sc2 的每个 token(8~11)都在 sc0 的每个 token(0~3)之后。 因果性自动满足,不需要 mask。

两种块的计算方式不同

对角块:逐行串行循环

token-parallel kernel 对每个 token 跑一个 program。 以 token 6(在 D11 中)为例:

1
2
3
4
      j=4  j=5  j=6  j=7
i=6 [ x x x . ]
^
这一行有3个元素要算

kernel 的循环:

1
2
3
j=4: Aqk[6,4] = <q6, k4 * decay> * scale     -- 算1个标量
j=5: Aqk[6,5] = <q6, k5 * decay> * scale -- 算1个标量
j=6: Aqk[6,6] = <q6, k6 * decay> * scale -- 算1个标量(对角)

每轮循环就是一次向量内积(沿 K 维 sum),产出一个标量。

而 token 4 只需要循环1次,token 7 需要循环4次。 每个 token 的工作量不同,但因为是 token 并行,所以每个 token 独立,没问题。

关键:这里没法用矩阵乘。因为:

  • 每个 program 只处理一行(一个 token)
  • 产出的是一行里的若干个标量
  • 矩阵乘需要 [M,K] x [K,N] 产出 [M,N],但这里 M=1(单行),退化成向量-矩阵乘,对 tensor core 来说效率很低

非对角块:天然的矩阵乘法

O20 对应 sc2(token 8~11)对 sc0(token 0~3),是一个 4x4 全满矩阵。

计算公式:

1
O20[i, j] = sum_k( q[i,k] * k[j,k] * 2^(g[i,k] - g[j,k]) )

可以改写成矩阵乘的形式:

1
O20 = (q_sc2 * gate_factor) @ (k_sc0 * gate_factor).T

即:

  • 左矩阵: [4, K] – sc2 的 4 个 token 的 gated query
  • 右矩阵: [K, 4] – sc0 的 4 个 token 的 gated key 的转置
  • 结果: [4, 4] – 完整的非对角块

对应代码(inter_solve_fused kernel 第 149 行):

1
b_Aqk10 += tl.dot(b_q1 * b_gqn, b_kgt)

这就是 [BC, BK] x [BK, BC] = [BC, BC]

tl.dot 调用 GPU 的 tensor core,这是 GPU 最擅长的操作。

对比

1
2
3
4
5
6
7
8
9
10
11
12
对角块 D11:                       非对角块 O20:

j=4 j=5 j=6 j=7 j=0 j=1 j=2 j=3
i=4 [ x . . . ] i=8 [ x x x x ]
i=5 [ x x . . ] i=9 [ x x x x ]
i=6 [ x x x . ] i=10 [ x x x x ]
i=7 [ x x x x ] i=11 [ x x x x ]

形状:下三角,不规则 形状:全满矩形
计算:逐 j 循环,标量内积 计算:一次矩阵乘 tl.dot
并行:每个 token 一个 program 并行:每个 chunk 一个 program
一次算完整个 4x4 块

如果强行把非对角块也用逐 j 循环来算:

  • O20 有 4x4=16 个元素
  • 每个 token 循环 4 次
  • 4 个 token 分别独立循环
  • 总共 16 次标量内积

如果用矩阵乘:

  • 1 次 tl.dot 就算完整个 4x4
  • tensor core 一次算 16 个结果
  • 快十几倍到几十倍

反过来,如果强行对角块也用矩阵乘:

  • 先算全满 4x4 矩阵,再 mask 掉上三角
  • 浪费了约一半计算
  • 而且对角块后面还要做 forward substitution(见下文),这一步无论如何是串行的

3.3 recompute_w_u_fwd

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def recompute_w_u_fwd(
k: torch.Tensor,
v: torch.Tensor,
beta: torch.Tensor,
A: torch.Tensor,
gk: torch.Tensor,
q: torch.Tensor | None = None,
cu_seqlens: torch.LongTensor | None = None,
chunk_indices: torch.LongTensor | None = None,
) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor | None, torch.Tensor | None]:
B, T, H, K, V = *k.shape, v.shape[-1]
HV = v.shape[2]
BT = A.shape[-1]
BK = 64
BV = 64

if chunk_indices is None and cu_seqlens is not None:
chunk_indices = prepare_chunk_indices(cu_seqlens, BT)
NT = triton.cdiv(T, BT) if cu_seqlens is None else len(chunk_indices)

w = torch.empty(B, T, HV, K, device=k.device, dtype=k.dtype)
u = torch.empty_like(v)
qg = torch.empty(B, T, HV, K, device=k.device, dtype=k.dtype) if q is not None else None
kg = torch.empty(B, T, HV, K, device=k.device, dtype=k.dtype)
recompute_w_u_fwd_kda_kernel[(NT, B*HV)](
q=q,
k=k,
qg=qg,
kg=kg,
v=v,
beta=beta,
w=w,
u=u,
A=A,
gk=gk,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
T=T,
H=H,
HV=HV,
K=K,
V=V,
BT=BT,
BK=BK,
BV=BV,
)
return w, u, qg, kg


@triton.heuristics({
'STORE_QG': lambda args: args['qg'] is not None,
'STORE_KG': lambda args: args['kg'] is not None,
'IS_VARLEN': lambda args: args['cu_seqlens'] is not None,
})
@fla_cache_autotune(
configs=[
triton.Config({}, num_warps=num_warps, num_stages=num_stages)
for num_warps in [2, 4, 8]
for num_stages in [2, 3, 4]
],
key=['H', 'HV', 'K', 'V', 'BT', 'BK', 'BV', 'IS_VARLEN'],
**autotune_cache_kwargs,
)
@triton.jit(do_not_specialize=['T'])
def recompute_w_u_fwd_kda_kernel(
q,
k,
qg,
kg,
v,
beta,
w,
u,
A,
gk,
cu_seqlens,
chunk_indices,
T,
H: tl.constexpr,
HV: tl.constexpr,
K: tl.constexpr,
V: tl.constexpr,
BT: tl.constexpr,
BK: tl.constexpr,
BV: tl.constexpr,
STORE_QG: tl.constexpr,
STORE_KG: tl.constexpr,
IS_VARLEN: tl.constexpr,
):
i_t, i_bh = tl.program_id(0), tl.program_id(1)
i_b, i_hv = i_bh // HV, i_bh % HV
i_h = i_hv // (HV // H)
if IS_VARLEN:
i_n, i_t = tl.load(chunk_indices + i_t * 2).to(tl.int32), tl.load(chunk_indices + i_t * 2 + 1).to(tl.int32)
bos, eos = tl.load(cu_seqlens + i_n).to(tl.int32), tl.load(cu_seqlens + i_n + 1).to(tl.int32)
T = eos - bos
else:
bos, eos = i_b * T, i_b * T + T

k += (bos * H + i_h) * K
v += (bos * HV + i_hv) * V
u += (bos * HV + i_hv) * V
w += (bos * HV + i_hv) * K
gk += (bos * HV + i_hv) * K
beta += bos * HV + i_hv
A += (bos * HV + i_hv) * BT
if STORE_QG:
q += (bos * H + i_h) * K
qg += (bos * HV + i_hv) * K
if STORE_KG:
kg += (bos * HV + i_hv) * K

p_b = tl.make_block_ptr(beta, (T,), (HV,), (i_t * BT,), (BT,), (0,))
b_b = tl.load(p_b, boundary_check=(0,))

p_A = tl.make_block_ptr(A, (T, BT), (HV*BT, 1), (i_t * BT, 0), (BT, BT), (1, 0))
b_A = tl.load(p_A, boundary_check=(0, 1))

for i_v in range(tl.cdiv(V, BV)):
p_v = tl.make_block_ptr(v, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
p_u = tl.make_block_ptr(u, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
b_v = tl.load(p_v, boundary_check=(0, 1))
b_vb = (b_v * b_b[:, None]).to(b_v.dtype)
b_u = tl.dot(b_A, b_vb)
tl.store(p_u, b_u.to(p_u.dtype.element_ty), boundary_check=(0, 1))

for i_k in range(tl.cdiv(K, BK)):
p_w = tl.make_block_ptr(w, (T, K), (HV*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
p_k = tl.make_block_ptr(k, (T, K), (H*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
b_k = tl.load(p_k, boundary_check=(0, 1))
b_kb = b_k * b_b[:, None]

p_gk = tl.make_block_ptr(gk, (T, K), (HV*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
b_gk = tl.load(p_gk, boundary_check=(0, 1)).to(tl.float32)
b_kb *= exp2(b_gk)
if STORE_QG:
p_q = tl.make_block_ptr(q, (T, K), (H*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
p_qg = tl.make_block_ptr(qg, (T, K), (HV*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
b_q = tl.load(p_q, boundary_check=(0, 1))
b_qg = b_q * exp2(b_gk)
tl.store(p_qg, b_qg.to(p_qg.dtype.element_ty), boundary_check=(0, 1))
if STORE_KG:
last_idx = min(i_t * BT + BT, T) - 1
o_k = i_k * BK + tl.arange(0, BK)
m_k = o_k < K
b_gn = tl.load(gk + last_idx * HV*K + o_k, mask=m_k, other=0.).to(tl.float32)
b_kg = b_k * tl.where((i_t * BT + tl.arange(0, BT) < T)[:, None], exp2(b_gn[None, :] - b_gk), 0)
p_kg = tl.make_block_ptr(kg, (T, K), (HV*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
tl.store(p_kg, b_kg.to(p_kg.dtype.element_ty), boundary_check=(0, 1))

b_w = tl.dot(b_A, b_kb.to(b_k.dtype))
tl.store(p_w, b_w.to(p_w.dtype.element_ty), boundary_check=(0, 1))

每个program处理一个(chunk, batch, hv_head)组合

这个kernel计算4样东西:

输出 公式 含义
u A @ (v × beta) 修正后的value
w A @ (k × beta × exp2(gk)) 修正后的key(用于更新隐状态)
qg q × exp2(gk) gate衰减后的query(可选,仅disable_recompute时存)
kg k × exp2(gn - gk) chunk末尾归一化后的key(用于chunk间递推)

其中 A 就是前一步求出的 Akk_inv,形状 [B, T, HV, BT]

block 视角

img

计算 u = A @ (v × beta),V维度分块, BV=64

  • 具体例子(BT=4, BV=2)

  • 1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    b_A = [[1,    0,    0,    0],
    [a10, 1, 0, 0],
    [a20, a21, 1, 0],
    [a30, a31, a32, 1]]

    b_vb = [[v00×b0, v01×b0],
    [v10×b1, v11×b1],
    [v20×b2, v21×b2],
    [v30×b3, v31×b3]]

    b_u = b_A @ b_vb

    u[0] = v0×b0
    u[1] = a10×v0×b0 + v1×b1
    u[2] = a20×v0×b0 + a21×v1×b1 + v2×b2
    u[3] = a30×v0×b0 + a31×v1×b1 + a32×v2×b2 + v3×b3

计算 w = A @ (k × beta × exp2(gk)),K维度分块, BK=64。和计算 u 的逻辑完全对称

4 chunk_gated_delta_rule_fwd_h

代码链接:flash-linear-attention/fla/ops/common/chunk_delta_h.py at main · fla-org/flash-linear-attention

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
@dispatch('common')
def chunk_gated_delta_rule_fwd_h(
k: torch.Tensor,
w: torch.Tensor,
u: torch.Tensor,
g: torch.Tensor | None = None,
gk: torch.Tensor | None = None,
initial_state: torch.Tensor | None = None,
output_final_state: bool = False,
chunk_size: int = 64,
save_new_value: bool = True,
state_v_first: bool = False,
cu_seqlens: torch.LongTensor | None = None,
cu_seqlens_cpu: torch.LongTensor | None = None,
chunk_indices: torch.LongTensor | None = None,
) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor | None]:
B, T, H, K, V, HV = *k.shape, u.shape[-1], u.shape[2]
BT = chunk_size

if chunk_indices is None and cu_seqlens is not None:
chunk_indices = prepare_chunk_indices(cu_seqlens, chunk_size)
# N: the actual number of sequences in the batch with either equal or variable lengths
if cu_seqlens is None:
N, NT, chunk_offsets = B, triton.cdiv(T, BT), None
else:
N, NT, chunk_offsets = len(cu_seqlens) - 1, len(chunk_indices), prepare_chunk_offsets(cu_seqlens, BT)
assert K <= 256, "current kernel does not support head dimension larger than 256."

if state_v_first:
h = k.new_empty(B, NT, HV, V, K)
final_state = k.new_zeros(N, HV, V, K, dtype=torch.float32) if output_final_state else None
else:
h = k.new_empty(B, NT, HV, K, V)
final_state = k.new_zeros(N, HV, K, V, dtype=torch.float32) if output_final_state else None

v_new = torch.empty_like(u) if save_new_value else None
def grid(meta): return (triton.cdiv(V, meta['BV']), N*HV)
chunk_gated_delta_rule_fwd_kernel_h_blockdim64[grid](
k=k,
v=u,
w=w,
v_new=v_new,
g=g,
gk=gk,
h=h,
h0=initial_state,
ht=final_state,
cu_seqlens=cu_seqlens,
chunk_offsets=chunk_offsets,
T=T,
H=H,
HV=HV,
K=K,
V=V,
BT=BT,
STATE_V_FIRST=state_v_first,
)
return h, v_new, final_state

@triton.heuristics({
'USE_G': lambda args: args['g'] is not None,
'USE_GK': lambda args: args['gk'] is not None,
'USE_INITIAL_STATE': lambda args: args['h0'] is not None,
'STORE_FINAL_STATE': lambda args: args['ht'] is not None,
'SAVE_NEW_VALUE': lambda args: args['v_new'] is not None,
'IS_VARLEN': lambda args: args['cu_seqlens'] is not None,
})
@fla_cache_autotune(
configs=[
triton.Config({'BV': BV}, num_warps=num_warps, num_stages=num_stages)
for num_warps in [2, 4]
for num_stages in ([2, 3, 4] if check_shared_mem('ampere') else [2, 1])
for BV in ([32, 64] if check_shared_mem('ada') else [32])
],
key=['H', 'HV', 'K', 'V', 'BT', 'STATE_V_FIRST'],
use_cuda_graph=USE_CUDA_GRAPH,
**autotune_cache_kwargs,
)
@triton.jit(do_not_specialize=['T'])
def chunk_gated_delta_rule_fwd_kernel_h_blockdim64(
k,
v,
w,
v_new,
g,
gk,
h,
h0,
ht,
cu_seqlens,
chunk_offsets,
T,
H: tl.constexpr,
HV: tl.constexpr,
K: tl.constexpr,
V: tl.constexpr,
BT: tl.constexpr,
BV: tl.constexpr,
USE_G: tl.constexpr,
USE_GK: tl.constexpr,
USE_INITIAL_STATE: tl.constexpr,
STORE_FINAL_STATE: tl.constexpr,
SAVE_NEW_VALUE: tl.constexpr,
STATE_V_FIRST: tl.constexpr,
IS_VARLEN: tl.constexpr,
):
i_v, i_nh = tl.program_id(0), tl.program_id(1)
i_n, i_h = i_nh // HV, i_nh % HV
if IS_VARLEN:
bos, eos = tl.load(cu_seqlens + i_n).to(tl.int32), tl.load(cu_seqlens + i_n + 1).to(tl.int32)
T = eos - bos
NT = tl.cdiv(T, BT)
boh = tl.load(chunk_offsets + i_n).to(tl.int32)
else:
bos, eos = i_n * T, i_n * T + T
NT = tl.cdiv(T, BT)
boh = i_n * NT

if STATE_V_FIRST:
b_h1 = tl.zeros([BV, 64], dtype=tl.float32)
if K > 64:
b_h2 = tl.zeros([BV, 64], dtype=tl.float32)
if K > 128:
b_h3 = tl.zeros([BV, 64], dtype=tl.float32)
if K > 192:
b_h4 = tl.zeros([BV, 64], dtype=tl.float32)
else:
b_h1 = tl.zeros([64, BV], dtype=tl.float32)
if K > 64:
b_h2 = tl.zeros([64, BV], dtype=tl.float32)
if K > 128:
b_h3 = tl.zeros([64, BV], dtype=tl.float32)
if K > 192:
b_h4 = tl.zeros([64, BV], dtype=tl.float32)

# calculate offset
h += (boh * HV + i_h).to(tl.int64) * K*V
v += (bos * HV + i_h).to(tl.int64) * V
k += (bos * H + i_h // (HV // H)).to(tl.int64) * K
w += (bos * HV + i_h).to(tl.int64) * K
if SAVE_NEW_VALUE:
v_new += (bos * HV + i_h).to(tl.int64) * V

if USE_INITIAL_STATE:
h0 = h0 + i_nh * K*V
if STORE_FINAL_STATE:
ht = ht + i_nh * K*V

# load initial state
if USE_INITIAL_STATE:
if STATE_V_FIRST:
p_h0_1 = tl.make_block_ptr(h0, (V, K), (K, 1), (i_v * BV, 0), (BV, 64), (1, 0))
else:
p_h0_1 = tl.make_block_ptr(h0, (K, V), (V, 1), (0, i_v * BV), (64, BV), (1, 0))
b_h1 += tl.load(p_h0_1, boundary_check=(0, 1)).to(tl.float32)
if K > 64:
if STATE_V_FIRST:
p_h0_2 = tl.make_block_ptr(h0, (V, K), (K, 1), (i_v * BV, 64), (BV, 64), (1, 0))
else:
p_h0_2 = tl.make_block_ptr(h0, (K, V), (V, 1), (64, i_v * BV), (64, BV), (1, 0))
b_h2 += tl.load(p_h0_2, boundary_check=(0, 1)).to(tl.float32)
if K > 128:
if STATE_V_FIRST:
p_h0_3 = tl.make_block_ptr(h0, (V, K), (K, 1), (i_v * BV, 128), (BV, 64), (1, 0))
else:
p_h0_3 = tl.make_block_ptr(h0, (K, V), (V, 1), (128, i_v * BV), (64, BV), (1, 0))
b_h3 += tl.load(p_h0_3, boundary_check=(0, 1)).to(tl.float32)
if K > 192:
if STATE_V_FIRST:
p_h0_4 = tl.make_block_ptr(h0, (V, K), (K, 1), (i_v * BV, 192), (BV, 64), (1, 0))
else:
p_h0_4 = tl.make_block_ptr(h0, (K, V), (V, 1), (192, i_v * BV), (64, BV), (1, 0))
b_h4 += tl.load(p_h0_4, boundary_check=(0, 1)).to(tl.float32)

# main recurrence
for i_t in range(NT):
i_t_int64 = i_t.to(tl.int64)
if STATE_V_FIRST:
p_h1 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (V, K), (K, 1), (i_v * BV, 0), (BV, 64), (1, 0))
else:
p_h1 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (K, V), (V, 1), (0, i_v * BV), (64, BV), (1, 0))
tl.store(p_h1, b_h1.to(p_h1.dtype.element_ty), boundary_check=(0, 1))
if K > 64:
if STATE_V_FIRST:
p_h2 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (V, K), (K, 1), (i_v * BV, 64), (BV, 64), (1, 0))
else:
p_h2 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (K, V), (V, 1), (64, i_v * BV), (64, BV), (1, 0))
tl.store(p_h2, b_h2.to(p_h2.dtype.element_ty), boundary_check=(0, 1))
if K > 128:
if STATE_V_FIRST:
p_h3 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (V, K), (K, 1), (i_v * BV, 128), (BV, 64), (1, 0))
else:
p_h3 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (K, V), (V, 1), (128, i_v * BV), (64, BV), (1, 0))
tl.store(p_h3, b_h3.to(p_h3.dtype.element_ty), boundary_check=(0, 1))
if K > 192:
if STATE_V_FIRST:
p_h4 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (V, K), (K, 1), (i_v * BV, 192), (BV, 64), (1, 0))
else:
p_h4 = tl.make_block_ptr(h + i_t_int64 * HV*K*V, (K, V), (V, 1), (192, i_v * BV), (64, BV), (1, 0))
tl.store(p_h4, b_h4.to(p_h4.dtype.element_ty), boundary_check=(0, 1))

p_w = tl.make_block_ptr(w, (T, K), (HV*K, 1), (i_t * BT, 0), (BT, 64), (1, 0))
b_w = tl.load(p_w, boundary_check=(0, 1))
if STATE_V_FIRST:
b_v = tl.dot(b_w, tl.trans(b_h1).to(b_w.dtype))
else:
b_v = tl.dot(b_w, b_h1.to(b_w.dtype))
if K > 64:
p_w = tl.make_block_ptr(w, (T, K), (HV*K, 1), (i_t * BT, 64), (BT, 64), (1, 0))
b_w = tl.load(p_w, boundary_check=(0, 1))
if STATE_V_FIRST:
b_v += tl.dot(b_w, tl.trans(b_h2).to(b_w.dtype))
else:
b_v += tl.dot(b_w, b_h2.to(b_w.dtype))
if K > 128:
p_w = tl.make_block_ptr(w, (T, K), (HV*K, 1), (i_t * BT, 128), (BT, 64), (1, 0))
b_w = tl.load(p_w, boundary_check=(0, 1))
if STATE_V_FIRST:
b_v += tl.dot(b_w, tl.trans(b_h3).to(b_w.dtype))
else:
b_v += tl.dot(b_w, b_h3.to(b_w.dtype))
if K > 192:
p_w = tl.make_block_ptr(w, (T, K), (HV*K, 1), (i_t * BT, 192), (BT, 64), (1, 0))
b_w = tl.load(p_w, boundary_check=(0, 1))
if STATE_V_FIRST:
b_v += tl.dot(b_w, tl.trans(b_h4).to(b_w.dtype))
else:
b_v += tl.dot(b_w, b_h4.to(b_w.dtype))
p_v = tl.make_block_ptr(v, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
b_v = tl.load(p_v, boundary_check=(0, 1)) - b_v

if SAVE_NEW_VALUE:
p_v = tl.make_block_ptr(v_new, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
tl.store(p_v, b_v.to(p_v.dtype.element_ty), boundary_check=(0, 1))

last_idx = min((i_t + 1) * BT, T) - 1
if USE_G:
m_t = (i_t * BT + tl.arange(0, BT)) < T
b_g_last = tl.load(g + (bos * HV + last_idx * HV + i_h).to(tl.int64)).to(tl.float32)
p_g = tl.make_block_ptr(g + (bos * HV + i_h).to(tl.int64), (T,), (HV,), (i_t * BT,), (BT,), (0,))
b_g = tl.load(p_g, boundary_check=(0,)).to(tl.float32)
b_v = b_v * tl.where(m_t, exp2(b_g_last - b_g), 0)[:, None]
b_g_last = exp2(b_g_last)
b_h1 *= b_g_last
if K > 64:
b_h2 *= b_g_last
if K > 128:
b_h3 *= b_g_last
if K > 192:
b_h4 *= b_g_last

if USE_GK:
o_k1 = tl.arange(0, 64)
b_gk_last1 = tl.load(gk + (bos + last_idx) * HV*K + i_h * K + o_k1, mask=(o_k1 < K), other=0.).to(tl.float32)
if STATE_V_FIRST:
b_h1 *= exp2(b_gk_last1)[None, :]
else:
b_h1 *= exp2(b_gk_last1)[:, None]
if K > 64:
o_k2 = 64 + o_k1
b_gk_last2 = tl.load(gk + (bos + last_idx) * HV*K + i_h * K + o_k2, mask=(o_k2 < K), other=0.).to(tl.float32)
if STATE_V_FIRST:
b_h2 *= exp2(b_gk_last2)[None, :]
else:
b_h2 *= exp2(b_gk_last2)[:, None]
if K > 128:
o_k3 = 128 + o_k1
b_gk_last3 = tl.load(gk + (bos + last_idx) * HV*K + i_h * K + o_k3, mask=(o_k3 < K), other=0.).to(tl.float32)
if STATE_V_FIRST:
b_h3 *= exp2(b_gk_last3)[None, :]
else:
b_h3 *= exp2(b_gk_last3)[:, None]
if K > 192:
o_k4 = 192 + o_k1
b_gk_last4 = tl.load(gk + (bos + last_idx) * HV*K + i_h * K + o_k4, mask=(o_k4 < K), other=0.).to(tl.float32)
if STATE_V_FIRST:
b_h4 *= exp2(b_gk_last4)[None, :]
else:
b_h4 *= exp2(b_gk_last4)[:, None]
b_v = b_v.to(k.dtype.element_ty)

p_k = tl.make_block_ptr(k, (K, T), (1, H*K), (0, i_t * BT), (64, BT), (0, 1))
b_k = tl.load(p_k, boundary_check=(0, 1))
if STATE_V_FIRST:
b_h1 += tl.trans(tl.dot(b_k, b_v))
else:
b_h1 += tl.dot(b_k, b_v)
if K > 64:
p_k = tl.make_block_ptr(k, (K, T), (1, H*K), (64, i_t * BT), (64, BT), (0, 1))
b_k = tl.load(p_k, boundary_check=(0, 1))
if STATE_V_FIRST:
b_h2 += tl.trans(tl.dot(b_k, b_v))
else:
b_h2 += tl.dot(b_k, b_v)
if K > 128:
p_k = tl.make_block_ptr(k, (K, T), (1, H*K), (128, i_t * BT), (64, BT), (0, 1))
b_k = tl.load(p_k, boundary_check=(0, 1))
if STATE_V_FIRST:
b_h3 += tl.trans(tl.dot(b_k, b_v))
else:
b_h3 += tl.dot(b_k, b_v)
if K > 192:
p_k = tl.make_block_ptr(k, (K, T), (1, H*K), (192, i_t * BT), (64, BT), (0, 1))
b_k = tl.load(p_k, boundary_check=(0, 1))
if STATE_V_FIRST:
b_h4 += tl.trans(tl.dot(b_k, b_v))
else:
b_h4 += tl.dot(b_k, b_v)

if STORE_FINAL_STATE:
if STATE_V_FIRST:
p_ht = tl.make_block_ptr(ht, (V, K), (K, 1), (i_v * BV, 0), (BV, 64), (1, 0))
else:
p_ht = tl.make_block_ptr(ht, (K, V), (V, 1), (0, i_v * BV), (64, BV), (1, 0))
tl.store(p_ht, b_h1.to(p_ht.dtype.element_ty), boundary_check=(0, 1))
if K > 64:
if STATE_V_FIRST:
p_ht = tl.make_block_ptr(ht, (V, K), (K, 1), (i_v * BV, 64), (BV, 64), (1, 0))
else:
p_ht = tl.make_block_ptr(ht, (K, V), (V, 1), (64, i_v * BV), (64, BV), (1, 0))
tl.store(p_ht, b_h2.to(p_ht.dtype.element_ty), boundary_check=(0, 1))
if K > 128:
if STATE_V_FIRST:
p_ht = tl.make_block_ptr(ht, (V, K), (K, 1), (i_v * BV, 128), (BV, 64), (1, 0))
else:
p_ht = tl.make_block_ptr(ht, (K, V), (V, 1), (128, i_v * BV), (64, BV), (1, 0))
tl.store(p_ht, b_h3.to(p_ht.dtype.element_ty), boundary_check=(0, 1))
if K > 192:
if STATE_V_FIRST:
p_ht = tl.make_block_ptr(ht, (V, K), (K, 1), (i_v * BV, 192), (BV, 64), (1, 0))
else:
p_ht = tl.make_block_ptr(ht, (K, V), (V, 1), (192, i_v * BV), (64, BV), (1, 0))
tl.store(p_ht, b_h4.to(p_ht.dtype.element_ty), boundary_check=(0, 1))

这个 kernel进行chunk间隐状态递推

chunk间的隐状态必须串行递推: h_0 → h_1 → h_2 → ... → h_{NT-1}

公式

其中 v_new = u - w @ h_c

每个program处理:

  • 一个序列(i_n)的一个HV head(i_h)
  • V维度的一个分块 (i_v × BV(i_v+1) × BV)
  • K维度全部(通过b_h1/b_h2/b_h3/b_h4覆盖最多256维)
  • 时间维度上串行遍历所有chunk (for i_t in range(NT))

关键: 因为chunk间有依赖(h_c依赖h_{c-1}),时间维度不能并行,必须在一个program内串行遍历。

隐状态 h 的形状是 [K, V],每个program处理V的一个BV块,所以需要持有 [K, BV] 大小的隐状态。

但K可能很大(最多256),不能一次放进一个2D tile里做矩阵乘。所以把K维度拆成最多4个64维的块:

形状 K维度范围
b_h1 [64, BV] 0 ~ 63
b_h2 [64, BV] 64 ~ 127 (K>64时)
b_h3 [64, BV] 128 ~ 191 (K>128时)
b_h4 [64, BV] 192 ~ 255 (K>192时)

这4块拼起来就是完整的 h: [K, BV].

ASCII示意:

1
2
3
4
5
6
7
8
9
10
11
12
完整h: [K=256, BV=32]

拆分成:
┌─────────────────────────────┐
│ b_h1: [64, 32] K=[0:64) │
├─────────────────────────────┤
│ b_h2: [64, 32] K=[64:128) │
├─────────────────────────────┤
│ b_h3: [64, 32] K=[128:192)│
├─────────────────────────────┤
│ b_h4: [64, 32] K=[192:256)│
└─────────────────────────────┘
  • 如果 K=128, 只用 b_h1 和 b_h2
  • 如果 K=64, 只用 b_h1

block 视角

img

img

完整的一次循环迭代示意

假设 K=64(只有b_h1), BV=32, BT=64, USE_GK=True

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
chunk c 的处理流程:

1. store h[c] = b_h1 # 存当前隐状态 [64, 32]

2. b_w = load w[c*64:(c+1)*64, :] # [64, 64]
b_wh = dot(b_w, b_h1) # [64, 64] x [64, 32] = [64, 32]

3. b_u = load u[c*64:(c+1)*64, :] # [64, 32]
b_v_new = b_u - b_wh # [64, 32]
store v_new[c] = b_v_new

4. gk_last = load gk[c*64+63, :] # [64]
b_h1 *= exp2(gk_last)[:,None] # h衰减 [64,32]

5. b_k = load kg[:, c*64:(c+1)*64] # [64, 64], K x BT 布局
b_h1 += dot(b_k, b_v_new) # [64,64] x [64,32] = [64,32]

--> b_h1 现在是 h_{c+1}, 进入下一次循环

5 chunk_gla_fwd_o_gk

代码链接:flash-linear-attention/fla/ops/gla/chunk.py at main · fla-org/flash-linear-attention

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def chunk_gla_fwd_o_gk(
q: torch.Tensor,
v: torch.Tensor,
g: torch.Tensor,
A: torch.Tensor,
h: torch.Tensor,
scale: float,
state_v_first: bool = False,
cu_seqlens: torch.LongTensor | None = None,
chunk_size: int = 64,
chunk_indices: torch.LongTensor | None = None,
):
B, T, H, K, HV, V = *q.shape, v.shape[2], v.shape[-1]
BT = chunk_size

if chunk_indices is None and cu_seqlens is not None:
chunk_indices = prepare_chunk_indices(cu_seqlens, chunk_size)
NT = triton.cdiv(T, BT) if cu_seqlens is None else len(chunk_indices)

# Please ensure zeros, since vllm will use padding v
o = torch.zeros_like(v)
def grid(meta): return (triton.cdiv(V, meta['BV']), NT, B * HV)
chunk_gla_fwd_kernel_o[grid](
q=q,
v=v,
g=g,
h=h,
o=o,
A=A,
cu_seqlens=cu_seqlens,
chunk_indices=chunk_indices,
scale=scale,
T=T,
H=H,
HV=HV,
K=K,
V=V,
BT=BT,
STATE_V_FIRST=state_v_first,
)
return o

@triton.jit(do_not_specialize=['T'])
def chunk_gla_fwd_kernel_o(
q,
v,
g,
h,
o,
A,
cu_seqlens,
chunk_indices,
scale,
T,
H: tl.constexpr,
HV: tl.constexpr,
K: tl.constexpr,
V: tl.constexpr,
BT: tl.constexpr,
BK: tl.constexpr,
BV: tl.constexpr,
STATE_V_FIRST: tl.constexpr,
IS_VARLEN: tl.constexpr,
):
i_v, i_t, i_bh = tl.program_id(0), tl.program_id(1), tl.program_id(2)
i_b, i_hv = i_bh // HV, i_bh % HV
i_h = i_hv // (HV // H)
if IS_VARLEN:
i_tg = i_t.to(tl.int64)
i_n, i_t = tl.load(chunk_indices + i_t * 2).to(tl.int32), tl.load(chunk_indices + i_t * 2 + 1).to(tl.int32)
bos, eos = tl.load(cu_seqlens + i_n).to(tl.int64), tl.load(cu_seqlens + i_n + 1).to(tl.int64)
T = eos - bos
NT = tl.cdiv(T, BT)
else:
NT = tl.cdiv(T, BT)
i_tg = (i_b * NT + i_t).to(tl.int64)
bos, eos = (i_b * T).to(tl.int64), (i_b * T + T).to(tl.int64)

m_s = tl.arange(0, BT)[:, None] >= tl.arange(0, BT)[None, :]

q += (bos * H + i_h) * K
g += (bos * HV + i_hv) * K
v += (bos * HV + i_hv) * V
o += (bos * HV + i_hv) * V
h += (i_tg * HV + i_hv).to(tl.int64) * K * V
A += (bos * HV + i_hv) * BT

b_o = tl.zeros([BT, BV], dtype=tl.float32)
for i_k in range(tl.cdiv(K, BK)):
p_q = tl.make_block_ptr(q, (T, K), (H*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
p_g = tl.make_block_ptr(g, (T, K), (HV*K, 1), (i_t * BT, i_k * BK), (BT, BK), (1, 0))
if STATE_V_FIRST:
p_h = tl.make_block_ptr(h, (V, K), (K, 1), (i_v * BV, i_k * BK), (BV, BK), (1, 0))
else:
p_h = tl.make_block_ptr(h, (K, V), (V, 1), (i_k * BK, i_v * BV), (BK, BV), (1, 0))

# [BT, BK]
b_q = tl.load(p_q, boundary_check=(0, 1))
# [BT, BK]
b_g = tl.load(p_g, boundary_check=(0, 1)).to(tl.float32)
# [BT, BK]
b_qg = (b_q * exp2(b_g)).to(b_q.dtype)
b_h = tl.load(p_h, boundary_check=(0, 1))
if i_k >= 0:
if STATE_V_FIRST:
b_o += tl.dot(b_qg, tl.trans(b_h).to(b_qg.dtype))
else:
b_o += tl.dot(b_qg, b_h.to(b_qg.dtype))
b_o *= scale
p_v = tl.make_block_ptr(v, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
p_o = tl.make_block_ptr(o, (T, V), (HV*V, 1), (i_t * BT, i_v * BV), (BT, BV), (1, 0))
p_A = tl.make_block_ptr(A, (T, BT), (HV*BT, 1), (i_t * BT, 0), (BT, BT), (1, 0))
# [BT, BV]
b_v = tl.load(p_v, boundary_check=(0, 1))
# [BT, BT]
b_A = tl.load(p_A, boundary_check=(0, 1))
b_A = tl.where(m_s, b_A, 0.).to(b_v.dtype)
b_o += tl.dot(b_A, b_v)
tl.store(p_o, b_o.to(p_o.dtype.element_ty), boundary_check=(0, 1))

img

  • 标题: KDA源码剖析之FLA
  • 作者: 鱿鱼圈
  • 创建于 : 2026-06-04 23:50:00
  • 更新于 : 2026-06-14 23:33:01
  • 链接: https://yuyanqi.com/2026/06/04/KDA源码剖析之FLA/
  • 版权声明: 本文章采用 CC BY-NC-SA 4.0 进行许可。
评论