夏 天 付格林 曲劭儒 羅中沛 任鵬舉

1（人機(jī)混合增強(qiáng)智能全國(guó)重點(diǎn)實(shí)驗(yàn)室（西安交通大學(xué)） 西安 710049）

2（視覺信息與應(yīng)用國(guó)家工程研究中心（西安交通大學(xué)） 西安 710049）

3（西安交通大學(xué)人工智能與機(jī)器人研究所 西安 710049）

（pengjuren@xjtu.edu.cn）

在科學(xué)計(jì)算、大數(shù)據(jù)分析和智能計(jì)算任務(wù)中，稀疏矩陣向量乘法（sparse matrix-vector multiplication，SpMV）是廣泛使用的計(jì)算核心.例如，在大數(shù)據(jù)推薦系統(tǒng)[1]和圖神經(jīng)網(wǎng)絡(luò)[2-4]中，通常會(huì)對(duì)圖鄰接矩陣進(jìn)行迭代的矩陣向量乘法來完成大規(guī)模的節(jié)點(diǎn)信息匯總和傳播，如在經(jīng)典的Pagerank 算法中，不同網(wǎng)頁之前的關(guān)聯(lián)度用稀疏矩陣表示，通過對(duì)稀疏矩陣做多次SpMV 迭代計(jì)算，能求解不同網(wǎng)頁的排序關(guān)系；而在現(xiàn)代科學(xué)計(jì)算任務(wù)的數(shù)值模擬計(jì)算中，也經(jīng)常使用迭代法來進(jìn)行大規(guī)模稀疏線性方程組的求解.在以上計(jì)算任務(wù)中，SpMV 構(gòu)成了主要的計(jì)算負(fù)載，其計(jì)算性能直接影響了整體算法的解算效率.因此，提升SpMV 的計(jì)算性能對(duì)于優(yōu)化許多科學(xué)計(jì)算和智能計(jì)算任務(wù)性能具有重要意義.

與此同時(shí)，基于超標(biāo)量和亂序技術(shù)的現(xiàn)代處理器架構(gòu)高度依賴計(jì)算與訪存模式的規(guī)則化和可預(yù)測(cè)性來充分發(fā)揮其計(jì)算性能.具體來說，現(xiàn)代處理器通常依賴準(zhǔn)確的分支預(yù)測(cè)技術(shù)來提升其流水線執(zhí)行效率，并且采用高效的數(shù)據(jù)預(yù)取技術(shù)來有效掩蓋訪存延遲，分支預(yù)測(cè)技術(shù)和數(shù)據(jù)預(yù)取技術(shù)均需要計(jì)算程序在控制跳轉(zhuǎn)和數(shù)據(jù)訪問行為上呈現(xiàn)高度的可預(yù)測(cè)性.此外，數(shù)據(jù)訪問的時(shí)間局部性與空間局部性直接影響緩存的使用效率，對(duì)處理器性能也至關(guān)重要.然而在SpMV 中，稀疏矩陣通常具有極高的稀疏率，非零元素的占比通常低于0.01%，并且其分布通常是隨機(jī)的，并無規(guī)律可循.這就導(dǎo)致在計(jì)算過程中對(duì)于向量元素的索引和訪問呈現(xiàn)出高度不規(guī)則性，缺乏訪存局部性，嚴(yán)重制約了緩存和訪存帶寬的使用效率[5].與此同時(shí)，由于矩陣中每行的非零元素?cái)?shù)量不同，也導(dǎo)致了分支預(yù)測(cè)和并行化計(jì)算的困難.這些因素導(dǎo)致了SpMV 在很多算法中構(gòu)成了嚴(yán)重的性能瓶頸.因此，SpMV 在現(xiàn)代處理器上的設(shè)計(jì)與優(yōu)化至今仍然是研究熱點(diǎn)和亟需解決的問題.

目前，對(duì)于SpMV 的研究成果已有很多，其中很多工作注重提升計(jì)算并行度以適配SIMD（singleinstruction-multiple-data）和SIMT 等并行計(jì)算技術(shù)[6-9].還有一些工作通過進(jìn)行矩陣的分塊計(jì)算來提升數(shù)據(jù)的局部性[10-14].然而，這些工作只局限在某個(gè)單一維度上，而缺乏對(duì)于整體算法的協(xié)同優(yōu)化，因此可能造成整體的優(yōu)化效果受到制約.例如，基于分塊方法的數(shù)據(jù)局部性提升經(jīng)常會(huì)導(dǎo)致程序復(fù)雜度的上升，可能會(huì)帶來指令數(shù)的增加，從而導(dǎo)致核心的計(jì)算指令密度下降、計(jì)算效率被削弱.此外，現(xiàn)有研究工作中對(duì)于分支預(yù)測(cè)和數(shù)據(jù)預(yù)取等可預(yù)測(cè)的計(jì)算特性如何影響性能仍然缺乏充分的討論和優(yōu)化.

本文通過研究，總結(jié)出制約SpMV 性能的3 個(gè)關(guān)鍵因素：高訪存延遲、低分支預(yù)測(cè)準(zhǔn)確率以及低代碼計(jì)算密度.針對(duì)這些問題，本文構(gòu)建了一種新型的SpMV 計(jì)算方法——pvSpMV（predictable vectorized SpMV）.該方法通過對(duì)稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)重組以及向量元素的位置重組，實(shí)現(xiàn)了程序控制流和訪存行為的高預(yù)測(cè)性，從而充分發(fā)掘亂序處理器流水線和數(shù)據(jù)預(yù)取的執(zhí)行效率；在此基礎(chǔ)上，通過對(duì)程序復(fù)雜度的降低和向量化程度的提升，顯著優(yōu)化代碼計(jì)算密度，從而實(shí)現(xiàn)對(duì)于SpMV 的全面提升和性能優(yōu)化.本文還基于SuiteSparse[15]稀疏矩陣測(cè)試集評(píng)估pvSpMV 性能，測(cè)試結(jié)果顯示，pvSpMV 可以在主流Intel 高性能處理器上獲得與主流商用計(jì)算庫MKL相比平均2.6 倍的加速比，相比于現(xiàn)有最先進(jìn)算法獲得平均1.3 倍的加速比.此外，pvSpMV 所需的預(yù)處理開銷較低，預(yù)期可以在迭代計(jì)算中被有效分?jǐn)?，從而不影響整體算法的性能提升.

本文的主要貢獻(xiàn)包括4 個(gè)方面：

1）構(gòu)建一種訪存模式可預(yù)測(cè)的SpMV 表示方法pvSpMV，可以高效適配數(shù)據(jù)預(yù)取器，大幅度提升取數(shù)帶寬，并且提升緩存的利用效率；

2）在此基礎(chǔ)上，分析SpMV 中的關(guān)鍵分支跳轉(zhuǎn)指令，并且提出一種分支跳轉(zhuǎn)可預(yù)測(cè)的稀疏矩陣數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)，提升處理器的執(zhí)行效率；

3）定義用于衡量程序效率的有效計(jì)算指令密度的性能指標(biāo)，并且以此指標(biāo)為基礎(chǔ)，設(shè)計(jì)一種采用極低程序復(fù)雜度完成的高效SIMD 計(jì)算方法，提升有效計(jì)算指令密度；

4）采用來自多種領(lǐng)域的稀疏矩陣開展了充分實(shí)驗(yàn)，并且與現(xiàn)有的先進(jìn)方法進(jìn)行對(duì)比，結(jié)果表明pvSpMV 可以有效緩解SpMV 問題中的多個(gè)關(guān)鍵性能瓶頸，具有較好的性能收益.

1 基礎(chǔ)知識(shí)

1.1 SpMV 概述

SpMV 的計(jì)算過程可以表示為等式Ax=y，其中A為參加運(yùn)算的稀疏矩陣，x為參加運(yùn)算的稠密向量，y為計(jì)算結(jié)果向量.在迭代計(jì)算中，A為方陣，每次迭代的x向量為上一次迭代中生成的y向量.在實(shí)際應(yīng)用中，矩陣A通常呈現(xiàn)出極高的稀疏性，為了避免存儲(chǔ)大量冗余的零元素，矩陣A通常采用稀疏壓縮格式進(jìn)行存儲(chǔ).常用的稀疏矩陣壓縮格式包括COO（coordinate list）,ELLPACK[16],CSR（column sparse row）,CSX[17],BRC[18]等.其中，CSR 格式是目前在科學(xué)計(jì)算領(lǐng)域使用最為廣泛的壓縮格式，本文基于這一格式進(jìn)行研究和設(shè)計(jì).

圖1 中給出基于CSR 格式的稀疏矩陣存儲(chǔ)案例以及對(duì)應(yīng)的SpMV 方法.可以看到，CSR 格式保存3個(gè)數(shù)組結(jié)構(gòu)，分別為：nnz用于保存所有非零元素值，colidx用于保存左右非零元素的列坐標(biāo)，rowptr用于記錄矩陣中每行非零元素的起始位置.在SpMV 的計(jì)算過程中，程序按照rowptr數(shù)組的索引，遍歷nnz數(shù)組中存儲(chǔ)的矩陣中每行非零元素值，并且與colidx索引的對(duì)應(yīng)x向量元素進(jìn)行乘累加操作，最終生成y向量的元素.對(duì)于x向量訪問遵循x[colidx[i]]的間接索引模式.由于稀疏矩陣A的非零元素通常隨機(jī)分布，所以x的訪存地址也是完全隨機(jī)的.

Fig.1 SpMV algorithm using CSR format圖1 基于CSR 格式的SpMV 算法

1.2 SpMV 相關(guān)工作

由于稀疏矩陣非零元素分布的不規(guī)則性，如何在SIMD 架構(gòu)的現(xiàn)代處理器上高效運(yùn)行SpMV 一直是研究熱點(diǎn).Liu 等人[7]基于ELLPACK 格式提出名為ESB 的編碼格式，通過列分塊改善訪存的空間局部性，提升SIMD 執(zhí)行效率.為了充分利用SIMD 單元，CSR5[8]和CVR[9]算法用矩陣多行的非零元素填充SIMD 單元，提升SIMD 單元執(zhí)行效率，但CSR5和CVR 需要額外復(fù)雜的預(yù)處理和后處理器操作，CSR5 需要加法樹設(shè)計(jì)與部分和存儲(chǔ)，CVR 需要在每次向量計(jì)算后動(dòng)態(tài)處理向量單元不同元素的輸出，這些額外的數(shù)據(jù)結(jié)構(gòu)和復(fù)雜的計(jì)算模式不利于計(jì)算效率的提升，因此如何設(shè)計(jì)低復(fù)雜度且高計(jì)算效率的SpMV 成為主要的研究挑戰(zhàn).

1.3 SpMV 的性能分析

為了進(jìn)一步分析SpMV 在現(xiàn)代處理器中的計(jì)算行為和性能瓶頸，圖2 列出一個(gè)典型SpMV 算法的標(biāo)量代碼實(shí)現(xiàn)以及在ARMv8 GCC 編譯器下生成的對(duì)應(yīng)的指令流示例，并且標(biāo)出其中的訪存指令、計(jì)算指令和跳轉(zhuǎn)指令.在圖3 中，我們使用Intel VTune Profiler工具測(cè)試多個(gè)不同領(lǐng)域和尺寸稀疏矩陣的SpMV 性能，并且記錄其在Intel Xeon Gold 6146 處理器上的性能分解.從圖3（a）中可以看到，訪存開銷和分支開銷構(gòu)成主要性能瓶頸.本節(jié)將逐一分析SpMV 在訪存、分支和計(jì)算效率上的性能問題.

Fig.2 Pseudo code of SpMV and instruction flow in ARMv8圖2 SpMV 偽代碼與ARMv8 指令流

Fig.3 Performance bottleneck analysis of SpMV on Intel Xeon CPU圖3 Intel Xeon CPU 上SpMV 的性能瓶頸分析

首先，現(xiàn)代處理器依賴數(shù)據(jù)局部性來提升緩存效率和降低訪存延遲，并且通過數(shù)據(jù)預(yù)取技術(shù)來進(jìn)一步隱藏訪存開銷.通常現(xiàn)代高性能商業(yè)處理器會(huì)在緩存配置硬件數(shù)據(jù)預(yù)取器，例如常見的鄰行預(yù)取器（next-line prefetcher）和步長(zhǎng)預(yù)取器（stride prefetcher）[19].這些預(yù)取器可以動(dòng)態(tài)檢測(cè)流式訪問和定步長(zhǎng)訪問等簡(jiǎn)單的訪存規(guī)律并進(jìn)行數(shù)據(jù)預(yù)取.

從圖2 中可知，SpMV 中對(duì)于矩陣數(shù)據(jù)（rowptr,colidx,nnz）和向量y的訪問均是順序串行的，因此可以高效地被步長(zhǎng)預(yù)取器識(shí)別并預(yù)取.相反地，對(duì)于向量x的訪問呈現(xiàn)出隨機(jī)的間接索引模式，通常無法被有效地識(shí)別和預(yù)取[20].此外，由于向量x訪問地址的隨機(jī)性，數(shù)據(jù)局部性極差，導(dǎo)致緩存使用效率大幅度下降，緩存缺失頻繁發(fā)生，以上問題均導(dǎo)致訪存成為SpMV 的主要性能瓶頸.圖3（b）和圖3（c）中分別列出了SpMV 的預(yù)取覆蓋率和訪存帶寬利用率.可以看到，由于對(duì)于x的訪問造成大量緩存缺失且無法被有效預(yù)取，導(dǎo)致整個(gè)計(jì)算過程中僅有少量數(shù)據(jù)訪問（不足30%）可以被有效預(yù)取.這也導(dǎo)致整體的訪存帶寬利用率低下，進(jìn)一步加劇訪存性能瓶頸.

其次，現(xiàn)代處理器通常使用分支預(yù)測(cè)技術(shù)來提升流水線的執(zhí)行效率.硬件分支預(yù)測(cè)器通過動(dòng)態(tài)記錄和分析不同分支跳轉(zhuǎn)指令的歷史行為，學(xué)習(xí)其跳轉(zhuǎn)規(guī)律以進(jìn)行預(yù)測(cè).對(duì)于具有較深流水線的高性能處理器，分支預(yù)測(cè)失敗需要清空錯(cuò)誤指令并且回退處理器狀態(tài)，這將嚴(yán)重降低處理器的性能.有研究表明，Intel 架構(gòu)的每次分支預(yù)測(cè)錯(cuò)誤將導(dǎo)致10～15 個(gè)周期的CPU 停滯[21].通常認(rèn)為，分支預(yù)測(cè)錯(cuò)誤率如果超過5%就會(huì)明顯制約處理器性能.由圖2 可以看出，SpMV 的指令流中存在3 個(gè)分支跳轉(zhuǎn)指令，分別用于判斷末尾行（BR1）、空行（BR2）和行遍歷結(jié)束（BR3）.其中，BR1 和BR2 具有規(guī)律性，可以比較準(zhǔn)確地進(jìn)行預(yù)測(cè).然而，BR3 的跳轉(zhuǎn)條件取決于每行的非零元素個(gè)數(shù)，考慮到稀疏矩陣的隨機(jī)分布，BR3 行為毫無規(guī)律，因此難以被正確預(yù)測(cè).因此，BR3 是SpMV 的關(guān)鍵分支跳轉(zhuǎn).圖3（d）列出了不同矩陣中BR3 分支預(yù)測(cè)情況，可以看到BR3 具有很高的失敗概率（20%～40%），這也導(dǎo)致了分支預(yù)測(cè)成為了SpMV 計(jì)算中的關(guān)鍵性能瓶頸之一.

此外，基于現(xiàn)代處理器的SpMV 優(yōu)化技術(shù)通常會(huì)采用SIMD 技術(shù)來提升計(jì)算并行度，從而優(yōu)化計(jì)算性能.SIMD 支持每條指令實(shí)現(xiàn)最多VL個(gè)元素的并行運(yùn)算.然而，由于稀疏矩陣的高度稀疏性和不規(guī)則性，SpMV 運(yùn)算中通常難以實(shí)現(xiàn)VL個(gè)元素的充分并行，因此需要通過增加遮蔽向量（mask vector）的方式降低其并行度.為了衡量SIMD SpMV 算法的運(yùn)算效率，我們定義有效計(jì)算密度（effective computation density）指標(biāo)D和浮點(diǎn)計(jì)算性能（floating point computing performance）指標(biāo)FPCP，其公式為：

其中，VInsnNum表示SpMV 中參與有效計(jì)算的SIMD指令；VRate表示向量利用率，即未被屏蔽的SIMD 向量元素的平均占比；TotalInsnNum表示程序執(zhí)行的指令總數(shù)；IPC表示單位周期數(shù)執(zhí)行指令數(shù)；Frequency表示處理器的執(zhí)行頻率.對(duì)于給定的工作負(fù)載，IPC相對(duì)穩(wěn)定，僅在較小范圍內(nèi)波動(dòng)，因此由式（2）可知，計(jì)算密度D與計(jì)算性能FPCP呈近似正相關(guān)關(guān)系.需要注意的是，這里的有效計(jì)算指的是對(duì)于每個(gè)非零元素的乘加運(yùn)算，因此其數(shù)量是由矩陣的非零元素個(gè)數(shù)決定的.通?；赟IMD 的研究主要關(guān)注向量利用率VRate這一指標(biāo)，然而通過式（1）可知，有效向量計(jì)算指令在整體指令流中的密度也非常重要.在一些算法中，為了追求高利用率而提升程序的復(fù)雜度反而妨礙了整體計(jì)算密度的提升.此外，分支預(yù)測(cè)錯(cuò)誤會(huì)造成大量無效指令被流水線發(fā)射，從而導(dǎo)致整體計(jì)算密度的下降.

圖3（e）比較了Scalar SpMV 算法和2 種SIMD 優(yōu)化算法CSR5[8]和CVR[9].可以看到，通過優(yōu)化，CSR5和CVR 均達(dá)到了很高的向量利用率，且遠(yuǎn)超于Scalar SpMV 算法.然而，通過圖3（f）的計(jì)算密度指標(biāo)對(duì)比可以看到，CSR5 和CVR 的整體計(jì)算效率僅稍好于Scalar SpMV.這是由于CSR5 和CVR 算法為了提升向量利用率，引入了額外的程序復(fù)雜度，從而導(dǎo)致整體指令數(shù)的上升，削減了向量指令的密度，由式（2）可知，計(jì)算指令密度下降會(huì)導(dǎo)致計(jì)算性能下降，不利于SpMV 的性能優(yōu)化.這一實(shí)驗(yàn)結(jié)果說明，為了有效提升計(jì)算性能，不應(yīng)該僅關(guān)注向量利用率，而是應(yīng)綜合考慮向量利用率與程序復(fù)雜度的權(quán)衡，盡可能降低與核心計(jì)算任務(wù)無關(guān)的指令數(shù)量.

通過對(duì)于SpMV 訪存、分支跳轉(zhuǎn)和計(jì)算并行度的分析，明確了核心性能瓶頸和指標(biāo).本文基于這些分析結(jié)果提出了pvSpMV 方法，通過構(gòu)建訪存和分支跳轉(zhuǎn)高度可預(yù)測(cè)的計(jì)算方法，結(jié)合低復(fù)雜度的向量化計(jì)算，實(shí)現(xiàn)了SpMV 計(jì)算方法的全面優(yōu)化.

2 訪存優(yōu)化方法

2.1 步長(zhǎng)預(yù)取器工作原理

步長(zhǎng)預(yù)取器是現(xiàn)代商業(yè)處理器中最常見的硬件預(yù)取器，它通過監(jiān)控處理器發(fā)送給緩存的訪存指令多個(gè)內(nèi)存地址，識(shí)別連續(xù)步長(zhǎng)訪問模式的訪存序列.建立步長(zhǎng)訪存模式后，步長(zhǎng)預(yù)取器的預(yù)取地址生成模塊（address generator）監(jiān)聽處理器發(fā)送給緩存訪存請(qǐng)求，根據(jù)訪問地址計(jì)算預(yù)取數(shù)據(jù)的內(nèi)存地址，計(jì)算公式為：

其中，load addr為觸發(fā)預(yù)取器生成預(yù)取請(qǐng)求的訪存地址；prefetch depth為預(yù)取深度，指一個(gè)訪問向量x的地址產(chǎn)生預(yù)取請(qǐng)求的數(shù)量，為硬件預(yù)取器參數(shù)；stride為訪存步長(zhǎng)，指連續(xù)訪問地址之間的差.預(yù)取器產(chǎn)生預(yù)取請(qǐng)求并通過預(yù)取隊(duì)列（prefetch queue）發(fā)送給下一級(jí)訪存結(jié)構(gòu)，在程序需要這個(gè)請(qǐng)求數(shù)據(jù)時(shí)，數(shù)據(jù)已經(jīng)被放入緩存，達(dá)到隱藏訪存延遲的效果.

在SpMV 運(yùn)算過程中，對(duì)向量x的訪存地址由于非零元分布不規(guī)則且不連續(xù)，因此步長(zhǎng)預(yù)取器不能有效預(yù)取向量x.而傳統(tǒng)的SpMV 優(yōu)化策略一般通過分塊的方式，使得單塊的數(shù)據(jù)尺寸小于緩存尺寸，改善訪存的空間局部性，忽視對(duì)向量非規(guī)則訪問模式的分析和優(yōu)化.本文從硬件步長(zhǎng)預(yù)取器的角度分析SpMV 運(yùn)行過程中的訪存效率問題，通過將非規(guī)則的向量訪問模式序列化，使得步長(zhǎng)預(yù)取器能有效預(yù)取向量，隱藏訪存延遲，優(yōu)化訪存性能.

2.2 向量訪存序列順序化

圖4 以4×8 矩陣為例，描述SpMV 運(yùn)算對(duì)向量x不同元素的訪問過程.B→C→H→A→D→G 是矩陣第1 次訪問向量元素的序列順序，由于緩存未加載向量元素，第1 次向量元素的訪問會(huì)產(chǎn)生緩存未命中，并從內(nèi)存獲取數(shù)據(jù).假設(shè)后續(xù)對(duì)同一元素的訪問由于時(shí)間局部性會(huì)在緩存命中，則訪問向量元素對(duì)應(yīng)的緩存未命中序列（cache miss sequence）為B→C→H→A→D→G，是順序訪問模式，能夠利用步長(zhǎng)預(yù)取器進(jìn)行數(shù)據(jù)預(yù)取.因此以對(duì)向量不同元素的第1 次訪問的先后訪問順序?qū)ο蛄縳進(jìn)行重排，緩存預(yù)取器的預(yù)取模式識(shí)別模塊（pattern detector）通過觀測(cè)向量x對(duì)應(yīng)的緩存未命中序列，識(shí)別向量x的連續(xù)訪存模式.

Fig.4 Illustration of SpMV memory access pattern圖4 SpMV 訪存模式示意圖

2.3 矩陣分塊策略

訪存優(yōu)化方法關(guān)鍵在于步長(zhǎng)預(yù)取器能識(shí)別訪問向量x對(duì)應(yīng)的緩存未命中這一順序序列.若重排后的向量x尺寸過大，導(dǎo)致部分已被訪問的向量元素被緩存替換，就會(huì)影響緩存對(duì)向量元素訪問模式的識(shí)別.因此要求后續(xù)對(duì)相同元素的訪存盡可能在緩存命中，步長(zhǎng)預(yù)取器觀測(cè)對(duì)向量元素的第1 次訪問產(chǎn)生的緩存未命中請(qǐng)求，構(gòu)建連續(xù)順序的訪問序列.

解決方案是從行方向?qū)ο∈杈仃嚽袎K，在塊內(nèi)分別做SpMV 運(yùn)算.切塊思路是讓單個(gè)矩陣塊內(nèi)對(duì)向量x元素的第1 次訪問被緩存的步長(zhǎng)預(yù)取器觀測(cè)，而其余次訪問在緩存命中，并且預(yù)留足夠的緩存空間給SpMV 運(yùn)算要用的CSR 數(shù)據(jù)結(jié)構(gòu)（value,rowptr,colidx），如圖4 所示.通過實(shí)驗(yàn)分析性能結(jié)果，將單個(gè)矩陣塊的重排后的向量x′尺寸設(shè)置為0.5 倍緩存尺寸.

2.4 基于Bitmap 的行重排策略

SpMV 在科學(xué)計(jì)算和圖計(jì)算中需要多次迭代，上一輪迭代的計(jì)算結(jié)果作為下一輪迭代的數(shù)據(jù)輸入.由于每次迭代都需要用本次計(jì)算的輸出向量y更新下一次運(yùn)算要使用的重排后的向量x′，因此如何減少更新向量x′的訪存開銷對(duì)性能優(yōu)化十分重要.

本文采用基于Bitmap 的行重排策略緩解向量x更新開銷過大問題，如圖5 所示.將稀疏矩陣的每一行分為不同的區(qū)域，將含有最多非零元素個(gè)數(shù)的區(qū)域?qū)?yīng)的Bitmap 位標(biāo)注為1，其余位標(biāo)注為0，得到矩陣每行對(duì)應(yīng)的Bitmap 值.并利用Bitmap 值以從大到小的順序?qū)仃囆羞M(jìn)行重排，目的是讓矩陣的非零元素分布更規(guī)律，近似為對(duì)角線分布，這樣不同塊的重排后的向量x′尺寸更小，向量x′的更新開銷減少.

Fig.5 Illustration of Bitmap row reorder strategy圖5 Bitmap 行重排策略示意圖

2.5 訪存優(yōu)化方法介紹

訪存優(yōu)化方法如圖6 所示.

Fig.6 Illustration of memory access optimization method圖6 訪存優(yōu)化方法示意圖

圖6 所示的主要步驟包括：

1）基于Bitmap 的行重排（Bitmap row reordering）.讓矩陣的非零元素呈近似對(duì)角線分布.

2）矩陣分塊（block partitioning）.從行方向?qū)⒕仃嚽蟹譃椴煌仃噳K，每個(gè)塊的非零元素存儲(chǔ)為0.5 倍二級(jí)緩存尺寸.

3）序列化向量生成（serial vector generation）.根據(jù)矩陣塊對(duì)向量x各分量的索引前后順序，構(gòu)建該矩陣塊對(duì)應(yīng)的重排后的向量x′.

本訪存優(yōu)化方法的研究思路是設(shè)計(jì)預(yù)處理算法，將向量x的不規(guī)則訪問模式轉(zhuǎn)化為對(duì)向量x′的順序訪問模式，借助處理器緩存的硬件步長(zhǎng)預(yù)取器實(shí)現(xiàn)對(duì)重排后向量x′訪問模式的識(shí)別，準(zhǔn)確預(yù)測(cè)程序要使用的向量x′數(shù)據(jù)的內(nèi)存地址并做數(shù)據(jù)預(yù)取，達(dá)到優(yōu)化SpMV 訪存瓶頸的效果.

3 分支預(yù)測(cè)優(yōu)化方法

如第2 節(jié)所述，分支預(yù)測(cè)誤判開銷是影響SpMV性能的重要來源之一，但現(xiàn)有研究尚未對(duì)SpMV 運(yùn)行過程中分支跳轉(zhuǎn)行為進(jìn)行細(xì)粒度分析和優(yōu)化，本文考慮稀疏矩陣的分支跳轉(zhuǎn)特征，提出行重排策略優(yōu)化SpMV 的分支預(yù)測(cè)開銷.

通過圖2 中的SpMV 指令流分析可知，對(duì)于稀疏矩陣每行非零元素的循環(huán)遍歷構(gòu)成了SpMV 的關(guān)鍵分支跳轉(zhuǎn).由于每行非零元素個(gè)數(shù)均不相同，這一分支跳轉(zhuǎn)的行為往往難以被正確預(yù)測(cè).圖7（a）中展示了典型稀疏矩陣web-Stanford的非零元素分布特征.可以看到，行長(zhǎng)度的分布呈現(xiàn)出冪律特征，即絕大部分的矩陣行的非零元素個(gè)數(shù)很少，僅有少量的行中非零元素?cái)?shù)量較多.這一特征被稱為Scale-Free，它是在圖結(jié)構(gòu)中廣泛存在的分布特征[22].圖7（a）中展示了基于BHR（branch history register）和PHT（pattern history table）技術(shù)的分支預(yù)測(cè)情況示例.可以看到，在行長(zhǎng)度隨機(jī)分布的情況下，現(xiàn)有分支預(yù)測(cè)技術(shù)難以進(jìn)行準(zhǔn)確預(yù)測(cè)，并且傾向于在每行循環(huán)遍歷的開始和最后迭代時(shí)產(chǎn)生預(yù)測(cè)錯(cuò)誤.我們還可以看到，在長(zhǎng)行的循環(huán)遍歷中，分支預(yù)測(cè)錯(cuò)誤的比例相對(duì)較低.然而，由于大部分基于圖結(jié)構(gòu)的稀疏矩陣符合Scale-Free 冪律分布特征，大部分行的非零元素?cái)?shù)量都很少，因此導(dǎo)致較高的分支預(yù)測(cè)錯(cuò)誤概率.

為了解決這一問題，我們采用矩陣行重排的方法.如圖7（b）所示，通過按照行長(zhǎng)度進(jìn)行排序，使得非零元素個(gè)數(shù)相同的行構(gòu)成一組（group），并且在內(nèi)存中相鄰排列，保證連續(xù)若干行的分支跳轉(zhuǎn)模式相同，從而提升分支預(yù)測(cè)的準(zhǔn)確率.這一方法可以廣泛適用于基于局部跳轉(zhuǎn)歷史的分支預(yù)測(cè)技術(shù)，對(duì)于傳統(tǒng)的BHR/PHT 分支預(yù)測(cè)器和更為先進(jìn)的TAGE 分支預(yù)測(cè)器[23]都可以達(dá)到很好的預(yù)測(cè)效果.此外，行重排的過程打亂了原有的矩陣分布，因此可能會(huì)破壞某些稀疏矩陣的原有結(jié)構(gòu)特征，從而造成計(jì)算效率的下降.為了抑制這一現(xiàn)象，我們將矩陣中連續(xù)的2 048行劃分為一個(gè)簇（bundle），并且僅在簇內(nèi)開展行重排，從而限制了亂序范圍.實(shí)驗(yàn)結(jié)果表明，采用2 048 的排序范圍已經(jīng)可以實(shí)現(xiàn)很好的分支預(yù)測(cè)優(yōu)化效果和接近100%的分支預(yù)測(cè)正確率.

4 SIMD 計(jì)算優(yōu)化方法

已有方案雖然能實(shí)現(xiàn)較高的SIMD 單元利用率，但是會(huì)引入額外的指令復(fù)雜度，不利于SpMV 計(jì)算性能的優(yōu)化.本文借鑒已有研究思路，提出簡(jiǎn)潔高效的SIMD 并行計(jì)算策略，提升SpMV 的計(jì)算性能.

為了進(jìn)一步提升SpMV 的計(jì)算性能，pvSpMV 采用SIMD 指令集實(shí)現(xiàn)了高效的并行計(jì)算.由式（1）可知，SIMD 計(jì)算效率優(yōu)化的關(guān)鍵在于提升有效計(jì)算密度，即提升非零元素乘加運(yùn)算在整個(gè)執(zhí)行指令流中的密度.為達(dá)成這一優(yōu)化目標(biāo)，需要提升SIMD 向量利用率，并且降低總體的指令發(fā)射數(shù)量.通過分支預(yù)測(cè)優(yōu)化技術(shù)，pvSpMV 已經(jīng)大幅度減少由于分支預(yù)測(cè)錯(cuò)誤而產(chǎn)生的錯(cuò)誤指令發(fā)射.在本節(jié)中，我們通過構(gòu)建簡(jiǎn)潔的計(jì)算核心進(jìn)一步降低程序復(fù)雜度，從而在保持高利用率基礎(chǔ)上減少與有效計(jì)算無關(guān)的指令數(shù)量.

具體來說，分支預(yù)測(cè)優(yōu)化后的矩陣格式呈現(xiàn)多簇的結(jié)構(gòu)，每簇中包含有多組，每組中包含有多個(gè)長(zhǎng)度相同的行.圖8 給出了一個(gè)簇內(nèi)部的結(jié)構(gòu)示例.假設(shè)SIMD 的向量長(zhǎng)度為VL，則可以在每個(gè)組中收集VL整數(shù)倍的行作為一個(gè)段（segment），并且將段中的多個(gè)行按照VL的長(zhǎng)度進(jìn)行轉(zhuǎn)置存儲(chǔ).每組中剩余的行統(tǒng)一存儲(chǔ)為碎片集合（fragment）.最終，一個(gè)簇會(huì)被轉(zhuǎn)化為若干個(gè)段和一個(gè)碎片集合.

Fig.8 Optimization of SIMD computation using two simple computing kernels圖8 基于2 種簡(jiǎn)化運(yùn)算核心的SIMD 計(jì)算優(yōu)化（VL=4）

基于這種存儲(chǔ)格式，pvSpMV 設(shè)計(jì)了2 類計(jì)算核心，分別用于對(duì)段和碎片集合進(jìn)行計(jì)算.如圖8 所示，對(duì)于段的計(jì)算采用簡(jiǎn)單的SIMD 乘累加方式，并行開展VL個(gè)行的計(jì)算.由于段的長(zhǎng)度被規(guī)定為VL的整數(shù)倍，因此可以充分使用SIMD 向量長(zhǎng)度而無需使用遮蔽向量.對(duì)于碎片集合，采用逐行的方式進(jìn)行計(jì)算.在每行的計(jì)算中，首先使用SIMD 指令計(jì)算VL整數(shù)倍的非零元素，然后使用標(biāo)量指令計(jì)算剩余非零元素.這樣就確保了每個(gè)SIMD 指令都可以充分利用向量資源，且無需使用遮蔽向量，降低了代碼的復(fù)雜度.而且，由于標(biāo)量計(jì)算的元素?cái)?shù)量很少（每行的標(biāo)量計(jì)算元素?cái)?shù)量不超過VL-1），因此標(biāo)量計(jì)算并不會(huì)對(duì)整體性能構(gòu)成影響.實(shí)驗(yàn)結(jié)果表明，在大多數(shù)稀疏矩陣中，使用標(biāo)量計(jì)算方法可以實(shí)現(xiàn)絕大多數(shù)元素的全向量長(zhǎng)度的SIMD 并行計(jì)算，僅有不足3%的元素需要使用標(biāo)量指令計(jì)算，因此充分發(fā)揮了SIMD 架構(gòu)的計(jì)算并行度.更重要的是，這種計(jì)算核心的構(gòu)造非常簡(jiǎn)單，代碼復(fù)雜度很低.采用這一優(yōu)化方法，可以大幅度提升有效計(jì)算密度，充分發(fā)揮SIMD 架構(gòu)的計(jì)算能力.

5 pvSpMV 整體框架與實(shí)現(xiàn)

5.1 整體框架

基于第2～4 節(jié)所述的面向訪存、分支預(yù)測(cè)和并行計(jì)算效率的優(yōu)化方法，我們提出pvSpMV 方法，通過構(gòu)建一種新型的稀疏矩陣存儲(chǔ)格式和SpMV 計(jì)算方法，實(shí)現(xiàn)計(jì)算性能的全面提升.

由于很多應(yīng)用場(chǎng)景中需要多次迭代計(jì)算SpMV來進(jìn)行問題求解，因此稀疏矩陣預(yù)處理的開銷可以被有效地分?jǐn)偟蕉啻蔚^程中，其對(duì)整體性能的影響可以忽略.pvSpMV 結(jié)合了上述多種優(yōu)化技術(shù)，通過一系列預(yù)處理方法進(jìn)行基于CSR 格式的稀疏矩陣存儲(chǔ)方法優(yōu)化，提升其在計(jì)算過程中的訪存和分支預(yù)測(cè)的可預(yù)測(cè)性.圖9 給出了pvSpMV 對(duì)于稀疏矩陣預(yù)處理的流程圖，具體劃分為6 個(gè)步驟.

Fig.9 Block diagram of preprocessing flow,computation method and multi-thread implementation of pvSpMV圖9 pvSpMV 的預(yù)處理流程、計(jì)算方法和多線程實(shí)現(xiàn)框圖

1）比特圖優(yōu)化.對(duì)于原始CSR 格式稀疏矩陣，采用基于比特圖的行重排方法進(jìn)行矩陣變換，提升相鄰行的非零元素分布特征，記錄行重排映射關(guān)系R1.

2）塊劃分.首先根據(jù)緩存尺寸確定每塊（block）子矩陣的向量長(zhǎng)度；然后對(duì)比特圖優(yōu)化后的矩陣進(jìn)行分塊，每塊子矩陣包含若干個(gè)連續(xù)的矩陣行，并且每塊子矩陣所索引的向量元素可以在緩存中全部存儲(chǔ).

3）簇劃分.在每塊子矩陣內(nèi)部進(jìn)行簇劃分，每個(gè)簇包含2 048 個(gè)相鄰的行.在簇的內(nèi)部，對(duì)2 048 個(gè)行進(jìn)行重排分組，每組的各個(gè)行具有相同的非零元素?cái)?shù)量，記錄行重排映射關(guān)系R2.

4）生成段和碎片集合.在每個(gè)簇的內(nèi)部，按照SIMD 向量長(zhǎng)度VL，在每組中收集數(shù)量為VL整數(shù)倍的行，并且轉(zhuǎn)置存儲(chǔ)為段；剩余的行組成碎片集合；記錄行重排映射關(guān)系R3.

5）串行化.對(duì)于每塊子矩陣，按照其中各個(gè)簇的計(jì)算和訪存模式，記錄輸入向量的串行化映射關(guān)系S1，并且修改子矩陣數(shù)據(jù)結(jié)構(gòu)中對(duì)應(yīng)位置colidx的值.

6）生成寫回順序.為了支持迭代的SpMV 運(yùn)算，每次迭代后的計(jì)算結(jié)果需要被重組為下次迭代所需的輸入向量順序，pvSpMV 記錄了每行的計(jì)算結(jié)果與下次迭代輸入向量中對(duì)應(yīng)位置的映射關(guān)系，記錄為W1[i]=S1[R3[R2[R1[i]]]]；并且還需要在最后一次迭代后，將運(yùn)算結(jié)果還原為原始向量順序，因此還需要還原的映射關(guān)系，記錄為W2=Rev（R3[R2[R1[i]]]），其中Rev（）表示將映射關(guān)系的鍵和值進(jìn)行對(duì)調(diào)的函數(shù).

綜上所述，經(jīng)過pvSpMV 預(yù)處理后的稀疏矩陣可表示為若干塊包含一系列簇的子矩陣組合，且每個(gè)簇的行均按照SIMD 計(jì)算模式進(jìn)行重組.此外，還會(huì)存儲(chǔ)3 個(gè)映射關(guān)系，分別為S1，W1和W2.

5.2 計(jì)算方法

根據(jù)5.1 節(jié)所述的稀疏矩陣格式，圖9 給出了基于pvSpMV 的多次SpMV 迭代計(jì)算流程.這一流程包含5 個(gè)步驟：

1）在第一次迭代中，輸入向量x0與原始的CSR格式稀疏矩陣進(jìn)行SpMV 運(yùn)算并生成結(jié)果向量y0；對(duì)于y0，使用S1映射對(duì)其進(jìn)行變換，將其轉(zhuǎn)換為符合串行訪問方式的元素順序.

2）使用轉(zhuǎn)換后的向量與各塊子矩陣進(jìn)行SpMV運(yùn)算，運(yùn)算結(jié)果寫入稠密排列的結(jié)果向量中.在這一過程中，每塊矩陣均構(gòu)成了對(duì)于向量元素的串行訪問，因此可以有效觸發(fā)數(shù)據(jù)預(yù)取，并且提升了緩存的局部性.

3）在各塊子矩陣均運(yùn)算結(jié)束后，對(duì)于運(yùn)算生成的結(jié)果向量按照W1進(jìn)行重排，使其符合串行化訪存優(yōu)化的向量元素順序，并且將轉(zhuǎn)換后的結(jié)果向量作為運(yùn)算向量，與優(yōu)化后的稀疏矩陣格式開始下次迭代.

4）重復(fù)步驟2 和步驟3，直到算法收斂.

5）在最后一次SpMV 計(jì)算迭代結(jié)束后，使用W2對(duì)結(jié)果向量進(jìn)行變換，使其元素順序恢復(fù)為原始的向量順序，至此多次迭代的SpMV 完成.

在這5 個(gè)步驟中，核心操作為步驟2 和步驟3，它們會(huì)經(jīng)歷多次迭代，是決定算法性能的關(guān)鍵計(jì)算步驟.其中，步驟2 的向量訪存模式在優(yōu)化后呈現(xiàn)出串行訪問特征，可以有效借助數(shù)據(jù)預(yù)取器實(shí)現(xiàn)高效地取數(shù)，此外，由于分支預(yù)測(cè)優(yōu)化和向量并行計(jì)算優(yōu)化，步驟2 的計(jì)算行為呈現(xiàn)出高度可預(yù)測(cè)性和極低的程序復(fù)雜度特征，可以充分發(fā)揮現(xiàn)代高性能處理器的計(jì)算能力，構(gòu)成了pvSpMV 的主要性能收益來源.與之相對(duì)應(yīng)的，步驟3 需要對(duì)下一次迭代所需的向量進(jìn)行賦值，這一過程需要對(duì)步驟2 生成的結(jié)果數(shù)據(jù)進(jìn)行亂序讀取，因此會(huì)帶來額外的性能開銷.因此，步驟2 的收益與步驟3 的開銷之間的權(quán)衡是很關(guān)鍵的.具體來說，步驟3 的開銷正比于迭代向量的大小，因此與稀疏矩陣的分塊數(shù)量直接相關(guān).本文經(jīng)過理論分析和實(shí)際測(cè)試，發(fā)現(xiàn)按照50%二級(jí)緩存的容量設(shè)置每塊的向量元素?cái)?shù)量是較為合理的.這一尺寸設(shè)計(jì)不僅可以滿足訪存優(yōu)化方法所需的緩存尺寸限制，而且可以避免過多分塊所導(dǎo)致的步驟3 開銷過大問題.實(shí)驗(yàn)表明，pvSpMV 算法中，對(duì)于步驟2 的訪存優(yōu)化所帶來的收益往往會(huì)超過步驟3 亂序?qū)懟厮鶎?dǎo)致的額外開銷，因此整體預(yù)期具有良好的性能優(yōu)化效果.

5.3 多線程計(jì)算方法

pvSpMV 具有較好的多核可擴(kuò)展性，可以通過在多線程中動(dòng)態(tài)分配塊和簇的方法，實(shí)現(xiàn)便捷的多線程計(jì)算擴(kuò)展.圖9 詳細(xì)展示了pvSpMV 的多線程計(jì)算方法.具體來說，在步驟2 中，可以根據(jù)矩陣大小和線程數(shù)量選擇不同的調(diào)度策略.為了確保多線程之間的動(dòng)態(tài)負(fù)載均衡，pvSpMV 在塊數(shù)達(dá)到線程數(shù)2 倍以上的情況下，會(huì)選擇以塊為粒度進(jìn)行多核的動(dòng)態(tài)調(diào)度；當(dāng)塊數(shù)較少時(shí)，pvSpMV 會(huì)采用以簇為粒度的方式進(jìn)行多核計(jì)算調(diào)度，從而保證多核的動(dòng)態(tài)負(fù)載均衡策略不會(huì)失效.在步驟3 中，pvSpMV將需要賦值的下次迭代所需向量按照1 024 個(gè)元素為粒度進(jìn)行了劃分，每個(gè)線程搶占式地進(jìn)行不同組1 024 個(gè)元素的賦值操作，并且在賦值過程中讀取共享的結(jié)果向量數(shù)據(jù).這樣的處理方法既可以確保線程間的負(fù)載均衡，又可以避免由于寫沖突所導(dǎo)致的“偽共享”（false sharing）問題.

6 性能評(píng)估與討論

6.1 實(shí)驗(yàn)設(shè)置

我們?cè)谏逃脁86 平臺(tái)評(píng)估設(shè)計(jì)的性能表現(xiàn)，表1描述了實(shí)驗(yàn)平臺(tái)微架構(gòu)參數(shù)，實(shí)驗(yàn)平臺(tái)為Intel Xeon Gold 6164 服務(wù)器，處理器架構(gòu)為Intel Skylake-X 多核架構(gòu)，包含一級(jí)數(shù)據(jù)緩存和二級(jí)緩存為各個(gè)處理器核心私有，容量分別為32 KB 和1 MB，三級(jí)緩存為多個(gè)處理器核心共享，三級(jí)緩存容量為24.75 MB.為了避免Intel 處理器在高頻工作下觸發(fā)DVFS 的自動(dòng)降頻技術(shù)，本文將處理器頻率固定為2.5 GHz，從而提升測(cè)試環(huán)境的穩(wěn)定性.另外，我們還在AMD 平臺(tái)測(cè)試方案的性能.

Table 1 Evaluation Platform Parameters表1 實(shí)驗(yàn)平臺(tái)微架構(gòu)參數(shù)

用于研究SpMV 性能的基準(zhǔn)方案選擇Intel Math Kernel Library[6]（MKL 版本為2020.2）中提供的SpMV計(jì)算核心實(shí)現(xiàn)，采用無優(yōu)化性能作為基準(zhǔn).pvSpMV代碼采用GCC 編譯器（GCC 11.3）進(jìn)行編譯，并且采用了Intel AVX-512 向量指令集.此外，本文還對(duì)比了其他先進(jìn)SpMV 技術(shù)的性能，包括CSR5[8]，CVR[9]，CSR2[24]，以及經(jīng)過優(yōu)化流程后的MKL SpMV 方法.

本文的測(cè)試數(shù)據(jù)選自公開測(cè)試集Suitesparse Matrix Collection[15]，并選取了其中與圖計(jì)算和科學(xué)計(jì)算相關(guān)的SNAP,Kamvar,Pajek,DIMACS10,Belcastro數(shù)據(jù)集.本文的全部測(cè)試矩陣集合包含84 個(gè)稀疏矩陣，覆蓋了社交網(wǎng)絡(luò)、路徑規(guī)劃、電網(wǎng)分析等多個(gè)應(yīng)用場(chǎng)景.

6.2 關(guān)鍵指標(biāo)評(píng)估

為了分析pvSpMV 對(duì)于分支預(yù)測(cè)和并行計(jì)算效率等性能瓶頸的優(yōu)化效果，我們?cè)诒竟?jié)中對(duì)pvSpMV的關(guān)鍵性能指標(biāo)進(jìn)行評(píng)估.為了測(cè)試各項(xiàng)指標(biāo)，在所有測(cè)試集中挑選了不同類型、不同尺寸的24 個(gè)矩陣，并且分別測(cè)試對(duì)比了每個(gè)矩陣在MKL[6],pvSpMV,CSR5[8],CVR[9]這4 種優(yōu)化方法上的指標(biāo)結(jié)果.

我們首先評(píng)估了分支預(yù)測(cè)優(yōu)化的效果，實(shí)驗(yàn)結(jié)果如圖10 所示.可以看出，MKL 和CSR5 均有較多的分支預(yù)測(cè)錯(cuò)誤率，CVR 由于計(jì)算模式相對(duì)簡(jiǎn)單，分支預(yù)測(cè)的難度更低.而pvSpMV則實(shí)現(xiàn)了分支預(yù)測(cè)的完全優(yōu)化，基本消除了分支預(yù)測(cè)錯(cuò)誤開銷.

Fig.10 Optimization evaluation of branch prediction圖10 分支預(yù)測(cè)優(yōu)化評(píng)估

之后，我們?cè)u(píng)估了SIMD 并行計(jì)算的優(yōu)化效果，具體分為2 項(xiàng)指標(biāo)，分別是：SIMD 指令的向量利用率，實(shí)驗(yàn)結(jié)果見圖11（a）；SIMD 指令的有效計(jì)算密度，實(shí)驗(yàn)結(jié)果見圖11（b）.通過這2 項(xiàng)指標(biāo)的對(duì)比可以看到，CSR5 和CVR 算法雖然可以實(shí)現(xiàn)較高的向量利用率，但是其整體的有效計(jì)算密度提升較為有限，這是由于它們計(jì)算程序的復(fù)雜度導(dǎo)致了額外指令數(shù).與之相比，pvSpMV 方法不僅可以達(dá)到很高的向量利用率，而且通過低復(fù)雜度的程序?qū)崿F(xiàn)，有效抑制了冗余代碼的數(shù)量，確保SIMD 指令的有效計(jì)算密度可以有效提升.

Fig.11 Evaluation of vector usage and computing density圖11 向量利用率和有效計(jì)算密度評(píng)估

最后，我們?cè)u(píng)估了pvSpMV 對(duì)訪存效率的優(yōu)化效果，實(shí)驗(yàn)結(jié)果如圖12 所示.二級(jí)緩存缺失數(shù)越高，則訪存效率越差.可以看到，相較不加優(yōu)化方法的基準(zhǔn)方案，pvSpMV 顯著降低二級(jí)緩存缺失數(shù)，對(duì)訪存效率有明顯提升，這也證明基于Bitmap 的行重排策略和向量x′訪問序列化對(duì)訪存效率的優(yōu)化效果.

Fig.12 Evaluation of L2 cache miss number圖12 二級(jí)緩存缺失數(shù)評(píng)估

6.3 整體性能評(píng)估

本節(jié)對(duì)pvSpMV 的整體性能進(jìn)行了全面的評(píng)估和測(cè)試.基準(zhǔn)性能為MKL 在不進(jìn)行優(yōu)化處理時(shí)的性能.圖13（a）中給出了pvSpMV 以及其他對(duì)比方法在全部85 個(gè)稀疏矩陣上的單線程計(jì)算性能和相比于基準(zhǔn)性能的加速比.可以看到，pvSpMV 在幾乎所有測(cè)試矩陣中均取得了最高的性能，并且相比于基準(zhǔn)方法的收益非?？捎^.具體來說，pvSpMV 的單線程平均算力達(dá)到了1.58 GFLOPS，相比于基準(zhǔn)性能的0.64 GFLOPS 取得了2.6 倍的性能提升，并且相比于CSR5 和CVR 分別取得了1.3 和1.58 倍的加速比.同時(shí)可以觀察到，隨著矩陣尺寸的增加，計(jì)算性能逐漸下降，這是由于訪存瓶頸的逐漸凸顯導(dǎo)致的.通過pvSpMV 的優(yōu)化，可以在大型矩陣的計(jì)算中取得超過1.50GFLOPS 的性能，這說明了本文采取的優(yōu)化方法可以有效緩解訪存瓶頸，提升整體計(jì)算性能.

Fig.13 Performance comparison between pvSpMV and other comparative methods over 85 sparse matrices圖13 pvSpMV 和其他對(duì)比方法在85 個(gè)稀疏矩陣上的性能對(duì)比

我們還評(píng)估了多線程的計(jì)算性能.圖13（b）給出了使用8 線程計(jì)算的性能測(cè)試結(jié)果.為了保證測(cè)試的公平性，我們關(guān)閉了Intel 處理器的SMT 技術(shù)，并且通過計(jì)算核心綁定的方式確保所有線程均運(yùn)行在同一個(gè)NUMA 節(jié)點(diǎn)上.通過實(shí)驗(yàn)結(jié)果可以看到，CVR和MKL 的多核可擴(kuò)展性較好，而pvSpMV 的性能提升則不足8 倍.這一現(xiàn)象的原因有2 點(diǎn)：1）pvSpMV通過提升訪存的可預(yù)測(cè)性來充分發(fā)揮數(shù)據(jù)預(yù)取器的取數(shù)能力，然而多線程模式下可能會(huì)以簇為粒度進(jìn)行多線程調(diào)度，這可能會(huì)破壞向量元素的訪存串行性，從而降低數(shù)據(jù)預(yù)取的收益；2）由于pvSpMV 大幅提升了計(jì)算效率和取數(shù)效率，導(dǎo)致在多線程下的訪存帶寬趨向飽和，因此導(dǎo)致整體性能收益不足8 倍.然而，即使在這種情況下，pvSpMV 仍然取得了最高的多線程性能，并且相比于基線性能實(shí)現(xiàn)了平均2.8倍的加速比，相比CSR5 和CVR 分別實(shí)現(xiàn)了1.3 倍和1.7 倍的加速比.

我們還在AMD 處理器平臺(tái)評(píng)估計(jì)算性能，采用256b 的AVX2 向量指令集，與MKL[6]和CSR2[24]方法比較性能，實(shí)驗(yàn)結(jié)果如圖14 所示.相較基準(zhǔn)方案，pvSpMV 實(shí)現(xiàn)了1.7 倍的性能加速比，證明pvSpMV能部署在多個(gè)不同的處理器平臺(tái)且能實(shí)現(xiàn)顯著的性能提升.

Fig.14 Performance comparison between pvSpMV and other comparative methods on AMD platform圖14 pvSpMV 和其它對(duì)比方法在AMD 處理器平臺(tái)的性能對(duì)比

為了進(jìn)一步分析多線程的性能問題，我們進(jìn)而挑選了典型稀疏矩陣coPapersCiteseer（包含32×106個(gè)非零元素），并且測(cè)試了從線程1 到線程12 不同情況下的性能結(jié)果，如圖15 所示.可以看到，隨著線程數(shù)趨近并超過8 線程，pvSpMV，CSR5 和CSR 均出現(xiàn)了性能增長(zhǎng)變緩，考慮到pvSpMV 有效地優(yōu)化了訪存效率和CPU 執(zhí)行效率，這一現(xiàn)象顯示出訪存帶寬可能已經(jīng)接近飽和狀態(tài).

Fig.15 Performance comparison of different number of threads圖15 不同線程數(shù)的性能對(duì)比

6.4 有效性分析

本文通過結(jié)合多種優(yōu)化手段實(shí)現(xiàn)了pvSpMV 的全面優(yōu)化.本節(jié)將通過一系列消融實(shí)驗(yàn)，逐步探討不同的優(yōu)化方法對(duì)于整體性能的效果.具體來說，為了分項(xiàng)測(cè)試各個(gè)優(yōu)化方案，我們?cè)O(shè)計(jì)了4 個(gè)對(duì)比實(shí)驗(yàn)方案：

1）pvSpMV-base.作為研究pvSpMV 的基線性能方案，這一方案采用最樸素的CSR 計(jì)算方法，采用標(biāo)量指令集，不采用訪存優(yōu)化、分支預(yù)測(cè)優(yōu)化或SIMD指令.

2）pvSpMV-v1.基于pvSpMV-base 版本，增加了基于bitmap 和分塊的訪存優(yōu)化，增加了矩陣訪存模式的串行化處理.

3）pvSpMV-v2.基于pvSpMV-v1 版本，增加了基于分簇和行重排的分支預(yù)測(cè)優(yōu)化方法.

4）pvSpMV-v3.基于pvSpMV-v2 版本，增加了基于分段和碎片集合的SIMD 指令優(yōu)化方法.這一版本是完整的優(yōu)化版本，也是正式的性能評(píng)估版本.

圖16 展示了在全部85 個(gè)測(cè)試矩陣上不同實(shí)驗(yàn)方案相比于基準(zhǔn)性能（pvSpMV-base）的加速比.可以看出，在采用訪存優(yōu)化技術(shù)后，pvSpMV-v1 相比于基準(zhǔn)的樸素CSR SpMV 版本已經(jīng)取得了明顯的收益，并且在使用分支預(yù)測(cè)優(yōu)化和SIMD 優(yōu)化技術(shù)后均可以獲得可觀的性能收益.這一實(shí)驗(yàn)結(jié)果表明，pvSpMV的各項(xiàng)優(yōu)化方法可以有效緩解相關(guān)性能瓶頸，共同完成了計(jì)算性能的全面提升.

Fig.16 Ablation experiments on pvSpMV optimization methods圖16 pvSpMV 各種優(yōu)化方法的消融實(shí)驗(yàn)

6.5 預(yù)處理性能

通常來說，在開始計(jì)算前，SpMV 的優(yōu)化方法需要對(duì)矩陣數(shù)據(jù)格式進(jìn)行一次性地預(yù)處理，通常需要將原始數(shù)據(jù)格式矩陣轉(zhuǎn)為特定的數(shù)據(jù)格式矩陣.預(yù)處理后的矩陣可以在后續(xù)多次的SpMV 迭代中重復(fù)使用.考慮到實(shí)際應(yīng)用場(chǎng)景中SpMV 的迭代次數(shù)可以達(dá)到幾十甚至上百次，因此通常認(rèn)為預(yù)處理的開銷可以被有效分?jǐn)?pvSpMV 遵循了類似的設(shè)計(jì)思路，通過引入預(yù)處理環(huán)節(jié)來完成矩陣的優(yōu)化.為了盡量降低預(yù)處理開銷，我們?cè)陬A(yù)處理環(huán)節(jié)采取極致優(yōu)化性能的策略，基于C 語言實(shí)現(xiàn)了高效的預(yù)處理流程.圖17 給出了各個(gè)優(yōu)化方法的單線程預(yù)處理的開銷，預(yù)處理開銷用單次基準(zhǔn)方案的運(yùn)行時(shí)間量化計(jì)算.可以看到CVR 的預(yù)處理開銷最小，平均需要4.9 倍的基準(zhǔn)方案運(yùn)行時(shí)間.MKL 的預(yù)處理開銷最大，平均需要21.7 倍的基準(zhǔn)方案運(yùn)行時(shí)間.pvSpMV 的預(yù)處理開銷略高于CVR，平均需要11 倍的基準(zhǔn)方案運(yùn)行時(shí)間.真實(shí)場(chǎng)景中SpMV 通常會(huì)迭代多次，并考慮到pvSpMV所提供的潛在性能收益，我們認(rèn)為這一預(yù)處理開銷是可以接受的.

Fig.17 Evaluation of preprocessing costs of different methods圖17 不同方法的預(yù)處理開銷評(píng)估

7 結(jié)論

SpMV 在數(shù)值計(jì)算、大數(shù)據(jù)分析、圖計(jì)算領(lǐng)域和智能計(jì)算中均是非常重要的計(jì)算核心，并且構(gòu)成了現(xiàn)代高性能處理器的主要性能瓶頸之一.本文通過對(duì)SpMV 在現(xiàn)代超標(biāo)量亂序處理器上的指令流組成和運(yùn)行情況的分析，總結(jié)出3 個(gè)關(guān)鍵的性能因素：訪存延遲、分支預(yù)測(cè)錯(cuò)誤開銷和有效計(jì)算密度.基于這一分析結(jié)果本文提出，通過對(duì)稀疏矩陣的存儲(chǔ)格式和計(jì)算方法進(jìn)行修改，可以構(gòu)建一種具有高度預(yù)測(cè)性且程序復(fù)雜度低的計(jì)算核心，充分發(fā)揮現(xiàn)代處理器數(shù)據(jù)預(yù)取、分支預(yù)測(cè)和SIMD 指令集的執(zhí)行效率，實(shí)現(xiàn)訪存瓶頸、分支錯(cuò)誤和計(jì)算效率的全面優(yōu)化.基于這一思路，本文提出了pvSpMV 這一新型的SpMV計(jì)算方法，采用了基于訪問串行化的訪存優(yōu)化技術(shù)、基于行重排的分支預(yù)測(cè)優(yōu)化技術(shù)以及具備高有效計(jì)算密度的SIMD 優(yōu)化技術(shù)，從而實(shí)現(xiàn)了SpMV 的顯著性能提升.測(cè)試結(jié)果顯示，pvSpMV 可以在主流Intel高性能處理器上獲得與主流商用計(jì)算庫MKL 相比平均2.6 倍的加速比，相比于現(xiàn)有最先進(jìn)算法獲得平均1.3 倍的加速比.此外，pvSpMV 所需的預(yù)處理開銷較低，預(yù)期具有很好的實(shí)際應(yīng)用價(jià)值.

作者貢獻(xiàn)聲明：夏天與付格林為共同第一作者，其中夏天提出論文思路，設(shè)計(jì)方法和撰寫論文；付格林負(fù)責(zé)實(shí)驗(yàn)設(shè)計(jì)、數(shù)據(jù)整理與分析、論文修改；曲劭儒負(fù)責(zé)代碼編寫和實(shí)驗(yàn)調(diào)優(yōu)；羅中沛負(fù)責(zé)實(shí)驗(yàn)平臺(tái)搭建；任鵬舉負(fù)責(zé)論文指導(dǎo)和修改.