在 solvcon 開發比 NumPy 更快的矩陣乘法

當一個不能直接表示成 BLAS view 的矩陣，會被 batch broadcasting 重複使用時， solvcon 的 Matmul 能比 NumPy 快。

最初的線索很反常：只加一個 ::-1，NumPy 的時間就從約 12 ms 變成 1998 ms。以下兩個 operands 會做 64 次 256 x 256 GEMM；LHS 的 values 沒有複製，只是沿最後一軸反向讀取：

Python

import numpy as np

rng = np.random.default_rng(0)
storage = rng.standard_normal((1, 256, 256))
lhs = storage[:, :, ::-1]
rhs = rng.standard_normal((64, 256, 256))
out = lhs @ rhs

lhs 的 shape 是 (1,256,256)，element strides 是 (65536,256,-1)； batch extent 1 broadcast 成 64，所以同一個 reversed matrix 被重用 64 次。

SimpleArray 會原樣保留這類 NumPy view 的 buffer、shape、 strides 與 logical origin。Slice 或 reverse 因而能直接流進後續 numerical kernels，不會自動 materialize 成 contiguous matrix。

Anomaly

Issue #1172 在 Apple M1、 Accelerate 單執行緒上，對同一個 (1,256,256) @ (64,256,256) 做了三組對照：

LHS layout	NumPy	solvcon prototype
dense	`11.960 ms`	`11.873 ms`
reversed columns	`1997.845 ms`	`12.099 ms`
reversed values，計時前排成 contiguous	`11.967 ms`	`11.950 ms`

Dense 時兩邊都約 12 ms。只把 columns reverse，NumPy 變成 1998 ms；保留相同 values 與 64 次 GEMM，但在計時前排成 contiguous，NumPy 又回到 12 ms。差距不在乘法規模、broadcast 次數或數值，只在資料排列。

Prototype 也只整理唯一一份 LHS，再讓 64 次 GEMM 共用，沒有改寫乘法。問題因此變成：同一筆資料整理成本要付很多次，還是只付一次？

Generic

Prototype 只證明方向。完整 Matmul 還要處理 DOT、GEVM、GEMV、GEMM、 right-aligned broadcasting，以及正、負或 step-two strides。solvcon 因此先做 generic route：plan 將四種 roles 正規化成 M x K @ K x N，kernel 直接沿 signed strides 走訪資料。PR #1208 完成這條正確性基準，也立刻暴露效能牆：大型 C-contiguous GEVM/GEMV 仍跑 scalar loop，曾比 NumPy 慢 4 到 14 倍。位址推導見 Generic。

BLAS

規則 layout 不該留在 scalar loop。若一個 core axis 是 unit stride，另一軸可作合法的 positive leading dimension，BLAS 就能直接讀取這個 view。

PR #1209 以 forced generic / CBLAS 對照量 crossover，再加入 direct BLAS。大型 compatible layouts 回到與 NumPy 相近；小 GEMM 在 min(M,N,K) < 8 時仍保留 generic。

但 core strides (256,-1) 的 reversed LHS 仍進不了 BLAS，step-two view 也一樣。只能每次走 generic，或先排成 BLAS 能讀的 layout。

Pack Once

關鍵是在哪裡 copy。Plan 已從 broadcasting 知道 operand 會被使用 B 次，因此可在 traversal 前把 B * T_generic 改成 T_pack + B * T_blas_call。

同一個 reversed matrix 走十次 generic contraction，或 pack 一次後走十次 BLAS 的成本比較 — `S=32` 的實際成本帳。一次整理 8 KiB LHS，十個 GEMM 共用 packed storage。

PR #1256 把這個判斷放進 selector：只有無法直接交給 BLAS、又由 shape broadcasting 證明會重用的 matrix 才可能 pack。Packing 完成後重建 plan，整個 batch 再走既有 BLAS route。 Dispatch 列出完整 selector。

Packing 本身是成熟技巧。這裡的關鍵是 array-level plan 在 traversal 前已同時知道 layout 與 reuse，因此只 materialize 必要的 supplied operands，而且每份只做一次。

Results

正式 benchmark 把 output batch 固定為十個 contractions，交叉測試 float32 / float64、LHS / RHS / cross broadcasting、negative / step-two inner strides，以及 S=8 到 S=256。S >= 16 時選擇 pack once；72 個因此改變 route 的 cases 全部比舊 generic route 與 NumPy 快。

reversed broadcast matrix 在不同尺寸的 NumPy 相對 solvcon 速度比 — PR #1256，float64、LHS broadcast、negative inner stride。Packing 在 `S=16` 啟用，之後六個尺寸皆由 solvcon 勝出。

在完整的 72 個 affected cases 中，相對舊 route 的 speedup 是 3.30x 到 549.52x，相對 NumPy 2.5.1 是 1.11x 到 2.97x。圖中的 S=32 是 0.0100 ms 對 NumPy 的 0.0244 ms，快 2.44x；S=256 則是 1.6779 ms 對 1.9245 ms，大矩陣同樣勝出。

S=8 與 S=12 保留 generic，也仍由 NumPy 領先。HPC 在這裡不是永遠 zero-copy，也不是永遠先 copy；是讓 shape 揭露 reuse、stride 限定可用 kernel，再由量測決定一次資料整理能不能回本。量測方法與 raw data 收錄在 Benchmarks。

To be continued

目前 Python 的 opt-in entry point 是 matmul_planned()；matmul() 與 @ 尚未切換。一般 API 的整合與後續 routes 仍在收尾。