稀疏矩陣運算 練習題 (Practice - Sparse Matrices and Basic Storage Formats (COO and CSR))

SpMV = A·X + Y 的 sparse matrix-vector multiply + accumulate (A 稀疏，X/Y 為 dense)。每個 nonzero 只做 2 ops (1 乘 1 加)，卻要讀 value(4B)+colIdx(4B)=8B 且矩陣無 data reuse → OP/B ≈ 0.25，屬 memory-bandwidth-bound，故 FLOPS 偏低。

反矩陣運算會產生大量額外 nonzero (fill-ins)，使反矩陣往往比原稀疏矩陣大很多，計算與儲存皆不切實際。改用 iterative method：當 A 為 positive-definite 時可用 conjugate gradient，反覆計算 A·X + Y 並用梯度修正 X 直到收斂。

三陣列：rowIdx[]、colIdx[]、value[] (各 nnz 長度)。平行單位是一個 thread 對應一個 nonzero。因為同一 row 的多個 nonzero 被分配到不同 thread，會同時更新同一個 y[row] → race condition → 必須 atomicAdd(&y[row], x[col]*value)。

用 rowPtrs[] 取代 rowIdx：rowPtrs[r] 是第 r 個 row 的 nonzero 在 value/colIdx 的起始 offset。平行單位是一個 thread 對應一個 row，每 thread 用區域 sum 累加該 row 的 dot product，最後寫入自己獨佔的 y[row] → 不同 thread 寫不同輸出 → 無需 atomic。

長度為 numRows + 1。rowPtrs[numRows] 儲存一個「不存在 row」的起點，其值 = 總 nonzero 數。它讓 kernel 能用 rowPtrs[row+1] 統一界定每一 row (含最後一 row) 的結尾 (exclusive)，無需特例處理。

兩步驟：(1) Padding — 找出 nonzero 最多的 row，其餘 row 尾端補 0 到相同長度成為矩形 (colIdx 一併補)；(2) Transposition — 以 column-major 排放 (等同轉置)，不再需要 rowPtrs。索引公式：i = t*numRows + row。

ELL 既能「給 row → 取該 row 全部 nonzero」(如 CSR)，又能「給 nonzero index i → 反推其 row/col」(如 COO)。col 直接由 colIdx[i] 取得；row 由 i = t*numRows + row 反推為 row = i % numRows (因 row < numRows，故 row % numRows == row)。

當少數 row 的 nonzero 數異常多時，所有 row 都被 pad 到該最大長度 → padding 爆量。Hybrid 作法：轉 ELL 前，把超長 row 尾端的多餘 nonzero 抽出放入獨立 COO，剩餘部分做 ELL，最後兩段結果相加。同時降低 ELL 的 padding 空間開銷與 control divergence (所有 thread 迭代次數變接近)。

nonzeros (row-major)：(0,0)=1,(0,2)=7,(1,2)=8,(2,1)=4,(2,2)=3,(3,0)=2,(3,3)=1。
COO：value=[1,7,8,4,3,2,1]，colIdx=[0,2,2,1,2,0,3]，rowIdx=[0,0,1,2,2,3,3]。
CSR：value、colIdx 同上；每 row nnz = 2,1,2,2 → rowPtrs=[0,2,3,5,7] (長度 numRows+1=5，末元素 7 = 總 nnz)。

設每列最大 nonzero 數為 k。
COO = 3z (rowIdx+colIdx+value)。CSR = 2z + (m+1) (value+colIdx + rowPtrs)。
ELL = 2·m·k (value+colIdx 皆 pad 成 m×k)；需 k，無法只由 m,n,z 求得 (依分佈而定)。
JDS = 2z + (k+1) + m (value+colIdx + iterPtr + row 陣列)；同樣缺 k (最長列長度)。

(1) Histogram：掃描 rowIdx[] 數出每個 row 的 nonzero 數 (per-row count)。(2) Prefix sum / scan：對 per-row counts 做 exclusive scan 得到 rowPtrs[] (各 row 的起始 offset)。之後依 row 重排 (或散佈) value/colIdx 即完成。

COO 中連續 thread 處理連續 nonzero (t0→value[0], t1→value[1]…)，同一時刻存取相鄰位址 → coalesced。CSR 中連續 thread 各負責一整 row，第一次迭代讀各 row 的起點 (value[0],[2],[5],[7])，彼此相距遠 → strided，一個 warp 的請求被拆成多個 memory transaction，浪費頻寬。由於 SpMV 是 memory-bound，這直接拖慢 CSR；ELL/JDS 用 column-major 排放正是為了把 CSR 的 strided 存取轉回 coalesced。

ELL：用 padding+transpose 取得 coalescing，但每 thread 仍依自己 row 的 nnz 迭代 → divergence 未解、空間因 padding 變差。Hybrid：抽走長尾入 COO → padding 大減、divergence 降低，代價是溢位列的 accessibility 變差。JDS：依列長排序使同 warp 列長相近 → 有效降低 divergence 且免 padding (空間更省)，但最不易增刪、且無 padding 故起點無法強制對齊 (coalescing 有時略遜 ELL)。因 SpMV 效能是 data-dependent——取決於 nonzero 分佈 (例如是否有少數超長列)——所以最佳格式須視資料而定，無單一通用最優。

核心訊息：compaction (壓縮去零) ↔ regularization (規則化) 的權衡；SpMV 為 memory-bound (OP/B≈0.25)；最佳格式 data-dependent。