適合向量化的稀疏矩陣存儲格式研究

1、分類號UDCTP391學(xué)號14063065密級公開工程碩士學(xué)位論文適合向量化的稀疏矩陣存儲格式研究碩士生姓名 謝佩珍工程領(lǐng)域 計算機(jī)技術(shù)研究方向 高性能計算指導(dǎo)教師 徐涵研究員國防科學(xué)技術(shù)大學(xué)研究生院二。一六年十二月論文書脊適合向量化的稀疏矩陣存儲格式研究國防科學(xué)技術(shù)大學(xué)研究生院Research on Sparse Matrix Storage FormatSuitable for VectorizationCandidate： Peizhen XieAdvisor： Prof. Xu HanA thesisSubmitted in partial fulfillment of the req

2、uirementsfor the professional degree of Master of Engineeringin Computer TechnologyGraduate School of National University of Defense TechnologyChangsha, Hunan, P.R.China(December, 2016)獨創(chuàng)性聲明本人聲明所呈交的學(xué)位論文是我本人在導(dǎo)師指導(dǎo)下進(jìn)行的研究工作及取得 的研究成果.盡我所知，除了文中特別加以標(biāo)注和致謝的地方外，論文中不包含 其他人已經(jīng)發(fā)表和撰寫過的研究成果，也不包含為獲得國防科學(xué)技術(shù)大學(xué)或其它 教育機(jī)構(gòu)的

3、學(xué)位或證書而使用過的材料.與我一同工作的同志對本研究所做的任 何貢獻(xiàn)均已在論文中作了明確的說明并表示謝意.學(xué)位論文題目：學(xué)位論文作者簽名：嘩1爪為日期： 由年d月注日學(xué)位論文版權(quán)使用授權(quán)書本人完全了解國防科學(xué)技術(shù)大學(xué)有關(guān)保留使用學(xué)位論文的規(guī)定.本人授權(quán) 國防科學(xué)技術(shù)大學(xué)可以保留并向國家有關(guān)部門或機(jī)構(gòu)送交論文的復(fù)印件和電子 文檔，允許論文被查閱和借閱；可以將學(xué)位論文的全部或部分內(nèi)容編入有關(guān)數(shù)據(jù) 庫進(jìn)行檢索，可以采用影印、縮或掃描等復(fù)制手段保存、匯編學(xué)位論文，（保密學(xué)位論文在解密后適用本授權(quán)書.）學(xué)位論文題田由總向jr代的普 疏I雕件杵修桃 或葉宿學(xué)位論文作者簽名:作者指導(dǎo)教師簽名:日期：E年止月

4、邱日期：目錄摘要 TOC o 1-5 h z HYPERLINK l bookmark41 o Current Document ABSTRACTi HYPERLINK l bookmark48 o Current Document 第一章緒論11.1研究背景11.2相關(guān)課題及研究現(xiàn)狀1 HYPERLINK l bookmark54 o Current Document 1.3論文組織結(jié)構(gòu)5 HYPERLINK l bookmark60 o Current Document 第二章稀疏矩陣向量乘優(yōu)化技術(shù)概述7 HYPERLINK l bookmark63 o Current Document

5、2.1優(yōu)化SpMV算法的時間復(fù)雜度7 HYPERLINK l bookmark66 o Current Document 2.2存儲優(yōu)化72.3硬件技術(shù)優(yōu)化8 HYPERLINK l bookmark69 o Current Document 2.4并行優(yōu)化92.5自適應(yīng)優(yōu)化102.6本章小結(jié)10 HYPERLINK l bookmark72 o Current Document 第三章 CSR格式的一維向量化存儲格式11 HYPERLINK l bookmark75 o Current Document 3.1動機(jī)11 HYPERLINK l bookmark78 o Current Doc

6、ument CSR(r)向量化存儲格式12CSR(r)存儲格式描述123.2.3 rmax 的確定143.2.4實驗結(jié)果及分析20 HYPERLINK l bookmark86 o Current Document CSR(r,ll,12)向量化存儲格式23 HYPERLINK l bookmark110 o Current Document CSR(r,ll,12)存儲格式描述23 HYPERLINK l bookmark113 o Current Document CSR(r,ll,12)格式存儲算法25vmax 數(shù)值的確定27 HYPERLINK l bookmark128 o Curr

7、ent Document 3.3.4實驗結(jié)果及分析29 HYPERLINK l bookmark101 o Current Document 3.4小結(jié)32 HYPERLINK l bookmark104 o Current Document 第四章CSR(r,l1,l2)格式的二維存儲格式33 HYPERLINK l bookmark107 o Current Document 4.1研究動機(jī)33BCSR(r,ll,12)存儲格式描述33BCSR(r,ll,12)存儲格式算法344.4實驗結(jié)果與分析374.4.1預(yù)處理時間374.4.2測試結(jié)果分析37 HYPERLINK l bookmar

8、k131 o Current Document 4.5小結(jié)39 HYPERLINK l bookmark134 o Current Document 第五章結(jié)束語40 HYPERLINK l bookmark137 o Current Document 5.1工作總結(jié)40 HYPERLINK l bookmark144 o Current Document 5.2工作展望41 HYPERLINK l bookmark149 o Current Document 致謝42 HYPERLINK l bookmark152 o Current Document 參考文獻(xiàn)43 HYPERLINK l

9、bookmark196 o Current Document 作者在學(xué)期間取得的學(xué)術(shù)成果46表 目 錄 TOC o 1-5 h z 表1.1稀疏矩陣A的COO存儲格式2表1.2稀疏矩陣A的CSR存儲格式3表1.3稀疏矩陣A的ELL存儲格式3表1.4稀疏矩陣A的DIA存儲格式4表3.1矩陣A12表3.2稀疏矩陣A的CSR(r)存儲格式13表3.3 rmax取不同值時稀疏矩陣A的CSR(r)存儲格式(rmax=(0,l,2,3,4,5,6) . 15 表3.4用于實驗評估的矩陣信息21表3.5矩陣B24表3.6稀疏矩陣B的CSR(r,ll,I2)存儲格式24表3.7 vmax取不同值時稀疏矩陣B的

10、CSR (r,ll,12)存儲格式(vmax=(0,l,2) 28 表4.1 CSR(r,ll,12)格式下的矩陣B33表4.2稀疏矩陣B的BCSR(r,ll,12)存儲格式34圖目錄 TOC o 1-5 h z 圖2.1計算存儲示意圖8圖2.2SIMD結(jié)構(gòu)示意圖9圖2.3流水線技術(shù)圖9圖3.1 CSR(r)(rmax=0,10)稀疏矩陣存儲格式性能圖19圖3.2CSR(r)存儲格式預(yù)處理時間成本22圖 3.3CSR、 CSR(r)、 MKL 性能對I：匕圖23圖 3.4CSR、 CSR(r)、 MKL 力口速對比圖23圖3.5 CSR(r,ll,12)(vmax=0,2)稀疏矩陣存儲格式性能

11、圖29圖3.6預(yù)處理成本圖30圖 3.7CSR、CSR(r)、CSR(r,ll,12)、MKL，性能對I:匕圖31圖 3.8CSR、CSR(r)、CSR(r,ll,12)、MKL 加速比圖31圖4.1 BCSR(r,ll,12)存儲格式預(yù)處理成本圖37圖 4.2 CSR、CSR(r)、CSR(r,ll,12)、BCSR(r,Il,12)、MKL 性能對比圖38圖 4.3CSR、CSR(r)、CSR(r,ll,12)、BCSR(r,ll,12)、MKL 加速比圖39摘要稀疏矩陣向量乘(SpMV)(y = A * x)廣泛用于科學(xué)計算和工程計算中，如大規(guī)模 線性代數(shù)系統(tǒng)的求解，粒子輸運模擬，流體動

12、力學(xué)偏微分方程的求解，天體物理 學(xué)和偏微分問題等。它被歸類為“七個小矮人”的成員，即被認(rèn)為在下一個十年 最重要的數(shù)值方法之一。因此對于稀疏矩陣向量乘(SpMV)及其優(yōu)化技術(shù)的研究有 助于提升解決相關(guān)領(lǐng)域問題的運算效率，有著巨大的研究價值與意義。由于稀疏矩陣含有大量零元素，這直接導(dǎo)致SpMV訪存的不規(guī)則性和差的浮 點性能oSpMV訪存的不規(guī)則性不僅使得計算平臺在進(jìn)行稀疏矩陣向量乘運算時很 難充分使用向量單元進(jìn)行加速，并且還會增加Cache未命中次數(shù)。由于稀疏矩陣 自身的特點，使得稀疏矩陣向量乘運算的實現(xiàn)對稀疏矩陣的存儲格式依賴十分嚴(yán) 重。針對目前主流的稀疏矩陣存儲格式CSR,本文提出了新的適合向

13、量化的稀疏 矩陣存儲格式CSR(r), CSR(r,ll,12), BCSR(r,ll,12)0本文的主要工作總結(jié)如下：查閱并研究國內(nèi)外現(xiàn)有優(yōu)化技術(shù)，從面向計算體系結(jié)構(gòu)的優(yōu)化方面入手， 論述、總結(jié)并歸納了在該方向上現(xiàn)有優(yōu)化技術(shù)的優(yōu)勢與不足，從而為本文的研究 提供了基本的研究方向。對目前主流存儲格式CSR進(jìn)行改進(jìn)，提出了新的適合向量化的稀疏矩陣 存儲格式 CSR(r), CSR(r,ll,12)o 與傳統(tǒng)格式 CSR 相比，CSR(r)、CSR(r,ll,12)的性 能分別提局了 45 % , 49 % o將CSR(r,ll,12)向量化存儲格式與塊存儲格式相結(jié)合，提出二維的稀疏矩陣 存儲格式B

14、CSR(r,ll,12)o與傳統(tǒng)格式CSR相比，BCSR(r,ll,12)的性能提高了 66%。對本文的上述方面研究作了總結(jié)性的概括，給出了本課題今后的研究方 向，展望并提出下一步工作。關(guān)鍵字：稀疏矩陣向量乘；行壓縮稀疏矩陣存儲格式；向量化；性能優(yōu)化ABSTRACTThe Sparse matrix-vector multiplication (SpMV) operation(y=A*x) is widely used in scientific and engineering calculations such as linear algebra, particle transport si

15、mulation, hydrodynamics partial differential equation solution, astrophysics and fractional problems. It has been categorized as a member of the seven dwarfs, i.e. the numerical methods which are believed to be important for at least the next decade. The study of sparse matrix vector multiplication

16、helps to improve the research in there related fields and solve the problem of computing efficiency. Therefore, this study is meaningful.In SpMV, sparse matrices are dominated by a large number of zeros which cause irregular memory access patterns and poor floating-point performance. This irregulari

17、ty results in hardly using the vector unit and increasing the number of cache miss. Achieving high performance of SpMV, which heavily relies on the storage format of sparse matrices in memory, requires to choose a compact data structure and code transformations especially for the calculations involv

18、ing large-scale sparse matrices. In this paper, series optimization methods of SpMV have been approached. In order to take full advantage of the SIMD (single instruction multiple data) acceleration technology in SpMV, we proposed three CSR-extended format called CSR(r), CSR(r,ll,12) and BCSR(r,ll,12

19、).The main work of this paper can be summarized as follows：With the study of the existing optimization techniques, this paper focuses on the aspect of computer architecture, summaries and presents four kinds of optimization strategy.In order to take full advantage of the SIMD (single instruction mul

20、tiple data) acceleration technology in SpMV, we proposed two CSR-extended format called CSR(r) and CSR(r,ll,I2). Compared with the traditional CSR format, the CSR (r) and CSR (r, 11, 12) respectively gains 45% and 49% performance improvement.Combined with the block structure, we proposed another CSR

21、-extended format called BCSR(r,ll,12). The performance of BCSR(r,11,12) increased 66% compared with the CSR format.With the summary, this paper shows the optimization direction of this paper and presents the future work.Keywords： sparse matrix-vector multiplication (SpMV); compressed sparse row (CSR

22、); Single Instruction Multiple Data (SIMD); performance optimization第一章緒論1.1研究背景稀疏矩陣向量乘法(SpMV)操作(y = A*x)廣泛用于科學(xué)計算和工程計算 中，如大規(guī)模線性代數(shù)系統(tǒng)的求解，醫(yī)學(xué)成像，粒子輸運模擬和天體物理學(xué)等。 稀疏矩陣向量乘在許多物理過程離散求解的過程中稀疏矩陣有非常重要的作用， 如在圖論里稀疏矩陣用于實現(xiàn)圖的變換，在解流體力學(xué)的問題上稀疏矩陣用來求 解偏微分方程等等?，F(xiàn)今，稀疏矩陣的研究已經(jīng)滲入到很多領(lǐng)域，如流體動力學(xué)、 圖形學(xué)、粒子輸運、攝影測繪學(xué)等方面研究中，都出現(xiàn)了幾百萬階的大規(guī)模稀疏

23、矩陣。因此對于稀疏矩陣向量乘這一應(yīng)用程序的優(yōu)化和加速的研究，有助于提升 解決相關(guān)領(lǐng)域問題的運算效率，有著十分巨大的科研價值與意義。在過去的幾十年中，稀疏矩陣矢量乘法(SpMV) 一直用于研究最稀疏的基礎(chǔ) 線性代數(shù)子程序庫中的(BLAS)程序。SpMV是指將大小為mXn的稀疏矩陣A 乘以大小為n的密集向量x,并獲得大小為m的密集向量y的一個核心計算程序。 它的實現(xiàn)原理十分簡單，并且可以通過添加簡單的一些編譯器指對其進(jìn)行優(yōu)化在 SpMV中，稀疏矩陣內(nèi)大量的零元素，這直接導(dǎo)致訪存的不規(guī)則性和差的浮點計算 性能。這種不規(guī)則性導(dǎo)致在進(jìn)行稀疏矩陣向量乘運算時很難充分使用向量單元并 增加了 Cache未命中

24、的次數(shù)。由于稀疏矩陣自身的特點，使得稀疏矩陣向量乘運 算的實現(xiàn)對稀疏矩陣的存儲格式依賴十分嚴(yán)重。要實現(xiàn)高性能的SpMV操作，往 往需要對稀疏矩陣選擇十分合適的存儲格式。尤其是對于涉及大規(guī)模稀疏矩陣的 計算，需要對稀疏矩陣存儲格式的選擇至關(guān)重要。因此在大規(guī)模稀疏矩陣向量乘 運算中，進(jìn)一步研究與改進(jìn)稀疏矩陣的存儲格式顯得尤為重要。然而向量單元作為現(xiàn)代處理器上的重要加速組件，目前的優(yōu)化技術(shù)并未很好 地結(jié)合向量化加速技術(shù)，充分發(fā)揮出計算機(jī)向量化加速部件的優(yōu)化潛能。向量單 元將計算機(jī)程序從每次處理單個操作數(shù)的標(biāo)量實現(xiàn)轉(zhuǎn)換為向量實現(xiàn)，同時對多對 操作數(shù)處理一個操作。并且大部分稀疏矩陣都具有良好的數(shù)據(jù)局部性

25、特征，可以 充分利用向量單元對SpMV這一核心程序進(jìn)行加速。因此本研究將以適合運用向 量化技術(shù)并行加速的稀疏矩陣存儲格式為研究對象，改進(jìn)現(xiàn)有存儲格式，提出新 的適合向量化的稀疏矩陣存儲格式。1.2相關(guān)課題及研究現(xiàn)狀稀疏矩陣含有大量的非零元，目前稀疏矩陣的存儲格式大多數(shù)并不采取同時存儲零元和非零元的方式進(jìn)行存儲。目前在實際應(yīng)用中廣泛使用的稀疏矩陣存儲 格式主要有以下四種：(1) Coordinate(COO)：稀疏矩陣的每一個元素都可以用一個三元組來表示， 即元素的數(shù)值及該元素對應(yīng)的行坐標(biāo)和列坐標(biāo)。COO存儲格式就是采取的這種方 式存儲稀疏矩陣。如表1.1所示，COO格式引入有了三個一維數(shù)組，這

26、三個數(shù)組的 名稱依次是row indices,column indices,valuesorow indices數(shù)組存儲的是每個非零元 所在的行坐標(biāo)，column indices存儲的是每個非零元所在的列坐標(biāo)，values每個非 零元的數(shù)值。對于非零元個數(shù)為nnz的稀疏矩陣來說，COO存儲格式的三個一維 數(shù)組(row indices,column indices,values)長度顯然為 nnz。設(shè)一維數(shù)組 row indices 和column indices,每個元素所占的存儲空間為I, values每個元素所占的存儲空間 為D,則存儲該稀疏矩陣總開銷為nnz (2I+D)0表1.1稀疏矩

27、陣A的C00存儲格式1100001111100001011100000011100000011100000011101000011111000011舉例矩陣A：COO表示:Compressed Sparse Row (CSR):CSR存儲格式是在COO存儲格式的基礎(chǔ) 上改進(jìn)的行壓縮存儲格式。如表1.2所示，CSR格式引入有了三個一維數(shù)組，這三 個數(shù)組的名稱依次是row offsets,column indices,valueso row offsets數(shù)組存儲的是每 行的行偏移，行偏移表示某一行的第一個元素在values里面的起始偏移位置。如 圖1.2中，第一行元素1是0偏移，第二行元素2是2

28、偏移，第三行元素5是4偏 移，第4行元素6是7偏移，行偏移的最后一個元素的數(shù)值為矩陣總的非零元素 個數(shù)。column indices存儲的是每個非零元所在的列坐標(biāo)，values每個非零元的數(shù) 值。CSR格式中列號和數(shù)值與COO格式是一致的，唯一的不同在于CSR格式用 存儲行偏移量的存儲方法代替了 COO格式中存儲每個非零元素行坐標(biāo)的存儲方 法。對于非零元個數(shù)為nnz大小為mXn的稀疏矩陣來說，CSR存儲格式的三個一 維數(shù)組(row offsets,column indices,values)長度分別為 m+l,nnz,nnz。設(shè)一維數(shù)組row offsets和column indices,每個

29、元素所需存儲空間為I, values每個元素所需存儲 空間為D,則存儲該稀疏矩陣總開銷為(m+l)*I+nnz*(I+D)。從存儲空間大小上來 說CSR格式要優(yōu)于COO格式。表1.2稀疏矩陣A的CSR存儲格式COO表示：Row Indices 4 11 17 24 1100001111100001011100000011100000011100000011101000011111000011舉例矩陣A：0167012712323434545605670167Column Indices1111111111111111111111111111Values ELLPACK(ELL)2： ELL存儲

30、格式用兩個和原始矩陣相同行數(shù)的矩陣來存 儲稀疏矩陣。如表1.3所示，ELL格式引入了兩個大小相同的二維數(shù)組，這兩個二 維數(shù)組的名稱依次是co山mn indices,values。column indicesc此二維數(shù)組存儲的 是元素的數(shù)值，其每行數(shù)組長度為該稀疏矩陣非零元素最多的那一行的非零元素 的個數(shù)。每行按順序?qū)⒚恳粋€非零元素的列號寫入此二維數(shù)組對應(yīng)的數(shù)組行，若 未填寫完成該行數(shù)組則用*進(jìn)行填寫。values此二維數(shù)組存儲的是每一個非零元的 數(shù)值，其每行數(shù)組長度與column indices二維數(shù)組一致。每行按順序?qū)⒚恳粋€非 零元素的數(shù)值寫入此二維數(shù)組對應(yīng)的數(shù)組行，若未填寫完成該行數(shù)組則用

31、*進(jìn)行填 寫。在E-LL格式中每個元素的行號即為它存儲的二維數(shù)組的行號。對于行非零元 素個數(shù)最多為mrnnz大小為mX n的稀疏矩陣來說，ELL存儲格式的引入的兩個 二維數(shù)組(co山mn indices, values)大小均為 mrnnz*m。假設(shè)二維數(shù)組cdumn indices每個元素所占的存儲空間為, values每個元素所占的存儲空間為D,則存儲該稀疏矩陣總開銷為mrnnz*m* (I+D)。1100001111100001011100000011100000011100舉例矩陣A：01670127123*234*345*Column Indices表1.3稀疏矩陣A的ELL存儲格式

32、00001110111*456*10000111111105671100001111110167Diagnoal(DIA)：DIA存儲格式采取的存儲稀疏矩陣對角線的方法來存儲 稀疏矩陣，此格式適用于對角性特征明顯的稀疏矩陣。如表1.4所示，DIA格式引 入了一個一維數(shù)組和一個二維數(shù)組，數(shù)組的名稱依次是diagnoal indices, values. 數(shù)組diagnoal indices記錄的是對角線編號，對角線編號的規(guī)矩是稀疏矩陣的主對 角線編號為零，往右上角走依次編號為1,2,3.n,往左下角走依次編號為 -1 ,-2,-3-rio 圖中 diagnoal indices 中的數(shù)值依次應(yīng)為

33、-2,0,1 0 diagnoal indices長度為稀疏矩陣中存在非零元素的對角線條數(shù)。values二維數(shù)組存儲的是每條對 角線非零元素的數(shù)值。從主對角線，第一個對角線是零忽略，第二個對角線是5, 6,第三個對角線是零忽略，第四個對角線是1, 2, 3, 4,第五個對角線是7, 8, 9,第六第七個對角線忽略。將每條對角線的非零元填入values二維數(shù)組中采用 的是行對行的規(guī)則，若該列數(shù)組未填寫完成則用*代為填寫。比如圖1.4中5和6 是分別在第三行第四行的，則在前面補上無效元素*。在DIA格式中每個數(shù)值的行 號即是其二維數(shù)組的行號，列號為行號加上對角線編號。對于存在非零元素的對 角線條數(shù)

34、為dn大小為mxn的稀疏矩陣來說，DIA存儲格式的引入的兩個二維數(shù) 組(diagnoal indices, values)大小依次為 dn, dn*m。假設(shè)二維數(shù)組 column indices每個元素所占的存儲空間為I, values每個元素所占的存儲空間為D,則 存儲該稀疏矩陣總開銷為(dn*l+dn*m*D)。表1.4稀疏矩陣A的DIA存儲格式Diagnoal Indices舉例矩陣A：-7-6-1016711000011Values11100001*1111*1111*01110000*111*00111000*111*00011100*111*00001110*111*1000011

35、1*1111*110000111111*針對以上四種稀疏矩陣主要存儲結(jié)構(gòu)，研究人員近年來做了大量研究去尋找 改進(jìn)和突破3-3%佛羅里達(dá)大學(xué)于3】中給出了可以用于檢驗稀疏矩陣性能的測試矩陣。目前對于加速或者優(yōu)化稀疏矩陣向量乘的研究主要有以下兩大類，一是降低 SpMV這一應(yīng)用程序本身的算法復(fù)雜度；二是根據(jù)體系結(jié)構(gòu)的特點來對SpMV操 作進(jìn)行加速。在降低SpMV這一應(yīng)用程序本身的算法復(fù)雜度上，文獻(xiàn)4和文獻(xiàn)5 均在這個方向上取得了突破。楊超，劉芳芳等研究人員提出了一種適合向量化的 稀疏矩陣存儲格式CSRL,該存儲格式挖掘稀疏矩陣非零元之間的數(shù)據(jù)局部性特 征，通過存儲內(nèi)部數(shù)據(jù)連續(xù)的非零元片段，采用SIM

36、D指令直接進(jìn)行存儲和計算， 使得SpMV這一應(yīng)用程序的性能得到了提升。Weifeng Liu, Brian Vinter設(shè)計出了 一種新的稀疏矩陣存儲格式CSR5,此格式在多個計算平臺上(dual-socket Intel CPUs, nVidia GPU, AMD GPU and an Intel Xeon Phi)均取得了十分不錯的實測 性能。Eun Jin Im與Vuduc】等人設(shè)計出了基于高速緩存-寄存器的塊數(shù)據(jù)結(jié)構(gòu)， 并實現(xiàn)了針對不同計算平臺自適應(yīng)調(diào)優(yōu)的算法。繼這項研究之后，張云泉等學(xué)者目 于提出了一種基于塊結(jié)構(gòu)自適應(yīng)調(diào)優(yōu)算法，可以針對不同的矩陣不同的計算平臺 給出最優(yōu)的塊結(jié)構(gòu)大小。

37、文獻(xiàn)10提出了一種數(shù)據(jù)結(jié)構(gòu)，此數(shù)據(jù)結(jié)構(gòu)通過壓縮數(shù)據(jù) 集來降低在進(jìn)行SpMV操作時對內(nèi)存帶寬的需求。文獻(xiàn)11,12提出了 SpMV這一 運算在內(nèi)存中進(jìn)行性能優(yōu)化的方法。文獻(xiàn)13提出了混合用精度的方式來提升 SpMV的性能。Picahel等學(xué)者在文獻(xiàn)14,15,16中充分利用稀疏矩陣非零元素具 有數(shù)據(jù)局部性這一特點，提出了啟發(fā)式重排技術(shù)，提升了 SpMV的性能。孫相征 等學(xué)者中提出了一種新的存儲結(jié)構(gòu)CRSD17, Blelloch等學(xué)者【照設(shè)計出了對角塊 存儲結(jié)構(gòu)HDB,均優(yōu)化了稀疏矩陣向量乘這一核心計算程序的性能。從上述稀疏矩陣向量乘的學(xué)術(shù)成果的總結(jié)中認(rèn)識到，稀疏矩陣向量乘可以從 兩個方面進(jìn)行性

38、能優(yōu)化和加速。一是降低SpMV這一應(yīng)用程序本身的算法復(fù)雜度； 二是根據(jù)體系結(jié)構(gòu)的特點來對SpMV操作進(jìn)行加速。而且學(xué)術(shù)上這兩個方向都有 豐碩的學(xué)術(shù)成果，稀疏矩陣向量乘這一和核心計算程序效果提升顯著。但是，對 于稀疏矩陣向量乘的研究很多需要繼續(xù)去進(jìn)行算優(yōu)化。首先針對降低時間復(fù)雜度， 我們?nèi)钥梢愿M(jìn)一步進(jìn)行算法優(yōu)化。針對體系結(jié)構(gòu)的特點提出的性能優(yōu)化技術(shù)很 難做到具有普適性。即對任何稀疏矩陣任何平臺都具有好的加速效果。比如很多 存儲的格式的提出，都依賴于特定計算平臺，特定的硬件技術(shù)，或者具有特定特 征的系數(shù)矩陣，這樣該項技術(shù)在實際應(yīng)用的稀疏矩陣無法做到都有好的加速效果。 因此對于稀疏矩陣向量乘這一核

39、心計算程序的加速和性能優(yōu)化仍有較大研究的空 間。1.3論文組織結(jié)構(gòu)在SpMV中，稀疏矩陣含有大量零元素，這直接導(dǎo)致訪存的不規(guī)則性和差的 浮點計算性能。這種不規(guī)則性導(dǎo)致在進(jìn)行稀疏矩陣向量乘運算時很難充分使用向 量單元并增加了 Cache未命中的次數(shù)。由于稀疏矩陣自身的特點，使得稀疏矩陣 向量乘運算的實現(xiàn)對稀疏矩陣的存儲格式依賴十分嚴(yán)重。本文針對目前主流的稀 疏矩陣存儲格式CSR,結(jié)合向量化技術(shù)，本文提出了三種新的適合向量化的稀疏 矩陣存儲格式CSR(r), CSR(r,ll,12), BCSR(r,ll,12)。這三種稀疏矩陣向量化格式， 充分挖掘了稀疏矩陣非零元數(shù)據(jù)的局部性特征，采用向化加速技

40、術(shù)，顯著提升了 SpMV這一核心計算程序的性能論文分在結(jié)構(gòu)上分為五章：第一章 緒論。本章介紹了研究課題的背景和意義，闡述了目前主流的的四種 稀疏矩陣存儲格式和目前國內(nèi)外在稀疏矩陣向量乘這一領(lǐng)域的研究現(xiàn)狀。此外本 章總結(jié)本文所做的工作和創(chuàng)新點。第二章稀疏矩陣向量乘優(yōu)化技術(shù)概述。此章著重介紹稀疏矩陣向量乘現(xiàn)有的 優(yōu)化技術(shù)，并確定課題的研究方向。第三章 基于CSR格式的一維向量化存儲格式設(shè)計。本章提出了兩種適合向量 化的稀疏矩陣存儲格式CSR(r)和CSR(r,ll,12),并詳細(xì)介紹了這兩種存儲格式存儲 算法，和實測性能。第四章CSR(r,ll,12)格式的二維存儲格式設(shè)計。本章在CSR(r,ll

41、,12)存儲格式的 基礎(chǔ)上，與塊存儲格式相結(jié)合，提出了 CSR(r,ll,12)格式的二維存儲格式 BCSR(r,ll,12),并詳細(xì)闡述此存儲格式，及其在SpMV上的實測性能。第五章總結(jié)全文，對本文所做的研究工作進(jìn)行總結(jié)，闡述了本文的主要科研 成果和創(chuàng)新點。并提供了在稀疏矩陣向量乘這個領(lǐng)域今后可能繼續(xù)進(jìn)行優(yōu)化的方 向，給出了下一步的工作規(guī)劃。第二章 稀疏矩陣向量乘優(yōu)化技術(shù)概述2.1優(yōu)化SpMV算法的時間復(fù)雜度Strassen19提出一種基于分而治之的快速矩陣乘法，該方法可以對矩陣進(jìn)行劃 分，并將稀疏矩陣盡可能把地將矩陣劃分為四個維度盡量相同的子矩陣塊。然后 運用將浮點乘法運算轉(zhuǎn)換為加法運算的

42、方式，降低子矩陣塊間SpMV算法的時間 復(fù)雜度，提升了 SpMV的性能。另外，Yuster等人在文獻(xiàn)20中設(shè)計出了一種通過 構(gòu)造子矩陣的乘法順序的SMP算法，該算法整體提升了 SpMV性能。但是，此類減低算法時間復(fù)雜度的算法也存在數(shù)值穩(wěn)定性地低這一缺陷。因 此對于某些對精度要求高，此類算法計算的結(jié)果無法滿足其需求。另外，實際應(yīng) 用中需要大量執(zhí)行浮點運算的大規(guī)模稀疏矩陣，由于相對誤差的累積，計算的結(jié) 果很難達(dá)到計算所需要精確度。降低SpMV算法的時間復(fù)雜度的算法可以有效提 升SpMV的性能，節(jié)省計算時間，但是某些計算精度要求很嚴(yán)或者是大量浮點計 算相對誤差累積影響到計算所需精度時，此類算法就無法

43、滿足所需的運算要求。2.2存儲優(yōu)化目前主流的計算平臺上，稀疏矩陣向量乘的運行效率不到機(jī)器浮點運算峰值 的百分之十。其中一個非常重要的原因就是SpMV算法中的浮點計算與存儲訪問 的比率太低，大量的時間花費在數(shù)據(jù)的存儲訪問中.為了提高SpMV的性能，目前 文獻(xiàn)中提出了一系列的優(yōu)化方法，如：寄存器分塊技術(shù)、cache分塊技術(shù)21,22壓 縮存儲技術(shù)等。這些方法的主要目的是減少訪存次數(shù)，以及通過數(shù)據(jù)重組提高 數(shù)據(jù)訪問的局部性，從而提升SpMV性能。由于CPU與存儲器性能差大且存儲層次的復(fù)雜,不同的存儲層次訪問時間差 距很大（圖2.1）0訪問寄存器的所需時間大約為1拍，訪問一級Cache的所需時間 大約

44、為3拍，訪問二級Cache的所需時間大約為10拍，訪問本地內(nèi)存的所需時間 大約為150拍。因此要提高SpMV應(yīng)用程序的性能十分困難，必須減少SpMV的 訪存次數(shù)開發(fā)稀疏矩陣的數(shù)據(jù)局部性，從而減少通信和訪存時間。Proceor-DRAM Memory Gap阡roc60%/yr (2X/1.5yr) Moores LawCPUProcessor-MemoryPerformance Gap:，(grows 50% I year)DRAM岫 9%/yr.(2X/ioyrs).圖2.1計算存儲不意圖寄存器1拍LI Cachei m 拍jL2 Cache 1。拍L3 Cache30拍，本地內(nèi)存15。拍非

45、本地內(nèi)存1500拍2.3硬件技術(shù)優(yōu)化在OSKI【3】與pOSKI【23】引入了有關(guān)底層部分的優(yōu)化操作。硬件技術(shù)優(yōu)化主要包 含兩種優(yōu)化技術(shù)，即向量化優(yōu)化，流水線優(yōu)化。大量現(xiàn)代計算平臺已在處理器中實現(xiàn)了 SIMD體系結(jié)構(gòu)和指令集，短向量 成為提升CPU單核性能的重要方法。SIMD （ Single Instruction Multiple Data,單指 令多數(shù)據(jù)）體系結(jié)構(gòu)（圖2.2）是普遍而高效的并行體系結(jié)構(gòu)，它能在不提升取指令、 指令解碼帶寬的基礎(chǔ)上提高執(zhí)行吞吐率【3%在微處理器體系結(jié)構(gòu)中扮演著非常重 要的角色。Intel Xeon E5-2670向量寄存器寬度達(dá)到256位。Intel Xeo

46、n Phi向量 寄存器的寬度則達(dá)到了 512位。向量化部件出現(xiàn)后，計算機(jī)可以對多個操作數(shù)進(jìn) 行矢量實現(xiàn)，即計算機(jī)可以同時對多個操作數(shù)進(jìn)行操作。Intel Xeon支持AVX指 令集，有兩 個128bit的讀端口，一個時鐘周期內(nèi)可以同時發(fā)射2條讀指令，讀取 4個雙精度浮點數(shù)，支持256位的向量計算.劉芳芳利用AVX指令集優(yōu)化SpM-V 這一核心程序性能有明顯的提升。圖2.2SIMD結(jié)構(gòu)示意圖流水線技術(shù)（圖2.3）通過將計算機(jī)指令分解為多步，并讓各步操作重疊，從而 實現(xiàn)指令級并行優(yōu)化。但是實現(xiàn)Pipeline優(yōu)化的主要依賴于代碼重構(gòu)和編譯優(yōu)化， 提升SpMV性能有限。第一拍第二拍第三拍第四拍第五拍

47、第六拍第七拍第八拍取指令求和11二|跚，取指令求和2 卜|1，|2.4并行優(yōu)化近30年來提高系統(tǒng)性能有三駕馬車，主頻，指令級并行，并行度。近年來前 兩駕馬車無力，提高性能全靠增加并行度。2007年銀河五號，速度為0.1P,具有2048 個結(jié)點，0.4096萬個核；2010年天河一號，速度為4.7P,有7168個結(jié)點，18.6368 萬個核；2013年天河二號，速度為55.6P,有16000個結(jié)點，312.0000萬個核。 由以上數(shù)據(jù)分析可得知，SpMV應(yīng)用程序要充分發(fā)揮多核計算機(jī)的強(qiáng)大的計算能 力,必須具備良好的并行性。pOSKI實現(xiàn)了多線程并行的SpMV優(yōu)化，優(yōu)化效果 顯著。但是，文獻(xiàn)19中

48、亦提出了一種基于多進(jìn)程的MPI優(yōu)化，由于節(jié)點間的通 訊開銷大，優(yōu)化效果平平。2.5自適應(yīng)優(yōu)化自適應(yīng)優(yōu)化技術(shù)是針對不同計算平臺不同稀疏矩陣從多種運算方法中選擇一 種使SpMV性能最優(yōu)的運算方法。文獻(xiàn)9提出一種可以針對不同計算平臺不同稀 疏矩陣找到該稀疏矩陣最適合的固定分塊大小的自適應(yīng)優(yōu)化方法。該方法首先通 過稠密矩陣乘向量的方式計算出不同固定大小的矩陣塊在該本臺上的執(zhí)行時間。 然后針對不同的稀疏矩陣用不同的固定大小的塊進(jìn)行劃分，得出該矩陣在固定塊 劃分的情況下得到的塊的數(shù)目。然后兩者相乘，選擇時間最短即性能最優(yōu)的塊大 小對稀疏矩陣進(jìn)行劃分，進(jìn)行計算。自適應(yīng)優(yōu)化技術(shù)一般情況下只是針對某一特 定的優(yōu)

49、化策略，而不是一種集成的調(diào)優(yōu)方案，因而這類優(yōu)化不具備普適性。2.6本章小結(jié)本章詳細(xì)介紹了五種在稀疏矩陣向量乘這一研究領(lǐng)域比較常用的五種優(yōu)化技 術(shù)。即降低時間復(fù)雜度的優(yōu)化技術(shù)、存儲優(yōu)化技術(shù)、硬件優(yōu)化、并行優(yōu)化和自適 應(yīng)優(yōu)化?？偨Y(jié)了這五類優(yōu)化技術(shù)的優(yōu)化策略，以及不足之處。SpMV在科學(xué)研究和 工程領(lǐng)域都有廣泛運用，因此很多學(xué)者對其進(jìn)行研究，在多個優(yōu)化方向上均有突 破性進(jìn)展。本章認(rèn)為稀疏矩陣向量乘的優(yōu)化可以繼續(xù)在五個方向上繼續(xù)深入，提 出更好的優(yōu)化算法。第三章CSR格式的一維向量化存儲格式3.1動機(jī)稀疏矩陣矢量乘法(SpMV)操作(y=A*x)廣泛用于科學(xué)和工程計算，如 線性代數(shù)，粒子運輸模擬，流體

50、動力學(xué)偏微分方程解，圖變換，天體物理學(xué)和分 數(shù)問題【2,7,24,33。在SpMV中，稀疏矩陣由大量的零來支配，這導(dǎo)致不規(guī)則的存儲 器訪問模式和差的浮點性能。這種不規(guī)則性導(dǎo)致幾乎不使用向量單元和增加高速 緩存未命中的數(shù)量。實現(xiàn)高性能的SpMV,其嚴(yán)重依賴于存儲器中的稀疏矩陣的存 儲格式，需要選擇緊湊的數(shù)據(jù)結(jié)構(gòu)和代碼變換，特別是對于涉及大規(guī)模稀疏矩陣 的計算。諸多學(xué)者對稀疏矩陣矢量乘法的優(yōu)化進(jìn)行了大量研究，并取得了許多研 究進(jìn)展【4,34,35】。存儲稀疏矩陣的常用格式之一是壓縮稀疏行(CSR)。顧名思義, 該方案壓縮稀疏矩陣并減少矩陣的存儲需求，并且通過僅對基于緩存的傳統(tǒng)微處 理器執(zhí)行必要的計

51、算來適當(dāng)?shù)貓?zhí)行。向量單元是現(xiàn)代處理器的重要加速部件，向量化計算是提升CPU單核性能的 重要方法。向量化部件出現(xiàn)后，計算機(jī)可以對多個操作數(shù)進(jìn)行矢量實現(xiàn)，即計算 機(jī)可以同時對多個操作數(shù)進(jìn)行操作。相比原來計算機(jī)完全采用標(biāo)量計算，矢量計 算計算機(jī)速度更快。然而，CSR格式采取的是存儲每一個非零元素的行坐標(biāo)，列 坐標(biāo)，數(shù)值來表示稀疏矩陣。稀疏矩陣在進(jìn)行S pMV(y = A*x)操作時，CSR格式 存儲的數(shù)值數(shù)組中的數(shù)據(jù)不具有連續(xù)性，且向量b是間接取址，這導(dǎo)致CSR格式 存儲的稀疏矩陣在進(jìn)行SpMV向量乘很難發(fā)揮出計算機(jī)底層向量化部件的加速的 性能。因此，本文可以從研究SpMV的向量化實現(xiàn)作為切入點，進(jìn)

52、一步提升SpMV 此重要核心程序的性能。通過矢量化技術(shù)，圍繞SpMV進(jìn)行了許多相關(guān)工作。串行庫OSKI和pOSKI 提供了一種基于BCSR的算法，該算法利用SIMD加速技術(shù)進(jìn)行優(yōu)化，作為對輸 入向量的不規(guī)則訪問進(jìn)行馴服，并利用其固有的重用。CSRD17采用對角結(jié)構(gòu)并 將矩陣分解為密集的對角矩陣。Xing Liu等【36】提供了一套基于x86的多核處理器 的ELLPACK格式的優(yōu)化技術(shù)。芳芳劉等印】提出了一種稱為CSRL的存儲方案， 其組合了相鄰的非零元素，這可以減少存儲器數(shù)據(jù)訪問的數(shù)量并且針對一些特殊 稀疏矩陣改進(jìn)向量單元上的SpMV性能。上述工作通過提供新的存儲格式來解決矢量化挑戰(zhàn)，但仍然是

53、一些限制。例 如，OSKI和pOSKI利用稀疏矩陣的特定屬性并且實現(xiàn)不同程度的效率。 ELLPACK存儲格式加速矩陣A的數(shù)據(jù)訪問，但往往產(chǎn)生向量x的不規(guī)則存儲器 訪問的高開銷。CSRD限于矩陣的對角線特性;劉芳芳提出的CSRL格式只適用于 具有足夠非零元素的稀疏矩陣，也就是說，在某些情況下，可能存在太多單個的 非零元素沒有向量化。3.2 CSR(r)向量化存儲格式3.2.1 CSR(r)存儲格式描述CSR(r)格式是在CSR格式上改進(jìn)的一種新的適合SpMV向量化運算的稀疏矩 陣存儲格式oCSR(r)格式通過采用的是存儲零元的方式連接稀疏矩陣同一行相鄰或 者相近的非零元素形成了具有連續(xù)性數(shù)據(jù)片段

54、。同時，將零元素添加到此連續(xù)數(shù) 據(jù)片段的末尾，以確保該數(shù)據(jù)片段能完全向量化。在該數(shù)據(jù)片段內(nèi)存儲的數(shù)據(jù)是 連續(xù)的，且可以完全向量化。這種存儲格式存儲的稀疏矩陣，在進(jìn)行SpMV操作 時，可以完全采用矢量計算，充分發(fā)揮了計算機(jī)底層向量化部件的加速作用。CSR(r) 存儲格式用存儲數(shù)據(jù)片段的方式代替了 CSR格式中存儲單個元素的存儲方式，更 好地挖掘了稀疏矩陣良好的數(shù)據(jù)局部性，使計算機(jī)向量化加速部件的在稀疏矩陣 進(jìn)行SpMV操作時取得更好的加速效果。CSR (r)格式的描述見下表3.2：表3.1矩陣A101000101100000001001000000000000011000100000000010

55、1000110000000100100000000000000000100100000000100010010000000000010010000000000010000111000001011000011010000000010000101000000000000101100000001011000001011000101100000101100101100000010110010110000001011表3.2稀疏矩陣A的CSR(r)存儲格式CSR(r)：0161112171117101511151413181018141818131121311212112121121Row Offset

56、s04610141618222428323638424650543.2.2 CSR(r)存儲算法介紹CSR(r)格式存儲的內(nèi)部數(shù)據(jù)連續(xù)的數(shù)據(jù)片段。該數(shù)據(jù)片段存儲稀疏矩陣內(nèi)部連 續(xù)的零元和所有的非零元素。同一行中相鄰或者相近的非零元素能否結(jié)合成一個 數(shù)據(jù)片段，需要通過兩個非零元素之間含有的零元素的個數(shù)來判斷。在CSR(r)格 式中本文定義用英文字母r表示兩個相鄰非零元素之間的零元素的數(shù)量，rmax是 r的最大允許值。與CSR存儲格式類似，CSR(r)格式引入了三個一維數(shù)組在存儲稀疏矩陣，這 三個一維數(shù)組一次命名為Values, Col(t), Row Offsets.Values存儲的是稀疏矩陣

57、內(nèi) 部的所有數(shù)據(jù)片段內(nèi)數(shù)據(jù)的數(shù)值。Col(t)依次記錄了稀疏矩陣內(nèi)部每個數(shù)據(jù)片段首 個非零元素的在稀疏矩陣內(nèi)的列坐標(biāo)，以及該完成該數(shù)據(jù)片段的SpMV操作所需 的向量化次數(shù)。Row Offsets記錄的是每行首個非零元素在Col(t)數(shù)組中的行偏移 (即某一行的第一個元素在col(t)里面的起始偏移位置)。存儲格式的詳細(xì)算法如 下所示：(1 )將0寫入Row Offsets數(shù)組內(nèi)，作為Row Offsets數(shù)組的第一個數(shù)。從左至右按行讀取仍未被讀到的稀疏矩陣中的數(shù)據(jù)，將遇到的首個非零 元素寫入Values數(shù)組，將該非零元定義為當(dāng)前非零元。然后則記錄下該非零元素 在稀疏矩陣中的列坐標(biāo)column

58、index,數(shù)據(jù)片段長度length記為1。查找與當(dāng)前非零元素距離最近的非零元并將該非零元定義為最近非零元 素，并記錄此兩個非零元素之間含有的零元素的數(shù)量r。若rv=rmax,則將這兩個相鄰非零元素之間的零元素和最近非零元按順國防科學(xué)技術(shù)大學(xué)研究生院工程碩士學(xué)位論文 序?qū)懭雟alues數(shù)組，將最近非零元定義為當(dāng)前非零元。然后修改數(shù)據(jù)片段長度的 值。數(shù)據(jù)片段長度length的數(shù)值變?yōu)楫?dāng)前數(shù)據(jù)片段長度length加上兩個非零元素 之間含有的零元素的數(shù)量r再加1。跳到(2)繼續(xù)搜索。若rrmax,在當(dāng)前連續(xù)數(shù)據(jù)片段的末尾附加零元素，直到該片段的長度 可被s整除(s依次向量化計算所需要的操作數(shù)的個數(shù))

59、。接下來修改數(shù)據(jù)片段長 度的值。數(shù)據(jù)片段長度length的數(shù)值變?yōu)楫?dāng)前數(shù)據(jù)片段長度length加上添加零元 的個數(shù)。然后將column index, length/s依次記入Col(t)數(shù)組。假設(shè)最近非零元不 是該稀疏矩陣此行的最后一個非零元，則跳到(1)繼續(xù)進(jìn)行搜索；假設(shè)最近非零元 為該稀疏矩陣此行最后一個非零元，則需要判斷該行是否為稀疏矩陣的最后一行。 假設(shè)此行不是稀疏矩陣的最后一行，則將此時Col(t)數(shù)組長度寫入Row Offsets數(shù) 組中；假設(shè)此行不是稀疏矩陣的最后一行，則將此時Col(t)數(shù)組長度寫入Row Offsets數(shù)組中，并結(jié)束搜索。在稀疏矩陣進(jìn)行SpMV(y = A *

60、 x)操作時，數(shù)據(jù)片段中內(nèi)的數(shù)據(jù)都是連續(xù)的， 且完全可以采用向量化計算完成SpMV操作。數(shù)據(jù)片段內(nèi)的數(shù)據(jù)對稠密向量x的 訪存是規(guī)則的、連續(xù)的。因此在稀疏矩陣進(jìn)行SpMV(y = A*x)操作時可以直接采 用向量化讀取和計算.另外由于數(shù)據(jù)采取的是以連續(xù)的數(shù)據(jù)片段存儲的方式，大大 減少了對于向量x的間接取址的次數(shù)，改善了因間接訪存帶來的訪存帶寬利用率 低這一問題。由于稀疏矩陣在進(jìn)行SpMV(y = A*x)操作時，可以完全采用SIMD 指令一次對多個操作數(shù)進(jìn)行計算，大大減少了浮點計算的時間。rmax的確定CSR(r)格式存儲的內(nèi)部數(shù)據(jù)連續(xù)的數(shù)據(jù)片段。該數(shù)據(jù)片段存儲稀疏矩陣內(nèi)部連 續(xù)的零元和所有的非