淑樺說得好
你放不下的 別人已經進去了
(X
哈哈哈 今天還是ai slob 抱歉 昨天根莖葉
昨天跟今天的稿子還沒修好 麻煩先擔待一下
2026/09/08 有小更新了 ++了一些技術細節
昨天我們看的是 single-block FWT。

如果 N <= 15，我們還可以很激進地把整個 chunk 放在一個 CUDA block 的 register/shared memory 範圍裡處理。這樣 global memory 只進出一次，中間大部分 butterfly 都在 register 或 warp shuffle 裡完成。

但 N 再往上就不行了。

N=16 的 vector 長度是：

2^16 = 65536 floats

FP32 大小是：

65536 * 4 bytes = 256 KB

這已經不可能靠一個 block 吃完整個 transform。

所以今天要切第二段：large-N multi-pass。

不能 single-block 之後，問題變成資料重排

FWT 的 butterfly stage 有一個很麻煩的地方：stage 越後面，pair 距離越遠。

stage 0:  距離 1
stage 1:  距離 2
stage 2:  距離 4
...
stage s:  距離 2^s

小 stride 可以在 thread 或 warp 裡處理。

中等 stride 可以靠 shared memory。

但大 stride 會跨 block。

一旦跨 block，就不可能用 __syncthreads() 解決。CUDA block 之間沒有一般意義上的同步，也不能互相直接交換 register。

所以做 large-N FWT 時，真正的問題變成：

怎麼把原本一條長 vector 的 butterfly，
拆成幾個可以由 block 各自完成的 pass？

Kronecker split：把一條 vector 看成矩陣

Hadamard 有一個很方便的結構：

H_N = H_a ⊗ H_b

如果 N = a + b，那長度 2^N 的 vector 可以看成一個矩陣：

rows = 2^a
cols = 2^b

接著整個 transform 可以拆成兩段：

先對每一段連續資料做 H_b
再沿著另一個維度做 H_a

以 N=20 為例，如果切成 15 + 5：

把長度 2^20 的向量看成
    32 列  ×  32768 行

Pass 1：每一個長度 32768 的連續 chunk，做完 15 層 FWT
Pass 2：每一條 column（高度 32）再做完剩下 5 層

Pass 1 的資料是連續的，昨天的 single-block kernel 可以直接拿來用。
Pass 2 的 pair 住在不同 chunk 裡，這才是今天要寫的 kernel。

切法會直接決定 Pass 2 有多重。如果第一段只吃 N=10，N=20 會剩下 10 層跨 chunk，column 高度是 1024，單一 warp 放不下。如果第一段吃到 N=15，剩下 5 層，column 高度剛好 32，一個 warp。

所以昨天把 N=15 做成 single-block，意義不只是 N=15 自己變快。它也決定 N=16 到 N=20 的第二段能不能走 warp shuffle。

Pass 1：連續 chunk，直接複用昨天的 kernel

Pass 1 的 layout 最友善。每個 chunk 是連續的 2^{N_INTRA} 個 float。N_INTRA = 15 時，chunk 是 32768 個 float，一個 block 處理一個 chunk，float4 load/store 很自然。

這段沒有新的 butterfly 數學。它就是把昨天的 kernel 對每個 chunk 各跑一次。帳單在兩次 pass 之間：Pass 1 寫回 global memory，Pass 2 再讀回來。

Pass 2：column 的 stride 才是帳單

把 vector 看成 rows × cols 之後，Pass 2 要對每個 column 做 FWT。同一個 column 裡，相鄰兩個 element 在 memory 的距離是 cols，也就是 2^{N_INTRA}。

N_INTRA = 15 時，這個 stride 是 32768 個 float，等於 128 KB。naive 地按 column 去讀，一個 warp 的 32 個 lane 會打到 32 個相距 128 KB 的地址。那不是 coalesced access，是在懲罰 memory system。

所以 Pass 2 要同時處理兩件事：

1. butterfly 本身怎麼做（warp 裡，還是 shared memory 裡）
2. 從 global memory 把資料撈進來時，transaction 長什麼樣子

下面兩顆 kernel，就是這兩件事的兩種答案。

N_INTER <= 5：一個 warp 吃一條 column

剩下層數不超過 5 時，column 高度最多 32。剛好一個 warp。

這時可以讓 256 個 thread 的一個 block，一次處理好幾條 column：

constexpr int ROWS = 1 << ROW_QUBITS;   // 剩下 5 層時 = 32
constexpr int COLS_PER_BLOCK = 256 / ROWS;

int col_in_block = threadIdx.x / ROWS;
int row_idx      = threadIdx.x % ROWS;
int col_idx      = blockIdx.x * COLS_PER_BLOCK + col_in_block;

float val = in[row_idx * cols + col_idx];

#pragma unroll
for (int stage = 0; stage < ROW_QUBITS; ++stage) {
    int stride = 1 << stage;
    float peer = __shfl_xor_sync(0xffffffff, val, stride);
    if ((row_idx & stride) == 0)
        val = val + peer;
    else
        val = peer - val;
}

out[row_idx * cols + col_idx] = val;

butterfly 這一段很乾淨：沒有 shared memory，沒有 bank conflict，exchange 全在 register。N=16 剩下 1 層、N=20 剩下 5 層，都可以走這條。

但它不是免費的。val 從哪裡來，才是成本。

ROW_QUBITS = 5 時，一個 warp 的 32 個 lane 對應同一條 column 的 32 個 row。它們讀的地址相差 cols。block 內計算再便宜，這筆 strided load/store 還是要付。

所以這顆 kernel 的設計目標很明確：把 block 內交換 壓到最低，承認 global memory 那一趟 避不掉。

N_INTER > 5：tiled shared memory，而且要向量化

剩下超過 5 層，column 高度超過 32，單一 warp 放不下。這時只好一次搬 TILE_COLS 條 column 進 shared memory，在 tile 裡做 row 方向的 butterfly，再寫回去。

舊的寫法 tile 太窄，一次只搬 16 條 column，bandwidth 很難打滿。後來改成一次 64 條，並且 load/store 用 float4：

constexpr int TILE_COLS = 64;
constexpr int SMEM_PITCH = TILE_COLS + 4;  // padding，避開 bank conflict

float4 v;
asm volatile("ld.global.nc.v4.f32 {%0,%1,%2,%3}, [%4];"
             : "=f"(v.x), "=f"(v.y), "=f"(v.z), "=f"(v.w)
             : "l"(ptr));

smem[r * SMEM_PITCH + c + 0] = v.x;
smem[r * SMEM_PITCH + c + 1] = v.y;
smem[r * SMEM_PITCH + c + 2] = v.z;
smem[r * SMEM_PITCH + c + 3] = v.w;

ld.global.nc 是 non-coherent load：這份資料用過就不會很快再讀，告訴 L1 不必當熱資料守著。

SMEM_PITCH = 64 + 4 則是為了 bank。shared memory 通常 32 個 bank。如果 pitch 剛好是 64，同一條 column 走下去很容易撞同一批 bank。多 pad 4 個 float，是讓相鄰 row 的同一條 column 錯開。

butterfly 本身還是最普通的 pairwise add/sub，只是改在 padded tile 上做：

float a = smem[r_i * SMEM_PITCH + c];
float b = smem[r_j * SMEM_PITCH + c];
smem[r_i * SMEM_PITCH + c] = a + b;
smem[r_j * SMEM_PITCH + c] = a - b;

這條路的 tuning 點不是 Tensor Core：

TILE_COLS 太小 → 一次搬的連續資料不夠，bandwidth 打不滿
TILE_COLS 太大 → shared memory 壓力上升，occupancy 掉
pitch 沒補好   → bank conflict
load 不連續    → sector efficiency 下降

large-N FWT 難寫，就是難在這裡。不是加減法算不動，是資料不再自然連續。

為什麼 N=16 會突然變慢

同一組 FP32、batch=128 的數字：

N=15 :  82.98 us    （還在 single-block）
N=16 : 584.56 us    （開始 2-pass）

不是多一層加減就突然慢七倍。是 N=16 開始必須把資料寫出再讀回。

state size = 128 × 65536 × 4 bytes  ≈ 33.5 MB
2-pass，每 pass 讀一次寫一次
traffic  ≈ 33.5 MB × 4  ≈ 134 MB

134 MB / 584 us  ≈ 229 GB/s

這已經很接近在搬資料，不是在等 ALU。33.5 MB 也超過很多筆電 GPU 的 L2，所以不能幻想「反正 cache 會幫你留著」。

所以 large-N 的下一步不是再塞更多計算，而是：少一個 pass、讓每次 load/store 更連續、把 corner turning 變便宜。

正確性也要分 pass 驗

小 N 對了，只代表 single-block kernel 對。large-N 還要驗 Pass 1 寫出去的 layout、Pass 2 的 column mapping、normalization 是最後乘一次還是每 pass 各乘一次、batch stride，以及 in-place 有沒有覆蓋。

我們踩過一個很典型的錯：reinterpret_cast<float4*> 之後，pointer 算術還用 float 的 index 去加。float4* 加 1 會跳 16 bytes，不是 4 bytes。結果 batch 之間 memory overlap。kernel 照跑、速度看起來也正常，output 已經被寫穿。

所以驗收要兩層：小 N 對 naive FWT 做 elementwise；中大 N 不能只看 norm 有沒有守住。norm 太寬了，蓋不住這種 overlap。

今天切出來的形狀

N <= 15
    一個 block 做完

N > 15
    Pass 1：每個連續 chunk 做完 N=15
    Pass 2：
        剩下 <= 5 層 → warp shuffle
        剩下  > 5 層 → float4 + padded shared memory tile

能不能把第一段推到 15，會直接決定第二段有多重。N_INTER <= 5 時，block 內交換可以幾乎免費；超過就開始跟 coalescing、bank conflict、tile 寬度打架。

下一篇回頭看 Tensor Core。Hadamard 只有 +1 / -1，看起來很適合 MMA。問題是：為了把它改寫成 GEMM，你得多搬幾次資料。

會贏的喔
新宿再次響起了五條悟的吟唱。“九綱”“偏光”“烏與聲明”“表裏之間”、宿儺明白自己再也沒有任何機會阻止茈的誕生了....
會贏的喔

(X
明天見~