Dispatch

Dispatch 不改變矩陣乘法的語意。MatmulPlan 已經驗證 shape、解析 broadcasting，並把 operands 表示成 M x K @ K x N。Dispatch 只回答下一個 問題：這份 plan 應該由哪一個 kernel 執行，開始前是否需要 packing？

答案是一個很小的 value：

C++

struct MatmulSelection
{
    MatmulKernel kernel = MatmulKernel::Generic;
    PackingState packing;
};

它在 batch traversal 前只計算一次。若要 packing，executor 先建立新的 row-major operands 並重建 plan；之後的每個 batch coordinate 都使用已固定的 kernel，不會在 inner loop 裡重新判斷 layout。

Facts

Selector 讀取的是 plan 中已經確定的 facts，不是 array 曾經經過哪些 slice 或 transpose operations：

fact	它回答的問題
role	這是 DOT、GEVM、GEMV 還是 GEMM？
M, N, K	一次 contraction 有多大？
batch domain	kernel 會被呼叫幾次？
core strides	operand 能否直接表示成 BLAS matrix 或 vector view？
shape broadcasting	同一份 logical data 是否會被多個 batch coordinates 重用？
dtype 與 backend	目前 build 是否有相應的 BLAS function？

這些 facts 依下列順序收斂成 route：

1. 沒有可用 BLAS，或 dtype 不受 BLAS route 支援，選 Generic。
2. 依 role 檢查 matrix view 與 vector stride 是否可由 BLAS 表示。
3. 若 layout 可直接使用，而且 workload 已跨過目前實測的 crossover boundary，選 direct BLAS。
4. 若 layout 不能直接使用，檢查不相容的 operand 是否能安全地 pack，且 packed storage 是否會被重用。成立時選 pack-once BLAS。
5. 其餘情況回到 Generic。它是完整的 signed-stride implementation，不是錯誤 路徑。

不同 role 可用的 routes 並不完全相同：

• DOT：兩個 vector strides 符合 policy 時 direct BLAS；目前不 packing。
• GEVM / GEMV：matrix 和 vector stride 都可表示時 direct BLAS；batched operation 若重用 non-positive-stride vector，才考慮 pack once。
• GEMM：兩個 matrices 都可表示時 direct BLAS；否則只考慮 packing 被 shape broadcasting 重用的不相容 matrix。

不符合上述條件的任何 role 都回到 Generic。

Layout

BLAS matrix view 需要 pointer、leading dimension 與 transpose flag。solvcon 接受 row-major-compatible 與 column-major-compatible 兩種 core layout；padding 可以存在，但其中一個 axis 必須是 unit stride。merged implementation 的判斷就是：

C++

if (column_stride == 1 && row_stride >= columns)
{
    return BlasMatrixView<value_type>{
        data, row_stride, BlasTranspose::None};
}
if (row_stride == 1 && column_stride >= rows)
{
    return BlasMatrixView<value_type>{
        data, column_stride, BlasTranspose::Transpose};
}
return std::nullopt;

第二種沒有搬動或 transpose data。它只是交換 BLAS 看待 axes 的方式。negative stride、step-two inner stride，或兩個 core axes 都不是 unit stride 時，這個 predicate 會回傳 nullopt。

Vector 的規則由各 role 決定。DOT 目前只直接接受兩邊都是 +1，或兩邊都是 -1；GEVM / GEMV 的 BLAS wrapper 可以使用 positive vector stride。batched GEVM / GEMV 遇到被重用的 non-positive-stride vector 時，才考慮先 pack vector。

Reuse

「某個 batch offset 是零」不一定代表可以 pack once。Planner 區分兩件事：

• shape broadcasting 產生的 mapping stride 0：例如 (1,M,K) 對 (B,K,N)，extent 1 讓同一個 matrix 合法地重用 B 次；
• 原始 view 在 extent 大於 1 的 batch axis 上有 physical stride 0：多個 logical indices alias 同一份 storage。

目前 matrix-packing route 只接受第一種。若不相容的 operand 並未被 shape broadcast，或原始 batch view 自己有 zero stride，就保留 Generic。exact predicate 如下：

C++

bool const lhs_supported =
    !required.lhs ||
    (m_plan.lhs_is_broadcast() && !m_plan.lhs_has_zero_batch_stride());
bool const rhs_supported =
    !required.rhs ||
    (m_plan.rhs_is_broadcast() && !m_plan.rhs_has_zero_batch_stride());

if (!required || !lhs_supported || !rhs_supported)
{
    return PackingState{};
}
return required;

required.lhs 或 required.rhs 只標記無法建立 BLAS matrix view 的 operand。 Packing 不會複製已經相容的一側。完成後必須重建 plan，因為 logical values 沒有 變，但 core strides 與 data pointer 已經改變：

C++

if (packing.lhs)
{
    m_packed_lhs.emplace(m_lhs.to_row_major());
}
if (packing.rhs)
{
    m_packed_rhs.emplace(m_rhs.to_row_major());
}

Array const & lhs = m_packed_lhs ? *m_packed_lhs : m_lhs;
Array const & rhs = m_packed_rhs ? *m_packed_rhs : m_rhs;
m_plan = MatmulPlan::make(lhs, rhs);

Packing 本身不是新的矩陣乘法技巧。這裡的設計重點是 plan 能在 batch traversal 以前，從 broadcasting mapping 算出哪些 supplied operands 會被重用；packing 因此按 supplied operand 執行，而不是按 output coordinate 重做。

Cost

Selector 無法從 shape 精確預測時間，但 plan 能先列出一張不含未知常數的工作帳。 對一個 batched GEMM，定義：

• C：output batch domain 中的 contraction 數量；
• U_lhs：supplied LHS 中會被走訪的 matrices 數量；
• U_rhs：supplied RHS 中會被走訪的 matrices 數量；
• R_lhs = C/U_lhs、R_rhs = C/U_rhs：每份 supplied matrix 的平均 reuse。

若 element size 是 s bytes，三條 route 的可數成本是：

• Generic：執行 C 次 scalar contractions，不配置 packed payload。
• Direct BLAS：呼叫 BLAS C 次，不配置 packed payload。
• Pack-once BLAS：先整理不相容的 supplied matrices，再呼叫 BLAS C 次。 LHS payload 是 s*U_lhs*M*K，RHS payload 是 s*U_rhs*K*N。

Payload 是 packed storage 的大小。Copy 至少還要讀一次 source、寫一次 destination， 所以實際 memory traffic 至少是 payload 的兩倍，另有 allocation、cache 與 BLAS call overhead。對 real-valued operands，C 個 GEMM 約做 2*C*M*N*K FLOPs。 這兩組可數資料不會直接給出秒數，但能先排除顯然不合理的 route，再用 profiling 決定 crossover。

Cross broadcasting 說明為什麼不能只用一個模糊的 batch size。對：

Text Only

lhs batch (2,1) x rhs batch (1,5) -> output batch (2,5)

工作帳是 C=10、U_lhs=2、U_rhs=5，所以每個 LHS reuse 5 次，每個 RHS reuse 2 次。若兩邊都不相容，pack-once route 建立 2+5=7 份 packed matrices， payload 是 s*(2*M*K + 5*K*N)，然後做十次 BLAS。它不是把 output batch size 10 誤當成兩邊各自的 reuse，也不會為十個 coordinates 各 pack 一對 operands。

本文的 S=32 LHS-broadcast case 則有 C=10、U_lhs=1、U_rhs=10，只有 LHS 不相容。Float64 LHS 的 packed payload 是 1*32*32*8 = 8192 bytes；copy 至少產生約 16 KiB read-plus-write traffic，後面接十個 GEMM，約 655,360 FLOPs。也就是每個 payload byte 後面約有 80 FLOPs（以最低 copy traffic 算則 約 40 FLOPs/byte）。這張 ledger 解釋 reuse 如何攤提 representation cost； 它仍不能取代實測，因為 scalar traversal、BLAS call、cache 與 allocator 的時間 不由 shape 單獨決定。

操作上可以先做三個判斷：layout 已相容時，packing 沒有價值；layout 不相容但 R=1 時，copy 沒有跨 contractions 攤提；layout 不相容且 R>1 時，再比較 packed bytes、後續 arithmetic 與實測 boundary。Generic 始終保留，因為小型 workload 或無法安全重用的 view 可能不值得 materialize。

Trace

以下 workload 讓一個 reversed LHS 在十個 GEMM 中重用：

Text Only

lhs  shape (1,32,32), element strides (1024,32,-1)
rhs  shape (10,32,32), C-contiguous
out  shape (10,32,32)
dtype float64, BLAS available

沿著 merged selector 逐步判斷：

1. Role 是 GEMM，M=N=K=32；工作帳是 C=10、U_lhs=1、 R_lhs=10。
2. RHS core strides (32,1) 可以 direct BLAS；LHS (32,-1) 無法通過 matrix-view predicate，所以 direct GEMM 不成立。
3. min(M,N,K)=32 跨過目前 packing boundary 16。LHS 又是 shape broadcasting reuse，且沒有 physical zero batch stride，因此 selection 是 {BlasGemm, pack lhs}。
4. to_row_major() 建立一份 LHS，重建 plan 後 core strides 成為 (32,1)； 十個 batch coordinates 隨後共用這份 storage 與固定的 BlasGemm kernel。

這個 trace 也指出 Generic 仍會被選中的鄰近 cases：min(M,N,K)<16、reversed matrix 沒有 shape-broadcast reuse、原始 batch view 使用 physical zero stride， 或 build / dtype 沒有 BLAS support。哪條 route 較快是 policy 問題；哪條 route 算得對仍由共同的 plan 保證。

Policy

下面的 table 只能用來預測 snapshot 95a68122 目前會選哪條 code route，不能 預測另一台機器的 portable optimum。這些 compile-time boundaries 由 PR #1209、 PR #1256 與 PR #1259 的 Accelerate profiling 形成。目前程式尚未依 backend 或 dtype 選擇不同 tuning table，所以這些數字描述 的是當前 implementation policy，不是 portable cost model，也不是 matmul semantics。

查看完整 boundary table

| route | boundary | | --- | --- | | DOT direct BLAS | `K >= 128`，兩個 vector strides 同為 `+1` 或同為 `-1` | | compact unbatched GEVM | BLAS-compatible matrix、positive vector stride；`K*N >= 729` | | other unbatched GEVM | 相同 layout 條件；`min(K,N) >= 32` | | unbatched GEMV | BLAS-compatible matrix、positive vector stride；`min(M,K) >= 32` | | GEMM direct BLAS | 兩個 matrix views 都相容；`min(M,N,K) >= 8` | | GEMM pack once | 每個不相容 matrix 都符合 shape-broadcast reuse 且無 physical zero batch stride；`min(M,N,K) >= 16` | | batched GEVM / GEMV direct BLAS | matrix view 相容、positive vector stride；matrix elements `>= 512` | | reused non-positive-stride vector: always pack | matrix view 相容；matrix elements `>= 4096` | | reused non-positive-stride vector: conditional | matrix elements `>= 1024` 且 batch size `>= 4` | | reused non-positive-stride vector: reuse floor | matrix elements `>= 576` 且 total output elements `>= 128` |

量測流程、raw data 與 crossover evidence 集中在 Benchmarks。這一頁保留 boundaries，是為了讓讀者能把 上面的 trace 與 merged code 對上，而不是把 constants 當成跨機器的建議值。

Source

本文對應 solvcon upstream commit 95a68122：

• execute() and role selection
• select_matrix_packing() and pack()
• execute_contractions() and BLAS calls
• make_matrix_view()