# Adapted from: https://github.com/vllm-project/vllm/tree/main/vllm/model_executor/layers/mamba/ops/ssd_chunk_state.py # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright contributors to the vLLM project # Copyright (c) 2024, Tri Dao, Albert Gu. # Adapted from https://github.com/state-spaces/mamba/blob/v2.2.4/mamba_ssm/ops/triton/ssd_chunk_state.py # ruff: noqa: E501 import math import torch import triton import triton.language as tl from .mamba_ssm import softplus @triton.jit def _chunk_cumsum_fwd_kernel( # Pointers to matrices dt_ptr, A_ptr, dt_bias_ptr, dt_out_ptr, dA_cumsum_ptr, # Matrix dimension batch, seqlen, nheads, chunk_size, dt_min, dt_max, # Strides stride_dt_batch, stride_dt_seqlen, stride_dt_head, stride_A_head, stride_dt_bias_head, stride_dt_out_batch, stride_dt_out_chunk, stride_dt_out_head, stride_dt_out_csize, stride_dA_cs_batch, stride_dA_cs_chunk, stride_dA_cs_head, stride_dA_cs_csize, # Meta-parameters DT_SOFTPLUS: tl.constexpr, HAS_DT_BIAS: tl.constexpr, BLOCK_SIZE_CHUNK: tl.constexpr, BLOCK_SIZE_H: tl.constexpr = 16, ): pid_b = tl.program_id(axis=0) # if dt is long, may cause problems, so use 64 bit # https://github.com/triton-lang/triton/issues/1058 pid_c = tl.program_id(axis=1).to(tl.int64) pid_h = tl.program_id(axis=2) dt_ptr += pid_b * stride_dt_batch + pid_c * chunk_size * stride_dt_seqlen dt_out_ptr += pid_b * stride_dt_out_batch + pid_c * stride_dt_out_chunk dA_cumsum_ptr += pid_b * stride_dA_cs_batch + pid_c * stride_dA_cs_chunk offs_h = pid_h * BLOCK_SIZE_H + tl.arange(0, BLOCK_SIZE_H) offs_c = tl.arange(0, BLOCK_SIZE_CHUNK) dt_ptrs = dt_ptr + ( offs_h[:, None] * stride_dt_head + offs_c[None, :] * stride_dt_seqlen ) A_ptrs = A_ptr + offs_h * stride_A_head dt_out_ptrs = dt_out_ptr + ( offs_h[:, None] * stride_dt_out_head + offs_c[None, :] * stride_dt_out_csize ) dA_cs_ptrs = dA_cumsum_ptr + ( offs_h[:, None] * stride_dA_cs_head + offs_c[None, :] * stride_dA_cs_csize ) chunk_size_limit = min(chunk_size, seqlen - pid_c * chunk_size) dt = tl.load( dt_ptrs, mask=(offs_h[:, None] < nheads) & (offs_c[None, :] < chunk_size_limit), other=0.0, ).to(tl.float32) if HAS_DT_BIAS: dt_bias = tl.load( dt_bias_ptr + offs_h * stride_dt_bias_head, mask=offs_h < nheads, other=0.0 ).to(tl.float32) dt += dt_bias[:, None] if DT_SOFTPLUS: dt = tl.where(dt <= 20.0, softplus(dt), dt) # As of Triton 2.2.0, tl.clamp is not available yet # dt = tl.clamp(dt, dt_min, dt_max) dt = tl.minimum(tl.maximum(dt, dt_min), dt_max) dt = tl.where( (offs_h[:, None] < nheads) & (offs_c[None, :] < chunk_size_limit), dt, 0.0 ) tl.store( dt_out_ptrs, dt, mask=(offs_h[:, None] < nheads) & (offs_c[None, :] < chunk_size), ) A = tl.load(A_ptrs, mask=offs_h < nheads, other=0.0).to(tl.float32) dA = dt * A[:, None] dA_cs = tl.cumsum(dA, axis=1) tl.store( dA_cs_ptrs, dA_cs, mask=(offs_h[:, None] < nheads) & (offs_c[None, :] < chunk_size), ) @triton.jit def _chunk_state_fwd_kernel( # Pointers to matrices x_ptr, b_ptr, states_ptr, dt_ptr, dA_cumsum_ptr, seq_idx_ptr, # Matrix dimensions hdim, dstate, chunk_size, batch, seqlen, nheads_ngroups_ratio, # Strides stride_x_batch, stride_x_seqlen, stride_x_head, stride_x_hdim, stride_b_batch, stride_b_seqlen, stride_b_head, stride_b_dstate, stride_states_batch, stride_states_chunk, stride_states_head, stride_states_hdim, stride_states_dstate, stride_dt_batch, stride_dt_chunk, stride_dt_head, stride_dt_csize, stride_dA_cs_batch, stride_dA_cs_chunk, stride_dA_cs_head, stride_dA_cs_csize, stride_seq_idx_batch, stride_seq_idx_seqlen, # Meta-parameters HAS_SEQ_IDX: tl.constexpr, BLOCK_SIZE_M: tl.constexpr = 16, BLOCK_SIZE_N: tl.constexpr = 16, BLOCK_SIZE_K: tl.constexpr = 16, ): pid_bc = tl.program_id(axis=1).to(tl.int64) pid_c = pid_bc // batch pid_b = pid_bc - pid_c * batch pid_h = tl.program_id(axis=2) num_pid_n = tl.cdiv(dstate, BLOCK_SIZE_N) pid_m = tl.program_id(axis=0) // num_pid_n pid_n = tl.program_id(axis=0) % num_pid_n b_ptr += ( pid_b * stride_b_batch + pid_c * chunk_size * stride_b_seqlen + (pid_h // nheads_ngroups_ratio) * stride_b_head ) x_ptr += ( pid_b * stride_x_batch + pid_c * chunk_size * stride_x_seqlen + pid_h * stride_x_head ) dt_ptr += pid_b * stride_dt_batch + pid_c * stride_dt_chunk + pid_h * stride_dt_head dA_cumsum_ptr += ( pid_b * stride_dA_cs_batch + pid_c * stride_dA_cs_chunk + pid_h * stride_dA_cs_head ) if HAS_SEQ_IDX: seq_idx_ptr += ( pid_b * stride_seq_idx_batch + pid_c * chunk_size * stride_seq_idx_seqlen ) offs_m = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M) offs_n = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N) offs_k = tl.arange(0, BLOCK_SIZE_K) x_ptrs = x_ptr + ( offs_m[:, None] * stride_x_hdim + offs_k[None, :] * stride_x_seqlen ) b_ptrs = b_ptr + ( offs_n[None, :] * stride_b_dstate + offs_k[:, None] * stride_b_seqlen ) dt_ptrs = dt_ptr + offs_k * stride_dt_csize dA_cs_last = tl.load(dA_cumsum_ptr + (chunk_size - 1) * stride_dA_cs_csize).to( tl.float32 ) dA_cumsum_ptrs = dA_cumsum_ptr + offs_k * stride_dA_cs_csize if HAS_SEQ_IDX: seq_idx_ptrs = seq_idx_ptr + offs_k * stride_seq_idx_seqlen chunk_size_limit = min(chunk_size, seqlen - pid_c * chunk_size) if HAS_SEQ_IDX: seq_idx_last = tl.load( seq_idx_ptr + (chunk_size_limit - 1) * stride_seq_idx_seqlen ) acc = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=tl.float32) for k in range(0, chunk_size_limit, BLOCK_SIZE_K): x = tl.load( x_ptrs, mask=(offs_m[:, None] < hdim) & (offs_k[None, :] < chunk_size_limit - k), other=0.0, ) b = tl.load( b_ptrs, mask=(offs_k[:, None] < chunk_size_limit - k) & (offs_n[None, :] < dstate), other=0.0, ).to(tl.float32) dA_cs_k = tl.load( dA_cumsum_ptrs, mask=offs_k < chunk_size_limit - k, other=0.0 ).to(tl.float32) if HAS_SEQ_IDX: seq_idx_k = tl.load( seq_idx_ptrs, mask=offs_k < chunk_size_limit - k, other=-1 ) dt_k = tl.load(dt_ptrs, mask=offs_k < chunk_size_limit - k, other=0.0).to( tl.float32 ) if not HAS_SEQ_IDX: scale = tl.exp(dA_cs_last - dA_cs_k) * dt_k else: scale = tl.where( seq_idx_k == seq_idx_last, tl.exp(dA_cs_last - dA_cs_k) * dt_k, 0.0 ) b *= scale[:, None] b = b.to(x_ptr.dtype.element_ty) acc += tl.dot(x, b) x_ptrs += BLOCK_SIZE_K * stride_x_seqlen b_ptrs += BLOCK_SIZE_K * stride_b_seqlen dt_ptrs += BLOCK_SIZE_K * stride_dt_csize dA_cumsum_ptrs += BLOCK_SIZE_K * stride_dA_cs_csize if HAS_SEQ_IDX: seq_idx_ptrs += BLOCK_SIZE_K * stride_seq_idx_seqlen states = acc.to(states_ptr.dtype.element_ty) states_ptr += ( pid_b * stride_states_batch + pid_c * stride_states_chunk + pid_h * stride_states_head ) offs_m = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M) offs_n = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N) states_ptrs = states_ptr + ( offs_m[:, None] * stride_states_hdim + offs_n[None, :] * stride_states_dstate ) c_mask = (offs_m[:, None] < hdim) & (offs_n[None, :] < dstate) tl.store(states_ptrs, states, mask=c_mask) @triton.jit def _chunk_state_varlen_kernel( # Pointers to matrices x_ptr, b_ptr, dt_ptr, dA_cumsum_ptr, chunk_states_ptr, cu_seqlens_ptr, states_ptr, initstates_ptr, # Matrix dimensions hdim, dstate, chunk_size, seqlen, nheads_ngroups_ratio, # Strides stride_x_seqlen, stride_x_head, stride_x_hdim, stride_b_seqlen, stride_b_head, stride_b_dstate, stride_dt_chunk, stride_dt_head, stride_dt_csize, stride_dA_cs_chunk, stride_dA_cs_head, stride_dA_cs_csize, stride_chunk_states_chunk, stride_chunk_states_head, stride_chunk_states_hdim, stride_chunk_states_dstate, stride_states_batch, stride_states_head, stride_states_hdim, stride_states_dstate, stride_init_states_batch, stride_init_states_head, stride_init_states_hdim, stride_init_states_dstate, # Meta-parameters HAS_INITSTATES: tl.constexpr, BLOCK_SIZE_M: tl.constexpr = 16, BLOCK_SIZE_N: tl.constexpr = 16, BLOCK_SIZE_K: tl.constexpr = 16, ): pid_b = tl.program_id(axis=1) pid_h = tl.program_id(axis=2) num_pid_n = tl.cdiv(dstate, BLOCK_SIZE_N) pid_m = tl.program_id(axis=0) // num_pid_n pid_n = tl.program_id(axis=0) % num_pid_n end_idx = tl.load(cu_seqlens_ptr + pid_b + 1) pid_c = (end_idx - 1) // chunk_size b_ptr += ( pid_c * chunk_size * stride_b_seqlen + (pid_h // nheads_ngroups_ratio) * stride_b_head ) x_ptr += pid_c * chunk_size * stride_x_seqlen + pid_h * stride_x_head dt_ptr += pid_c * stride_dt_chunk + pid_h * stride_dt_head dA_cumsum_ptr += pid_c * stride_dA_cs_chunk + pid_h * stride_dA_cs_head chunk_states_ptr += ( pid_c * stride_chunk_states_chunk + pid_h * stride_chunk_states_head ) if HAS_INITSTATES: # if there are init states provided, we differentiate between states (which # are boundary conditions at a chunk boundary) and initstates (which are boundary # conditions when a new example in a cont batch starts) initstates_ptr += pid_h * stride_init_states_head offs_m = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M) offs_n = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N) offs_k = tl.arange(0, BLOCK_SIZE_K) x_ptrs = x_ptr + ( offs_m[:, None] * stride_x_hdim + offs_k[None, :] * stride_x_seqlen ) b_ptrs = b_ptr + ( offs_n[None, :] * stride_b_dstate + offs_k[:, None] * stride_b_seqlen ) dt_ptrs = dt_ptr + offs_k * stride_dt_csize dA_cs_last = tl.load( dA_cumsum_ptr + (end_idx - pid_c * chunk_size - 1) * stride_dA_cs_csize ).to(tl.float32) dA_cumsum_ptrs = dA_cumsum_ptr + offs_k * stride_dA_cs_csize chunk_size_limit = end_idx - pid_c * chunk_size start_idx = tl.load(cu_seqlens_ptr + pid_b) start_idx_cur = tl.maximum(start_idx - pid_c * chunk_size, 0) acc = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=tl.float32) for k in range(0, chunk_size_limit, BLOCK_SIZE_K): x = tl.load( x_ptrs, mask=(offs_m[:, None] < hdim) & (offs_k[None, :] < chunk_size_limit - k) & (offs_k[None, :] >= start_idx_cur - k), other=0.0, ) b = tl.load( b_ptrs, mask=(offs_k[:, None] < chunk_size_limit - k) & (offs_n[None, :] < dstate) & (offs_k[:, None] >= start_idx_cur - k), other=0.0, ).to(tl.float32) dA_cs_k = tl.load( dA_cumsum_ptrs, mask=offs_k < chunk_size_limit - k, other=0.0 ).to(tl.float32) dt_k = tl.load(dt_ptrs, mask=offs_k < chunk_size_limit - k, other=0.0).to( tl.float32 ) scale = tl.where( (offs_k >= start_idx_cur - k) & (offs_k < chunk_size_limit - k), tl.exp(dA_cs_last - dA_cs_k) * dt_k, 0.0, ) b *= scale[:, None] b = b.to(x_ptr.dtype.element_ty) acc += tl.dot(x, b) x_ptrs += BLOCK_SIZE_K * stride_x_seqlen b_ptrs += BLOCK_SIZE_K * stride_b_seqlen dt_ptrs += BLOCK_SIZE_K * stride_dt_csize dA_cumsum_ptrs += BLOCK_SIZE_K * stride_dA_cs_csize # If the sequence starts after the last chunk idx, we don't need to add the contribution from the last chunk # If HAS_INITSTATES==True need to consider two possibilities # - if start_idx < pid_c * chunk_size, then we need to take the past_states_ptrs # - if state_idx >= pid * chunk_size, then we need to insert initstates if (start_idx < pid_c * chunk_size) or (HAS_INITSTATES): # first chunk dA_cs_boundary = 0.0 # default if not HAS_INITSTATES: past_states_ptrs = chunk_states_ptr + ( offs_m[:, None] * stride_chunk_states_hdim + offs_n[None, :] * stride_chunk_states_dstate ) else: # - this seems repetitive, buts its to help the compiler if start_idx < pid_c * chunk_size: past_states_ptrs = chunk_states_ptr + ( offs_m[:, None] * stride_chunk_states_hdim + offs_n[None, :] * stride_chunk_states_dstate ) else: past_states_ptrs = initstates_ptr + ( pid_b * stride_init_states_batch + offs_m[:, None] * stride_init_states_hdim + offs_n[None, :] * stride_init_states_dstate ) # need to adjust the boundary if start_idx > pid_c * chunk_size: dA_cs_boundary = tl.load( dA_cumsum_ptr + (start_idx - pid_c * chunk_size - 1) * stride_dA_cs_csize ).to(tl.float32) past_states = tl.load( past_states_ptrs, mask=(offs_m[:, None] < hdim) & (offs_n[None, :] < dstate), other=0.0, ).to(tl.float32) scale = tl.exp(dA_cs_last - dA_cs_boundary) acc += past_states * scale states = acc.to(states_ptr.dtype.element_ty) states_ptr += pid_b * stride_states_batch + pid_h * stride_states_head offs_m = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M) offs_n = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N) states_ptrs = states_ptr + ( offs_m[:, None] * stride_states_hdim + offs_n[None, :] * stride_states_dstate ) c_mask = (offs_m[:, None] < hdim) & (offs_n[None, :] < dstate) tl.store(states_ptrs, states, mask=c_mask) def _chunk_cumsum_fwd( dt, A, chunk_size, dt_bias=None, dt_softplus=False, dt_limit=(0.0, float("inf")) ): batch, seqlen, nheads = dt.shape assert A.shape == (nheads,) if dt_bias is not None: assert dt_bias.shape == (nheads,) nchunks = math.ceil(seqlen / chunk_size) dt_out = torch.empty( batch, nheads, nchunks, chunk_size, device=dt.device, dtype=torch.float32 ) dA_cumsum = torch.empty( batch, nheads, nchunks, chunk_size, device=dt.device, dtype=torch.float32 ) grid_chunk_cs = lambda META: ( batch, nchunks, triton.cdiv(nheads, META["BLOCK_SIZE_H"]), ) with torch.cuda.device(dt.device.index): _chunk_cumsum_fwd_kernel[grid_chunk_cs]( dt, A, dt_bias, dt_out, dA_cumsum, batch, seqlen, nheads, chunk_size, dt_limit[0], dt_limit[1], dt.stride(0), dt.stride(1), dt.stride(2), A.stride(0), dt_bias.stride(0) if dt_bias is not None else 0, dt_out.stride(0), dt_out.stride(2), dt_out.stride(1), dt_out.stride(3), dA_cumsum.stride(0), dA_cumsum.stride(2), dA_cumsum.stride(1), dA_cumsum.stride(3), dt_softplus, HAS_DT_BIAS=dt_bias is not None, BLOCK_SIZE_CHUNK=triton.next_power_of_2(chunk_size), ) return dA_cumsum, dt_out def _chunk_state_fwd( B, x, dt, dA_cumsum, seq_idx=None, states=None, states_in_fp32=True ): batch, seqlen, nheads, headdim = x.shape _, _, nchunks, chunk_size = dt.shape _, _, ngroups, dstate = B.shape assert nheads % ngroups == 0 assert B.shape == (batch, seqlen, ngroups, dstate) assert dt.shape == (batch, nheads, nchunks, chunk_size) assert dA_cumsum.shape == dt.shape if seq_idx is not None: assert seq_idx.shape == (batch, seqlen) if states is not None: assert states.shape == (batch, nchunks, nheads, headdim, dstate) else: states_dtype = torch.float32 if states_in_fp32 else B.dtype states = torch.empty( (batch, nchunks, nheads, headdim, dstate), device=x.device, dtype=states_dtype, ) grid = lambda META: ( triton.cdiv(headdim, META["BLOCK_SIZE_M"]) * triton.cdiv(dstate, META["BLOCK_SIZE_N"]), batch * nchunks, nheads, ) with torch.cuda.device(x.device.index): _chunk_state_fwd_kernel[grid]( x, B, states, dt, dA_cumsum, seq_idx, headdim, dstate, chunk_size, batch, seqlen, nheads // ngroups, x.stride(0), x.stride(1), x.stride(2), x.stride(3), B.stride(0), B.stride(1), B.stride(2), B.stride(-1), states.stride(0), states.stride(1), states.stride(2), states.stride(3), states.stride(4), dt.stride(0), dt.stride(2), dt.stride(1), dt.stride(3), dA_cumsum.stride(0), dA_cumsum.stride(2), dA_cumsum.stride(1), dA_cumsum.stride(3), *( (seq_idx.stride(0), seq_idx.stride(1)) if seq_idx is not None else (0, 0) ), HAS_SEQ_IDX=seq_idx is not None, ) return states def chunk_state_varlen( B, x, dt, dA_cumsum, cu_seqlens, chunk_states, initial_states=None ): total_seqlen, nheads, headdim = x.shape _, nchunks, chunk_size = dt.shape _, ngroups, dstate = B.shape batch = cu_seqlens.shape[0] - 1 cu_seqlens = cu_seqlens.contiguous() assert nheads % ngroups == 0 assert B.shape == (total_seqlen, ngroups, dstate) assert dt.shape == (nheads, nchunks, chunk_size) assert dA_cumsum.shape == dt.shape assert chunk_states.shape == (nchunks, nheads, headdim, dstate) if initial_states is not None: assert initial_states.shape == (batch, nheads, headdim, dstate) states = torch.empty( batch, nheads, headdim, dstate, dtype=chunk_states.dtype, device=chunk_states.device, ) grid = lambda META: ( triton.cdiv(headdim, META["BLOCK_SIZE_M"]) * triton.cdiv(dstate, META["BLOCK_SIZE_N"]), batch, nheads, ) with torch.cuda.device(x.device.index): _chunk_state_varlen_kernel[grid]( x, B, dt, dA_cumsum, chunk_states, cu_seqlens, states, initial_states, headdim, dstate, chunk_size, total_seqlen, nheads // ngroups, x.stride(0), x.stride(1), x.stride(2), B.stride(0), B.stride(1), B.stride(2), dt.stride(1), dt.stride(0), dt.stride(2), dA_cumsum.stride(1), dA_cumsum.stride(0), dA_cumsum.stride(2), chunk_states.stride(0), chunk_states.stride(1), chunk_states.stride(2), chunk_states.stride(3), states.stride(0), states.stride(1), states.stride(2), states.stride(3), *( ( initial_states.stride(0), initial_states.stride(1), initial_states.stride(2), initial_states.stride(3), ) if initial_states is not None else (0, 0, 0, 0) ), HAS_INITSTATES=initial_states is not None, ) return states