Generic Matmul

Generic route 是 Matmul 的語意基準。Planner 把每個合法 input 翻成同一份 位址規則；scalar loop、BLAS 與 packing route 因而共用相同的 rank、 broadcasting 與 layout contract。

Roles

兩個 operands 各自可能是 vector 或 matrix。公開介面因此有四種結果形狀， executor 卻只需要一種 contraction：M x K @ K x N。

• DOT：(K,) @ (K,) -> (1,)，plan 使用 M=1, N=1。
• GEVM：(K,) @ (...,K,N) -> (...,N)，plan 使用 M=1。
• GEMV：(...,M,K) @ (K,) -> (...,M)，plan 使用 N=1。
• GEMM：(...,M,K) @ (...,K,N) -> (...,M,N)。

M=1 或 N=1 是 loop bound，不是新配置出來的 axis。建立 public output shape 時，vector 缺少的 axis 會再被拿掉；DOT 則保留 SimpleArray 既有的 (1,) 結果。

以下是 planner 真正做 role normalization 的核心：

C++

bool const lhs_vector = lhs.ndim() == 1;
bool const rhs_vector = rhs.ndim() == 1;

ssize_t const rows = lhs_vector ? 1 : lhs.shape(lhs.ndim() - 2);
ssize_t const columns = rhs_vector ? 1 : rhs.shape(rhs.ndim() - 1);
ssize_t const inner_size = lhs.shape(lhs.ndim() - 1);

Broadcast

考慮一個兩邊都需要 broadcasting 的 GEMM：

lhs       (2, 1, 3, 4)
rhs       (1, 5, 4, 6)

batch     (2, 1) x (1, 5)  ->  (2, 5)
core          (3, 4) x (4, 6)  ->  (3, 6)
output                            (2, 5, 3, 6)

Batch axes 從右側對齊。planner 不會 expand operands，而是替共同的 (2,5) 座標空間建立三組 offset mappings：

buffer	mapped batch strides	offset at (b0,b1)
output	(90,18)	90*b0 + 18*b1
lhs	(12,0)	12*b0
rhs	(0,24)	24*b1

所以 (b0,b1)=(1,3) 對應 output、lhs、rhs offsets 144、12、72。 lhs 的第二個 mapped stride 是 0，表示五個 b1 coordinates 都重用同一個 matrix；rhs 的第一個 stride 也一樣。

這個 0 是 planner 產生的 broadcasting mapping，不代表 input buffer 的 physical stride 原本就是零。建構 mapping 的分支也刻意保留這個差別：

C++

if (operand_extent == domain.extent(domain_axis))
{
    strides[domain_axis] = operand_stride;
}
else if (domain.extent(domain_axis) > 1)
{
    is_broadcast = true; // mapped stride stays zero
}

Signed Strides

Batch mapping 先選出一對 matrices，core strides 再決定 contraction 中每一步 往哪裡讀。考慮從 shape (1,4,4) 的 base 建立 base[:, :, ::-1]，再把它 broadcast 給兩個 RHS matrices：

lhs shape (1,4,4), core strides (4,-1)
rhs shape (2,4,3), core strides (3,1)
batch mappings: lhs (0), rhs (12), output (12)

Batch stride `0` 重用同一個 LHS；core stride `-1` 則讓每一步沿原始 buffer 反向移動。兩種 stride 解決的是不同問題。

對 output C[1,2,1]，batch mapping 讓 LHS 留在 base 0，RHS 前進到 base 12，output block 前進到 base 12。接著 (i,j)=(2,1) 由四個 core strides 展開：

A(i,k) = A0 + lhs_base + i*lhs_row_stride + k*lhs_inner_stride
B(k,j) = B0 + rhs_base + k*rhs_inner_stride + j*rhs_column_stride

C[19] = A0[8] * B0[13]
      + A0[7] * B0[16]
      + A0[6] * B0[19]
      + A0[5] * B0[22]

Negative stride 本身不等於 cache-hostile。這個 LHS 仍然逐一讀取相鄰 elements， 只是方向相反。真正的 performance gap 是目前的 BLAS matrix-view predicate 無法表示 negative core stride，因此這個 layout 會落到沒有顯式 blocking、 SIMD kernel 或 vendor BLAS 的 scalar fallback。Packing 的目的不是修復 locality，而是把同一組 values 轉成 BLAS 能直接消費的 representation。

Array 另外說明 buffer + shape + strides + logical origin 如何表示 transpose、slice 和 reverse view。Matmul planner 只消費這份 descriptor， 不必知道 view 經過哪些 Python slicing operations。

Kernel

走到 hot loop 時，rank、broadcasting 和 layout 判斷都已經消失。每一個 output cell 只剩 signed offsets 與 multiply-add：

C++

value_type total{};
ssize_t lhs_offset = lhs_row_base;
ssize_t rhs_offset =
    rhs_base + column * m_plan.rhs_column_stride();

for (ssize_t inner = 0; inner < m_plan.inner_size(); ++inner)
{
    total += m_lhs_data[lhs_offset] * m_rhs_data[rhs_offset];
    lhs_offset += m_plan.lhs_inner_stride();
    rhs_offset += m_plan.rhs_inner_stride();
}

K=0 時，value-initialized accumulator 會把 zero 寫回每個非空 output cell； M=0、N=0 或 empty batch 則不進入對應的外層 loop。這條 scalar route 不做 blocking 或 packing，它的價值是讓所有 optimized routes 共用一個可檢查的 semantic baseline。

Tests

#1208 用 NumPy differential tests 交叉組合四種 roles、right-aligned broadcasting、empty dimensions，以及 C-contiguous、Fortran、negative-stride 和 step-two views。之後每增加一條 optimized route，都必須對同一份 public contract 產生相同 shape 與 values。

Source

• MatmulPlan
• LoopDomain and MappedOffsetCursor
• execute_generic()