SELL-C-σ 和 SELL-P

Sparse Matrix Aug 15, 2021

在前面"稀疏矩陣 (Sparse Matrix) 介紹"有介紹了 DENSE、CSR、COO、ELL 這些比較常見的矩陣類型，而這一篇這將會介紹 ELL 的改進版， SELL-C-σ 以及 SELL-P。

本篇主要基於下面兩篇再加上我對於GPU的理解來敘述這兩種矩陣差異：

SELL-C-σ: Kreutzer, M., Hager, G., Wellein, G., Fehske, H., & Bishop, A. R. (2014). A unified sparse matrix data format for efficient general sparse matrix-vector multiply on modern processors with wide SIMD units
SELL-P: Anzt, H., Tomov, S., & Dongarra, J. (2014). Implementing a Sparse Matrix Vector Product for the SELL-C/SELL-C-σ formats on NVIDIA GPUs

ELL(PACK) 的優缺點

ELL 的優點是在於它的矩陣向量乘其實很合乎 SIMD 或 SPMD 的架構，可以很直覺地將每一個 thread (執行緒) 分配給每一列，記憶體上的讀取也都是連續地一組接著一組，但是為了合乎這個條件，ELL需要將每一列的儲存空間填到一樣長。

假設，原先矩陣有一列有 10 個 nonzeros (非零元素)，但其他列都只有 1 ~ 3 個 nonzeros，但 ELL 就得要將每一列都填滿 5 個，因此造成空間的浪費，以及其實計算到後面時，雖然還是符合良好的記憶體讀取，但只有其中一列有實際在工作，這樣效率其實也會很低。

當然，如果矩陣每一列都有著差不多的數量，那 ELL 就十分合適，但如果沒這麼合適的話有沒有其他方式去利用 ELL 優點，但卻又減少他缺點的影響呢？

那就是引入 sliced (切片) 的概念。將每幾列視為一小組，每一小組再用 ELL 去處理，因為只考慮一小組，所以某些列的影響只會限縮所在該組中，減輕 ELL 的缺點，而又由於 SIMD 有寬度 (simd-width) 上限或者 GPU 的 thread-block 或者 warp 也有一定的使用範圍，所以每一小組要含多少列，其實也會跟這些值有關。

下圖可以看到原先 ELL 會讓每一列都需要填滿 5 個，但換成 SELL-2 (每兩列一組) ，最長的列只會影響到他所在的那一組。由於切片的緣故，每一組的長度不唯一了，存法也有點不太一樣，是將一組的 ELL 存完後，再存下一組的 ELL，不能像 ELL 只存 value (值) 以及 col_idx (行座標)，還需要紀錄每一組的期始點 (block_row_ptr) ，每一組的結束點就可以利用下組的起始點來使用，而最後要額外紀錄結束點，說到這裡大家可能會覺得有點熟悉，他跟 CSR 的 row_ptr 概念一樣，只是現在紀錄的不是每一列的起始點，而是每一組的起始點。 Sliced ELL 也有需要付出的代價，當原先矩陣的列數不能被[切片所需列數]來整除的話，就要補上所缺的列數，或者要在程式內多加判斷。

SELL-C-σ

S 是指稱 Sliced
ELL 是指 ELLPACK
C 是 chunk-size (每幾列作為一組)
σ 是 sorting-space (每幾列視為同一組來排序)

C 這邊給定的就是 slice (切片) 的大小，那 σ 則是一次要將幾列排序，而排序方法是將非零元素多的往前排，而這樣讓相鄰的列有較相近的工作量要處理，如果 σ 是 1 就等同於不排序的狀況，而 SELL-C-σ 該篇中是直接可以用數字來表示相對應的值。例如： SELL-2-4 就是每 2 列為一組並對於每 4 列來排序。

下圖是範例， SELL-C (σ = 1) 的部分在前面有介紹了，未提及 σ 參數的話就是一般未排序的狀況，而最右邊則是以每四列下去做排序，那由於將數量多的往前排，所以浪費的空間也相對比較少，但這裡就需要多存一個列的排序資料，以便描述原始矩陣，SELL-C-σ 比 SELL-C 多考慮排序的狀況，其他並無太大差別

SELL-P

S 是指稱 Sliced
ELL 是指 ELLPACK
P 是 padding (每列的單位長度，需要補零至其倍數)

該篇當中是以 t 來表示每列的單位長度，需要補零至其倍數，當 t = 2 時，即使某一組最長的是 5，最終長度會使用 6 來符合 t = 2 (長度表示非零元素個數)

\[ 每一組的長度 = \lceil\frac{該組最長長度}{t}\rceil \times t \]

這裡就不像 SELL-C-σ 有直接用數字放在相對應的位置來表示，會直接額外寫出參數該用甚麼，例如 b = block size (slice 的大小) 以及剛剛說到的 t，其實這邊也給出了每一組的最小組成單位 b * t。

下圖是 SELL-P 的範例，在 t = 2 的情形下，當最長的長度為奇數時，SELL-P 就會補零 (padding) 補到湊滿偶數，也可以看到綠色方框 (b * t = 2 * 2) 為最小組成單位。SELL-P 則沒有考慮排序進去，所以跟 SELL-C 的差別，就是多存了一些零以及 block_row_ptr 會因此不一樣，如果剛好每組最長的都可以被 t 整除的話，那就會一樣了。

那為何這邊要多一個參數補零呢？因為它不單只有一條 thread 在同一列，補零到某個值的倍數，讓多個 thread 在同一列也還是有同樣分量的工作，減少分支的情形發生。

下圖是 SELLP 怎麼計算的，不會只安排一條 thread 在同一列，每次會推移一個 thread-block (或 warp) 將部分加總起來 (partial sum)，最後再把同一列的 partial sum 加起來輸出，補零的功用就可以在一些末端時，不用額外再用一些判斷去處理。而 b * t 的設計就要根據硬體架構下去設計，例如如果是以 warp 來作為單位的話可能就會取 b * t = 32 等

總結

此篇介紹了不同種的 Sliced ELL，這些都有一些參數需要根據硬體、程式設計或者矩陣類型來進行調整，Sliced ELL 利用了 CSR 的 row_ptr 概念紀錄 block_row，將原先矩陣進行切片來利用 ELL 的優點，並減少 ELL 在一些狀況下的缺點。最後一張圖大略描述他們怎麼儲存矩陣。