丁文隆 汪承寧,2 童 薇,2

1(華中科技大學計算機科學與技術(shù)學院 武漢 430074)

2(武漢光電國家研究中心(華中科技大學) 武漢 430074)

在科學與工程領(lǐng)域中，許多復雜的模型都會用線性系統(tǒng)Ax=b的形式表達[1]，其中，A通常是一個龐大的高精度浮點稀疏矩陣[2].在這樣的系統(tǒng)上進行求解會消耗大量時間與計算資源[3].目前主流的求解大規(guī)模稀疏線性系統(tǒng)的方法是克利洛夫(Krylov)子空間方法[1].該方法的求解過程涉及大量的矩陣向量乘法(matrix-vector multiplication, MVM)計算，所以改進MVM的計算模式是節(jié)省運算能耗與運算資源的關(guān)鍵途徑.

近年來，關(guān)于用憶阻陣列進行原位MVM運算的有關(guān)加速器被不斷地提出.最先提出的是一類用于執(zhí)行機器學習和圖形處理任務(wù)的憶阻加速器，它們通過向憶阻陣列施加電壓，在陣列上原位地執(zhí)行MVM運算[4-10].然而，這些憶阻加速器只提供8～16 b的計算，顯然不支持一些以高精度浮點數(shù)為主的計算應(yīng)用.于是，最新的研究提出了將IEEE-754雙精度浮點數(shù)部署在陣列上的方式[11].該研究將浮點數(shù)的53 b尾數(shù)(包括前導1)分位片地部署到不同的陣列上，最后通過乘加縮減樹整合不同位片的結(jié)果，實現(xiàn)高精度浮點數(shù)運算.

然而，目前仍然沒有任何工作提出在一個系統(tǒng)內(nèi)，能夠為不同精度的應(yīng)用提供計算能耗優(yōu)化的方法.本工作致力于提出一個基于憶阻陣列的模擬MVM運算系統(tǒng)，既能夠為精度較高的應(yīng)用提供無損的浮點數(shù)計算模式，又能夠讓較低精度的應(yīng)用執(zhí)行低能耗開銷計算.具體地，本文工作的貢獻主要有4個方面：

1) 設(shè)計了一種自選尾數(shù)壓縮機制，對于某些低精度的求解應(yīng)用，可以在陣列運算中選擇性地忽略激活若干權(quán)值較小的低位陣列，從而保證在滿足具體任務(wù)的求解精度的前提下，減少運算陣列以及外圍電路的能耗.同時提出了一種動態(tài)的對齊位優(yōu)化機制，摒棄原有的靜態(tài)對齊位的設(shè)定模式，根據(jù)矩陣實際指數(shù)范圍的要求來設(shè)置參與運算的對齊位陣列，減少冗余對齊位帶來的能耗.

2) 針對于尾數(shù)壓縮和對齊位優(yōu)化策略，提出了一種葉子結(jié)點數(shù)可變的流水乘加縮減樹的結(jié)構(gòu).解決了激活運算的陣列數(shù)不確定所帶來的偏移策略失效以及無法進行流水的問題.這種新型的乘加縮減樹結(jié)構(gòu)能夠在激活任意數(shù)量的尾數(shù)位和對齊位陣列的情況下，都正常地進行流水求和運算.

3) 基于現(xiàn)有工作中的分塊映射與片位部署的思想以及上述所提到的優(yōu)化改進策略，設(shè)計完成整個模擬憶阻MVM原位運算系統(tǒng).

4) 將帶優(yōu)化策略的MVM原位模擬運算系統(tǒng)集成到高性能線性代數(shù)的算法框架中進行求解計算.評估在不同的求解精度下，系統(tǒng)中關(guān)鍵能耗指標的減少程度.

1 背景知識

由于憶阻器陣列存內(nèi)運算具有高能效、高并行性等優(yōu)勢，受到國內(nèi)外學者廣泛地關(guān)注.到目前為止，用憶阻器陣列進行MVM運算的方法已經(jīng)被廣泛地研究，包括從用于機器學習模型的憶阻加速器的研究到應(yīng)用于科學計算的憶阻加速器的研究.

1.1 用于機器學習和圖形處理的憶阻加速器

近年來，隨著機器學習工作的流行性日益增加，學術(shù)界[12]和工業(yè)界[13]都提出了許多有關(guān)于機器學習的專用加速器的建議.這些加速器主要為機器學習模型提供MVM運算，通過加速這一過程的運算，極大程度提升機器學習模型的訓練速度.

Fig. 1 Illustration of memristive crossbars computation圖1 憶阻陣列計算示意圖

Fig. 2 Array organization in a memristive cluster圖2 憶阻存儲簇中的陣列組織

2) 位切片(bit slicing)技術(shù).由于模擬器件的精度受到非理想因素的限制，憶阻存儲單元無法精確地表達映射值[14].針對這個問題，現(xiàn)有的工作提出了位切片技術(shù)，這種技術(shù)通過將矩陣元素按照位片映射到多個陣列上來減少對器件精度的依賴[4-7].一個位切片的例子，這個矩陣將被映射到3個二進制存儲陣列中：

(1)

注意，在實際計算中，對于到來的電壓Vi采取這種位切片技術(shù).

3) 乘加縮減樹.如圖1所示,為了將各位片的點積結(jié)果進行求和運算，需要乘加縮減樹進行整合[4-5]，得到MVM運算的最終結(jié)果.乘加縮減樹中有2類操作：①加法操作，即不同位片得到的結(jié)果求和；②左移操作，父結(jié)點整合左右子樹時，高位片子樹需要進行左移操作才能和低位子樹對齊相加.一般地，乘加縮減樹都是葉子結(jié)點數(shù)為2的整數(shù)次方的滿二叉樹，在這種樹結(jié)構(gòu)下，整合移位時將高位子樹左移2i-1b即可(i為父結(jié)點到葉子結(jié)點的距離).

1.2 用于科學計算的憶阻加速器

1.1節(jié)提到的憶阻MVM加速器是針對機器學習模型來設(shè)計的，只提供8～16 b的計算[4-9].雖然這個運算精度大概率不會影響機器學習任務(wù)的推理準確度，但是對于高精度的科學計算求解來說，這種運算精度是遠遠不夠的[15].

1) 浮點數(shù)的憶阻陣列部署.為了解決科學計算中的精度問題，F(xiàn)einberg等人[11]首次提出了在憶阻加速器上部署科學計算.在IEEE-754標準中，雙精度浮點數(shù)由53 b的二進制尾數(shù)(包括前導1)、11 b的二進制指數(shù)和一個符號位表示[16].這項工作首先提出了如何將53 b的尾數(shù)部署到憶阻陣列上進行MVM運算，它運用了多個存儲陣列，每一個陣列代表矩陣中浮點數(shù)的一個位片，形成一個存儲簇；然后再將不同的陣列通過乘加縮減樹連接起來，如圖2所示：

存儲簇用于計算向量的一個位片與矩陣元素的乘積.在完成了位片映射后，進行3階段的運算：1)電壓加載在到模塊A中的陣列上，列電流存儲在采樣—保持緩沖區(qū)中，該電流代表當前向量位片和對應(yīng)的矩陣位片每一行的點積之和;2)在模塊B中，選擇器選擇采樣—保持緩沖區(qū)中代表相同行不同位片的電信號，經(jīng)過ADC轉(zhuǎn)化為數(shù)字信號后，送入乘加縮減樹相加;3)模塊C中乘加縮減樹對模塊B中輸入的不同位片結(jié)果進行整合運算.

經(jīng)過這3個階段，得到了矩陣的1行與向量位片的點積結(jié)果.為獲得最終結(jié)果，需要循環(huán)執(zhí)行第2階段與第3階段，獲得所有矩陣行與向量位片的乘積.

2) 異構(gòu)陣列集合.Feinberg等人[11]的工作還提出了一種異構(gòu)矩陣分塊思想.采用不同大小的塊去捕捉稀疏矩陣中非0元素的密集區(qū)域，并為每個大小的塊設(shè)定一個非0元素閾值，若達到閾值就將該分塊映射到存儲簇上[11].對比使用單一陣列大小的工作[4-5,9]，這種異構(gòu)設(shè)計極大程度上增加了陣列非0元素密度，減少了陣列映射開銷，保證了陣列的并行性和能源效率.

圖3展示了矩陣元素分塊的一種示例.其中，數(shù)據(jù)來源是SuiteSparse矩陣數(shù)據(jù)集合[2]中的685_bus矩陣.其中，圖3(a)展示了該矩陣非0元素的分布情況，圖3(b)展示了分塊之后的結(jié)果.分別采用了32×32，16×16，8×8，4×4這4種不同大小的矩陣塊對原矩陣進行分塊，閾值分別為128，32，8，2.

Fig. 3 Example of matrix blocking (685_bus matrix)圖3 矩陣元素分塊示例(685_bus矩陣)

Fig. 4 Example of alignment bits deployment in cluster圖4 存儲簇中對齊位部署示例

3) 存儲簇中對齊位的部署.完成了矩陣分塊之后，需要將每一塊的浮點數(shù)都映射到存儲簇中.然而同一矩陣塊中元素的二進制指數(shù)不盡相同.為了使得不同元素的尾數(shù)按照指數(shù)對齊，需要以塊中指數(shù)最大的元素作為基準，對較小指數(shù)的元素實行向右偏移部署的策略，然后在空余位補充0.

圖4展示了存儲簇中浮點數(shù)格式部署的一個例子.需要部署10.5，6.5，0.3這3個浮點數(shù)，其對應(yīng)的二進制指數(shù)分別為3，2，-2.部署時以10.5對應(yīng)的指數(shù)作為基準，6.5和0.3的部署分別向右偏移1 b和5 b，并將空余位置填充0.

4) 定點硬件上執(zhí)行浮點運算.在系統(tǒng)實現(xiàn)過程中，圖4所述的利用填充0實現(xiàn)的對齊步驟是在浮點數(shù)映射過程中實現(xiàn)的，該過程處于預處理階段.在陣列運算過程中，不會再出現(xiàn)動態(tài)的移位對齊等步驟，整個運算是在定點硬件上進行的.所以，同一個存儲簇中的陣列運算，僅涉及到定點運算.通過二進制分片策略與預處理填充對齊位的策略，實現(xiàn)了在定點硬件上執(zhí)行浮點運算.而在整合不同存儲簇的運算結(jié)果時，會涉及到定點數(shù)轉(zhuǎn)浮點數(shù)，浮點數(shù)相加等純浮點運算.故所構(gòu)造的系統(tǒng)是一種以定點硬件運算為主，純浮點運算為輔的系統(tǒng).

2 觀察與動機

之前的工作所使用的都是固定位數(shù)的浮點數(shù)尾數(shù)部署方式，這種部署模式可能會產(chǎn)生額外能耗開銷.本節(jié)將討論這個問題以及對應(yīng)的優(yōu)化策略.

2.1 為不同的應(yīng)用優(yōu)化計算能耗開銷

在科學與工程領(lǐng)域中，雖然多數(shù)的線性模型的矩陣元素都是浮點數(shù)，但是不同應(yīng)用所需的求解精度各不相同.例如，求解網(wǎng)頁排序算法(PageRank)、求解線性回歸等大規(guī)模數(shù)理統(tǒng)計模型時，所需要的求解精度都不高[17-18].已有的研究顯示，在執(zhí)行機器學習、圖形處理任務(wù)時，只使用8～16 b精度就能達到令人滿意的準確率[4-9].但是，對于一些由偏微分方程離散化得到的大型稀疏性方程組而言，需要的求解精度比較高[1]，其代表領(lǐng)域有航天航空領(lǐng)域以及量子力學領(lǐng)域等.然而，執(zhí)行高精度的求解運算必然會消耗更多的能量，這對于求解精度需求較低的線性系統(tǒng)來說是不劃算的.

現(xiàn)有工作中，高精度浮點數(shù)采用的是固定位片數(shù)的部署形式[11]，在每一個存儲簇中，有53 b的二進制尾數(shù)陣列、64 b的對齊位陣列以及9 b的校驗位陣列.在實際運算時，運算陣列與其外圍電路的能耗都與激活計算的陣列個數(shù)呈正相關(guān)，也就是說，隨著參與運算的陣列個數(shù)增多，計算能耗也會劇烈增加.然而，上述的浮點數(shù)部署模式是針對于通用IEEE-754雙精度浮點數(shù)格式設(shè)計的，對于所有的雙精度浮點應(yīng)用都能夠保證計算精度.對于低精度的求解應(yīng)用而言，并不需要激活所有的53 b尾數(shù)、64 b對齊位進行計算就能夠達到允許精度范圍內(nèi)的解，所以按照原有固定方式計算必然會產(chǎn)生不必要的能耗.

2.2 尾數(shù)壓縮與對齊位優(yōu)化策略

為了解決2.1節(jié)所述的為不同應(yīng)用優(yōu)化計算能耗開銷的問題，在本節(jié)中將論述一種尾數(shù)壓縮和對齊位優(yōu)化策略，動態(tài)地決定實際運算時激活的尾數(shù)以及對齊位陣列數(shù)量.

1) 尾數(shù)壓縮策略.在實際計算時，可以根據(jù)具體的應(yīng)用精度，自行選擇參與運算的尾數(shù)位數(shù)，從而大幅度減少運算陣列帶來的能量消耗.例如，如果選擇25 b尾數(shù)參與運算，在實際運算時只會激活高25 b尾數(shù)陣列進行運算，而不會對剩下的低位陣列(低28 b的陣列)施加電壓，也就不會產(chǎn)生這些陣列對應(yīng)的運算陣列和外圍電路能耗.

2) 對齊位優(yōu)化策略.在之前的有關(guān)于位片型憶阻陣列科學計算的系統(tǒng)設(shè)計中[11]，對齊位是固定的64 b.然而，并不是所有的應(yīng)用都需要多如64 b的指數(shù)范圍.為了節(jié)省計算開銷，采取一種動態(tài)的對齊位優(yōu)化策略，只部署指數(shù)范圍內(nèi)的對齊位.例如，一個矩陣塊中非0元素的二進制指數(shù)范圍是1～20，即只部署19 b的對齊位陣列即可，這樣便會一定程度上減少對齊位冗余帶來的陣列計算消耗.

2.3 葉子結(jié)點數(shù)可變的流水乘加縮減樹

由于本文所提出的系統(tǒng)應(yīng)用了尾數(shù)壓縮和對齊位優(yōu)化策略，導致不同應(yīng)用的計算中，尾數(shù)和對齊位的位數(shù)之和是不固定的.在這種情況下，進行計算的乘加縮減樹中葉子結(jié)點的個數(shù)是可變的.所以，為了解決尾數(shù)壓縮和對齊位優(yōu)化策略帶來的葉子結(jié)點數(shù)可變的問題，需要設(shè)計一種新型乘加縮減樹結(jié)構(gòu)，來完成存儲簇中不同位片結(jié)果的整合.

1)乘加縮減樹中的層級流水概念.假設(shè)目前存儲簇擁有128個憶阻陣列，陣列大小為n×n，即對應(yīng)的乘加縮減樹為7層.每加載一次向量，該樹就需要計算n次，以得到矩陣塊每一行和向量相乘的結(jié)果.然而，由于樹結(jié)點之間的計算可以并行操作，且同層之間沒有數(shù)據(jù)依賴，可以把乘加縮減樹的每一層看作一個流水部件，構(gòu)建一個7級的流水線.假設(shè)每個結(jié)點的計算時間為單位時間，這樣計算n個結(jié)果的時間為6+n.

2) 葉子結(jié)點數(shù)可變的流水.在以往固定位片的陣列部署模式下，葉子結(jié)點數(shù)是確定的，且數(shù)量往往是2的整數(shù)次方(若不為2的整數(shù)次方，則往往將葉子結(jié)點數(shù)補滿到2的整數(shù)次方).在此情況下，所構(gòu)造的樹為滿二叉樹.現(xiàn)在由于葉子結(jié)點數(shù)可變，會導致整棵樹不為原有的滿二叉樹形式，導致2個問題：①由于不是滿二叉樹，各葉子結(jié)點到根結(jié)點的距離不盡相同.若采用流水策略，各葉子結(jié)點數(shù)值無法同時到達根結(jié)點，產(chǎn)生計算錯誤.②若所計算的乘加縮減樹不為滿二叉樹，則前面背景介紹中論述的高位子樹偏移2i-1(i為父結(jié)點到葉子結(jié)點的距離)與低位子樹對齊相加的策略就會失效，使用該方法來計算高位子樹的偏移量會帶來錯誤.為了解決這2個問題，需要提出葉子結(jié)點數(shù)可變的流水乘加縮減樹(在3.3.3節(jié)詳細論述).

3 基于異構(gòu)陣列的自選尾數(shù)壓縮系統(tǒng)設(shè)計

本節(jié)將展示一個基于異構(gòu)陣列集合，擁有尾數(shù)壓縮和對齊位優(yōu)化策略，面向科學計算的MVM系統(tǒng).

3.1 系統(tǒng)總體結(jié)構(gòu)設(shè)計

如圖5所示，系統(tǒng)主要由預處理模塊以及陣列運算和整合模塊構(gòu)成.

1) 預處理模塊.預處理部分主要負責處理原始矩陣，并提供陣列映射方案.具體地，該部分將對原始矩陣依次進行分塊、分片、執(zhí)行尾數(shù)壓縮與對齊位優(yōu)化策略，最終獲得能直接映射在憶阻陣列上的0-1矩陣，并為后續(xù)計算提供陣列激活策略.

2) 陣列運算和整合模塊.陣列運算和整合部分主要負責對于到來的向量x，通過陣列計算以及輔助處理子模塊的整合，獲得最終Ax的結(jié)果.具體地，需要先將向量位片按時序加載在存儲簇上，然后在陣列上執(zhí)行矩陣向量乘法運算，再在輔助處理子模塊中對不同存儲簇的結(jié)果以及未分塊元素進行整合，最后得到矩陣向量的乘法結(jié)果b.

Fig. 5 Overall structure design of the system圖5 系統(tǒng)總體結(jié)構(gòu)設(shè)計

Fig. 6 Design of the pre-processing step圖6 預處理過程設(shè)計

3.2 預處理模塊的詳細設(shè)計

圖6展示了預處理模塊主要任務(wù)的設(shè)計流程：

1) 矩陣分塊與分片.由用戶輸入最大塊的邊長L和最大塊中允許的最少非0元素個數(shù)閾值p.分別運用邊長為L，L/2，L/4，L/8的塊來采集矩陣中非0元素密集的區(qū)域，其對應(yīng)的非0元素閾值分別為p，p/4，p/16，p/64.在分塊過程中，先用邊長為L的塊去鋪排整個矩陣，根據(jù)閾值p去捕獲達到密度的塊，若密度達不到，則用4個邊長為L/2的塊去鋪排這個邊長為L的矩陣塊，捕獲非0元素數(shù)量超過閾值p/4的區(qū)域；以此類推，用L/4，L/8的塊去鋪排.若有非0元素沒有被任何塊捕獲，則其會進入未分塊元素集合，等待輔助處理子模塊統(tǒng)一整合處理.

對于每個浮點數(shù)而言，需要獲得其IEEE-754雙精度浮點數(shù)表示形式(或者其科學記數(shù)法表示形式)，進而提取其53 b尾數(shù)(包括前導1)、其指數(shù)以及其符號位.對于同一個矩陣塊，需要對其中所有的非0元素執(zhí)行此操作，獲得一個具有(sign,significand,expval)三元組元素的矩陣，其中sign表示符號位；significand表示有效數(shù)，即帶前導1的共53 b尾數(shù)；expval表示對應(yīng)二進制指數(shù).

2) 尾數(shù)壓縮與對齊位優(yōu)化.為實現(xiàn)自選尾數(shù)壓縮策略，可以在系統(tǒng)上增加一個尾數(shù)選擇壓縮的接口mantissa_bits，其意義是用戶設(shè)定的尾數(shù)保留位數(shù)，保留高mantissa_bits.

為實現(xiàn)冗余對齊位的優(yōu)化，需要在矩陣分塊之后，對矩陣塊中每個非0元素的指數(shù)進行簡單處理，獲得該塊中最大的二進制指數(shù)maxexp和最小的二進制指數(shù)minexp，本系統(tǒng)參考矩陣塊中指數(shù)的實際范圍來進行對齊位陣列的激活.同時，本系統(tǒng)同時還會給出一個可供用戶設(shè)置的對齊位數(shù)閾值max_alig.對齊位數(shù)的設(shè)置形式：

alig_bits=min(maxexp-minexp,max_alig).

(2)

在一般情況下，對齊位數(shù)alig_bits=maxexp-minexp，這樣做就可以用最少的對齊位數(shù)覆蓋整個矩陣塊的指數(shù)范圍，相對于固定的對齊位方案而言減少了不必要的計算能耗.給定用戶設(shè)置的閾值max_alig的目的是為了防止在矩陣塊中有個別指數(shù)極小的元素，導致maxexp和minexp差距過大，需要用到過多的對齊位，白白耗費陣列資源.當然，如果具體的問題對求解精度要求不高但對于節(jié)省能耗的要求較高，用戶也可以根據(jù)節(jié)省能耗和求解精度的具體要求選擇較小的對齊位閾值，這一點是靈活的.

如果矩陣塊中有元素的指數(shù)超過了閾值所限制的范圍，則在預處理步驟中這些元素會重新被歸類到未分塊元素中，而不會在后續(xù)進行陣列映射.尾數(shù)壓縮和對齊位優(yōu)化的過程都在預處理步驟進行，即在陣列映射前完成.可以通過這2個步驟確定每一個存儲簇中具體需要激活使用的陣列數(shù)量.

圖7展示了尾數(shù)壓縮與對齊位優(yōu)化策略在預處理模塊中整體的執(zhí)行流程：

Fig. 7 Design of mantissa compaction and alignment optimization圖7 尾數(shù)壓縮與對齊位優(yōu)化的設(shè)計

3) 映射.對于每一個矩陣塊，需要將其映射到對應(yīng)的存儲簇上.假設(shè)矩陣塊大小為N×N，一共需要映射的陣列個數(shù)為M個.則對應(yīng)的存儲簇由M個N×N的憶阻存儲陣列組成，將得到的矩陣塊尾數(shù)二進制串按照從高位到低位映射到這M個存儲陣列上即可.

3.3 陣列運算與整合模塊的詳細設(shè)計

本節(jié)主要闡述，在一個已經(jīng)完成映射的憶阻異構(gòu)陣列集合中，從加載向量x，到計算出最后的結(jié)果Ax，這一完整過程的設(shè)計.

3.3.1 陣列運算與整合模塊結(jié)構(gòu)設(shè)計

圖8展示了陣列計算和整合模塊的設(shè)計結(jié)構(gòu).其中，在每個存儲簇的輸入緩沖區(qū)中，對于到來的向量x進行分片和對齊，并以時序展開，依次輸入存儲簇中進行運算.存儲簇通過陣列點積的計算以及乘加縮減樹的計算，將向量x的當前位片值與當前矩陣塊的點積結(jié)果存在輸出緩沖中.接下來的工作由輔助處理子模塊完成，包括整合同一存儲簇不同向量位片的結(jié)果，整合不同存儲簇及未分塊元素的結(jié)果.全局緩存區(qū)用于存儲必要中間值以及最后的結(jié)果.

Fig. 8 Crossbar computation and integration module structure圖8 陣列運算與整合模塊的結(jié)構(gòu)

3.3.2 負數(shù)的處理

在矩陣A和向量x中，會出現(xiàn)負浮點數(shù)，本節(jié)將對這2種負數(shù)的處理方式進行論述.

1) 矩陣塊中的負數(shù)處理.為了處理矩陣A映射到存儲陣列上的負數(shù)值，可以將一個矩陣塊中的正數(shù)值和負數(shù)值分開映射(以三元組中的符號位元素作為判斷標準)，然后在后續(xù)整合模塊的時候相減即可.這樣，同一存儲簇中就要維護2個陣列集合，分別對應(yīng)該矩陣塊的正非0元素和負非0元素.在存儲簇c中，正存儲陣列集合和負存儲陣列集合分別為cpos和cneg.

2) 向量的負數(shù)處理.對于向量x來說，可以用負電壓來代表x中的負數(shù)元素.但是這會導致一個問題，即產(chǎn)生的電流值會有負數(shù)，進而經(jīng)過ADC轉(zhuǎn)換之后會有負數(shù)值進入乘加縮減樹.然而本工作中實現(xiàn)的乘加縮減樹沒有集成關(guān)于減法(或補碼)的運算，這會導致其無法處理負數(shù)值.為處理這個問題，可以給進入乘加縮減樹的電流轉(zhuǎn)換值都增加一個數(shù)值為N的偏移，其中N為該存儲陣列的大小.這樣便能夠?qū)㈦娏鬓D(zhuǎn)換的取值范圍從-N～N轉(zhuǎn)換到0～2N.然后，由于在整合模塊中正數(shù)陣列的值要和負數(shù)陣列的值一一相減，所以增加的偏移值會在運算中剛好抵消，而不用額外操作來消除偏移值.

3.3.3 葉子結(jié)點數(shù)可變的流水乘加縮減樹的設(shè)計

如2.3節(jié)所述，相對于普通的葉子結(jié)點數(shù)固定的滿二叉樹形式，葉子結(jié)點數(shù)可變的樹形式存在2個問題：移位策略失效、流水過程計算出錯.下面將針對這2個問題提出正確的樹的結(jié)構(gòu)設(shè)計.

1) 解決移位策略失效的問題

在乘加縮減樹是滿二叉樹時，假設(shè)當前結(jié)點到葉子結(jié)點的距離為i，現(xiàn)在該結(jié)點需要整合左右2棵子樹的元素，即將右子樹(規(guī)定右子樹為高位子樹)的值左移2i-1位，和左子樹相加即可.例如，某個葉子結(jié)點數(shù)為8的滿二叉樹，這棵樹共有3層.考慮其最高位元素，它相對于最低位的元素需要左移7 b.它在從葉子結(jié)點到達根結(jié)點前需要經(jīng)歷i=1，i=2，i=3這3層中對應(yīng)結(jié)點的共3次移位，總共移動21-1+22-1+23-1=1+2+4=7 b，剛好符合移位要求.對于編號為奇數(shù)的葉子結(jié)點也同理，例如第5個元素，其需要左移4 b.它在到達根結(jié)點之前需要經(jīng)歷i=3層結(jié)點的1次移位(因為在i=1與i=2層中其為左子樹不需移位)，即總共移動23-1=4 b，符合要求.

現(xiàn)在由于應(yīng)用了尾數(shù)壓縮與對齊優(yōu)化策略，葉子結(jié)點的數(shù)量不會固定為2的整數(shù)次方，乘加縮減樹不再為滿二叉樹，會導致原有移位策略失效.若按照傳統(tǒng)方法，從根結(jié)點開始構(gòu)造該1棵葉子結(jié)點為6的盡量平衡的二叉樹的話，結(jié)果如圖9(a)所示.其中，結(jié)點A的實際偏移量應(yīng)為5，而按照原有偏移策略算得的偏移量為22-1+23-1=2+4=6，產(chǎn)生錯誤.而如圖9(b)模式構(gòu)造樹，符合原有的計算模式，其偏移為21-1+23-1=1+4=5，偏移值正確.

Fig. 9 Comparison of trees with wrong and correct shift computation modes圖9 擁有錯誤與正確移位計算模式的樹的對比

當前的目標是設(shè)計擁有正確移位模式的乘加縮減樹的構(gòu)造策略.即需要對任意數(shù)量的葉子結(jié)點，找到一種通用策略，構(gòu)造如圖9(b)所示的正確移位模式的樹.受到偏移模式以及一些從葉子結(jié)點開始構(gòu)造的樹(如霍夫曼樹)的啟發(fā)，可以選擇由葉子結(jié)點開始構(gòu)造這樣1棵具有正確模式的計算樹.具體地，從葉子結(jié)點層開始，兩兩相鄰的結(jié)點為一組，構(gòu)造下一層父結(jié)點；若本層有落單結(jié)點(如果結(jié)點數(shù)為奇數(shù)則最后一個結(jié)點為落單結(jié)點)，則該落單結(jié)點和新構(gòu)造的父結(jié)點一起組成一個新的集合，繼續(xù)構(gòu)造接下來一層的父結(jié)點.如此循環(huán)地執(zhí)行該構(gòu)造策略，當新的集合中的結(jié)點數(shù)為1時，樹的構(gòu)造停止，該結(jié)點為根結(jié)點.

如此構(gòu)造乘加縮減樹，便可用i的值計算出這些結(jié)點所在的實際層數(shù)，保證每個結(jié)點在運算中能獲得正確的偏移量，到達根結(jié)點的時候獲得了正確的總移位值.圖10左圖用這種方法構(gòu)造了1棵有11個葉子結(jié)點的具有正確移位計算模式的乘加縮減樹.可以看到，在i=0這一層中，結(jié)點A落單，它與i=1層的5個結(jié)點一起，構(gòu)造了i=2這一層的結(jié)點，然而，在這一過程中，結(jié)點A所對應(yīng)的i值不變，仍然為0.同理，對于結(jié)點B，其為i=2這一層的落單結(jié)點，它與i=3這一層的結(jié)點一起構(gòu)造了i=4這一層的結(jié)點，同時，結(jié)點B對應(yīng)的i值仍然不變，即i=2.

所以，對于結(jié)點A的計算過程而言，其到達根結(jié)點C需要經(jīng)過結(jié)點B和結(jié)點C的移位計算，其中在結(jié)點B處移位22-1=2 b，在結(jié)點C處移位24-1=8 b，一共移位2+8=10 b，符合要求.

2) 解決無法進行流水操作的問題

想要在改進后的乘加縮減樹中執(zhí)行層級流水策略，會因為各葉子結(jié)點到根結(jié)點的距離不同而導致計算錯誤.如圖10左圖所示，葉子結(jié)點D到根結(jié)點C的路徑長為4，而葉子結(jié)點A到根結(jié)點C的路徑長為2.由于每進行一次流水操作都要對葉子結(jié)點加載一批新的數(shù)值，所以，當結(jié)點A經(jīng)過2次流水操作到達結(jié)點C時，同一時刻輸入的結(jié)點D值還沒到達結(jié)點C；當結(jié)點D值到達結(jié)點C時，此時從結(jié)點A到達結(jié)點C的值已經(jīng)是2個周期之后加載在結(jié)點A的值.所以，在計算過程中始終面臨著數(shù)據(jù)達到不同步的錯誤.

① 結(jié)點維護先進先出隊列.為了解決無法流水的問題，可以讓每一個結(jié)點都維護一個先進先出隊列，用隊列的長度來彌補某些葉子結(jié)點到根結(jié)點的距離比其他葉子結(jié)點短的問題.圖10右圖展示了圖10左圖中乘加縮減樹對應(yīng)的隊列結(jié)構(gòu).可以發(fā)現(xiàn)，由于結(jié)點A，B，C這一計算路徑長度為2，這個值要小于滿二叉樹部分的計算路徑(例如從結(jié)點D經(jīng)歷3個中間結(jié)點到結(jié)點C，其路徑長為4)，所以需要結(jié)點A和結(jié)點C的隊列長度為2，其余結(jié)點的隊列長度為1.當父結(jié)點需要整合左右子樹的值時，只需要從左右子樹取出其隊首的值進行操作即可，在圖10右圖中實線代表該子樹需要移位后才能進入父結(jié)點進行相加操作，虛線代表可以直接進入父結(jié)點進行加法操作.具體地，若當前需要計算的父結(jié)點的取值隊列為cur，左右子樹的取值隊列分別為left和right，則計算整合的過程描述為

cur.push(left.pop+right.pop?2i-1).

(3)

Fig. 10 Example of the correct pipelined multiply-and-add tree圖10 正確形式的流水乘加縮減樹的結(jié)構(gòu)示例

圖10右圖描述了第6次向樹的葉子結(jié)點中加載數(shù)值時，樹各結(jié)點中隊列的存儲情況.其中，ti代表該結(jié)點隊列中存儲值在數(shù)據(jù)流中的時間順序，例如t4代表該值是由第4次加載在葉子結(jié)點中的數(shù)據(jù)流經(jīng)過乘加縮減樹對應(yīng)位置計算得到的值.可以看到，在結(jié)點A和結(jié)點C中分別需要存儲t4，t5和t2，t3的值，以保證在i=2層計算時，結(jié)點B要從結(jié)點A中獲得t4時間順序的值，且保證在i=4層計算時，結(jié)點C要從結(jié)點B中獲得t2時間順序的值.

Fig. 11 Auxiliary processing module workflow圖11 輔助處理子模塊工作流程

② 維護正確的隊列長度與存儲值.現(xiàn)在設(shè)計的關(guān)鍵是要讓所有葉子結(jié)點的計算路徑上都存儲相應(yīng)長度的隊列，這樣層級流水才能取得正確的值.由此，在流水線沒有充滿的時候，不能夠一次性地激活所有流水層級，而是一級一級地激活.以圖10為例，第1次加載數(shù)據(jù)時，需要激活層i=0計算；第2次加載數(shù)據(jù)時，激活層i=1，即此時i=0和i=1這2層參與計算，直到最后第5次加載數(shù)據(jù)時，激活流水所有層級計算；并且在隨后的計算中，都是每一次加載數(shù)據(jù)，然后流水的所有層級一起計算.這樣一來，在激活i=1的時候，只有i=0和i=1這2個層級在計算，結(jié)點A按照順序儲存了2個待加值；同理，在激活層i=2和層i=3后，結(jié)點B按順序儲存了結(jié)果，符合預期.

3.3.4 輔助處理子模塊的運算設(shè)計

圖11展示了該過程中輔助處理子模塊的處理步驟與全局緩沖區(qū)中的內(nèi)容變化：

經(jīng)過了陣列運算與乘加縮減樹整合之后，每個存儲簇的緩沖輸出中都儲存了矩陣塊與向量位片的乘積結(jié)果，下面稱這個結(jié)果為關(guān)于該向量位片的部分和.

輔助處理子模塊將對這些結(jié)果進行整合，獲得最后的Ax結(jié)果.具體解釋：

1) 設(shè)該存儲簇c存儲的矩陣為Ac，其中正數(shù)與負數(shù)存儲陣列集合為cpos和cneg，向量第j位片對應(yīng)的部分和為Ac×xj.將不同位片的部分和相加，得到當前矩陣塊正數(shù)元素以及負數(shù)元素陣列與向量x的乘積結(jié)果bcpos與bcneg.

4) 將未分塊元素列表unblock_list中的元素都整合到向量b的中間結(jié)果中，得到b的最終結(jié)果.

4 基于異構(gòu)陣列的自選尾數(shù)壓縮系統(tǒng)實現(xiàn)

本工作將用軟件模擬實現(xiàn)所設(shè)計的憶阻MVM系統(tǒng).具體地，將用Python程序模擬實現(xiàn)硬件的存儲結(jié)構(gòu)、運算框架的構(gòu)造以及數(shù)據(jù)在系統(tǒng)中的計算.本節(jié)將對系統(tǒng)中關(guān)鍵的存儲結(jié)構(gòu)、系統(tǒng)處理流程的具體實現(xiàn)作出詳細的介紹.

4.1 硬件結(jié)構(gòu)的實現(xiàn)

本系統(tǒng)需要模擬的硬件主要包含憶阻運算陣列(存儲簇)以及乘加縮減樹結(jié)點.本系統(tǒng)將用Python的類(class)中的成員變量來模擬各個存儲結(jié)構(gòu)所包含的存儲項，若在運算過程中需要請求使用硬件資源，則只需要用軟件編程方法實例化對應(yīng)的類即可.具體地，本系統(tǒng)存儲結(jié)構(gòu)需要實現(xiàn)3種類.首先是有關(guān)陣列的2個類，分別是單個憶阻陣列和相同大小憶阻陣列組成的存儲簇，表1展示了其模擬存儲結(jié)構(gòu)(類中需要儲存的成員)；其次是乘加縮減樹結(jié)點硬件所對應(yīng)的類，表2展示了其模擬存儲結(jié)構(gòu).

Table 1 Data Storage Structure of a Storage Cluster and a Single Memristive Crossbar

Table 2 Storage Structure of Nodes in Multiply-and-Add Trees

4.2 系統(tǒng)處理流程的實現(xiàn)

在實現(xiàn)了模擬憶阻存儲結(jié)構(gòu)的基礎(chǔ)上，可以對憶阻陣列進行映射與計算.下面將介紹整個系統(tǒng)處理過程的實現(xiàn)要點.

4.2.1 整體實現(xiàn)流程

本節(jié)將從總體角度介紹如何用程序模擬數(shù)據(jù)在憶阻系統(tǒng)中的計算.圖12顯示了整個系統(tǒng)的實現(xiàn)流程，下面將對照該圖論述數(shù)據(jù)在系統(tǒng)中的計算以及系統(tǒng)的總體實現(xiàn)結(jié)構(gòu).

Fig. 12 Overall system implementation process圖12 系統(tǒng)整體實現(xiàn)流程

首先，從文件中讀取原始矩陣A(具體地，用Python的scipy.io庫讀取原始“.mat”文件，具體數(shù)據(jù)信息在5.1.2節(jié)中介紹).然后，進入預處理過程，該過程得到了經(jīng)過尾數(shù)壓縮與對齊位優(yōu)化后的異構(gòu)二進制矩陣塊，以及各矩陣塊對應(yīng)的信息.對每一個矩陣塊，都實例化對應(yīng)的存儲簇和乘加縮減樹結(jié)點的Python類，以達到用軟件方法模擬實現(xiàn)硬件資源的目的.預處理過程(包括分塊分片、尾數(shù)壓縮和對齊位優(yōu)化過程)都用Python源代碼實現(xiàn)，封裝成一個Python程序模塊，預處理具體數(shù)據(jù)時調(diào)用該模塊即可.最后進入運算過程，分為3個步驟：1)陣列運算，將向量x與憶阻陣列中0-1矩陣進行點乘，運算結(jié)果存儲在對應(yīng)變量中(用變量模擬圖8所示的緩沖輸出).2)進行乘加縮減樹運算，按照3.3.3節(jié)中的描述，編程實現(xiàn)每一個結(jié)點的移位和加法操作，同時維護結(jié)點正常的隊列長度，保證數(shù)據(jù)流在樹結(jié)構(gòu)中的正常路徑.3)輔助處理子模塊中運算的過程.主要包括部分和的求和、定點數(shù)轉(zhuǎn)浮點數(shù)、不同存儲簇結(jié)果與未分塊元素的整合，用變量來模擬圖8中的全局緩存區(qū).用Python源程序?qū)崿F(xiàn)整個運算過程，封裝成一個模塊，接收已完成數(shù)據(jù)映射的憶阻存儲結(jié)構(gòu)以及輸入向量x，返回矩陣向量乘法的結(jié)果.

4.2.2 預處理過程與陣列運算的實現(xiàn)

預處理過程主要包括分塊、分片、尾數(shù)壓縮和對齊位優(yōu)化、矩陣映射4個步驟，完成了映射之后，便可以加載向量元素進行陣列運算.

1) 分塊信息表block_list和未分塊信息表unblock_list的構(gòu)造與維護.在矩陣分塊過程中，按照塊由大到小循環(huán)遍歷整個矩陣，這個過程中需要記錄分塊元素信息，存在block_list中.表中元素為元組(x,y,size)，分別代表該分塊的左上角在原矩陣中的行數(shù)值、列數(shù)值以及塊的大小.如果遍歷到最小塊仍然滿足不了分塊閾值要求，或者遍歷完整個矩陣后有沒有處理到的邊角元素時，需要記錄其中的非0元素(即未分塊元素)，存在unblock_list中.該表中元素為元組(x,y,val)，分別代表該元素在原矩陣中的行數(shù)值、列數(shù)值以及對應(yīng)的浮點數(shù)值.

2) 矩陣塊處理映射過程.當block_list構(gòu)造完成之后，便可以進行相應(yīng)處理將每個矩陣塊映射到對應(yīng)的存儲簇上.首先，根據(jù)(x,y,size)提取當前需要處理的映射塊A′=A[x:x+size,y:y+size].對該矩陣塊A′進行分片、尾數(shù)壓縮與對齊位優(yōu)化，獲得一個提供最終映射信息的矩陣A″，其矩陣元素為(符號位，壓縮與對齊后的有效數(shù)二進制串，指數(shù))三元組.通過這個信息矩陣不難將矩陣塊按位片映射在表1所示的存儲簇結(jié)構(gòu)上，存儲關(guān)鍵數(shù)據(jù)信息.需注意，處理過程中要根據(jù)“符號位”的數(shù)值為1或0，選擇將其映射在負數(shù)陣列集合還是正數(shù)陣列集合中.

3) 陣列運算.完成映射后，需要根據(jù)當前存儲簇結(jié)構(gòu)存儲的y和size值，提取需要相乘的向量片段x′=x[y:y+size].然后，對x′進行分片與對齊，按照位片從高到低加時序地載到存儲簇的電壓端.對于每一次向量位片的加載，經(jīng)過流水乘加縮減樹的運算整合，輸出其與矩陣塊的乘積，即該向量位片對應(yīng)的部分和.

4.2.3 部分和求和的實現(xiàn)

1) 部分和求和計算.每加載一次x位片，乘加縮減數(shù)就輸出一次部分和.在此，采用動態(tài)策略進行部分和的求和運算，即每產(chǎn)生1次部分和，就進行1次加法運算，每一次求和計算：

(4)

其中，cpos(cneg)為當前進行部分和求和的存儲簇c對應(yīng)的正數(shù)(負數(shù))陣列集合，j代表加載第j個x位片，向量b則為該系統(tǒng)最后需要求得的結(jié)果，即Ax的結(jié)果.bcpos,j(bcneg,j)則代表當前存儲簇中正數(shù)(負數(shù))陣列集合在第j個x位片的作用下得到的部分和結(jié)果；bcpos,jsum(bcneg,jsum)代表前j個部分和結(jié)果相加得到的最終結(jié)果.

2) 部分和求和的提前終止策略.在上述過程中，隨著輸入的x位片指數(shù)降低，其在結(jié)果中的權(quán)值降低.事實上，在后續(xù)二進制轉(zhuǎn)浮點的操作中，結(jié)果中多余的低位片數(shù)值將會被截斷.本工作利用文獻[11]提到的方法，應(yīng)用部分和求和提前終止策略，設(shè)定結(jié)果最多保留位數(shù)為m(一般m=53)，則當結(jié)果bcpos,jsum(bcneg,jsum)的高m位穩(wěn)定不變時，求和過程結(jié)束.記結(jié)束時，當前存儲簇c中cpos和cneg對應(yīng)的部分和求和結(jié)果是bcpos和bcneg.

求和終止需要滿足3個條件：1)后續(xù)求和區(qū)域與高m位沒有重疊.2)高m位后面出現(xiàn)一個值為0的阻隔位，用于吸收后續(xù)求和元素產(chǎn)生的進位.3)值為0的阻隔位出現(xiàn)之后，需要再計算1次，若阻隔位沒有變成1，則求和才能終止.這一步是為了避免阻隔位翻轉(zhuǎn)讓后續(xù)計算向高m位產(chǎn)生進位.部分和求和終止時的形式如圖13所示：

Fig. 13 Partial sum early termination strategy圖13 部分和求和提前終止策略

4.2.4 定點數(shù)轉(zhuǎn)浮點數(shù)

expvalc=len_c+min_x+min_c-mantissa_bits,

(5)

4.2.5 不同存儲簇結(jié)果與未分塊元素的整合

4.2.2～4.2.4節(jié)所述的陣列計算、部分和的求和計算，均是針對同一個存儲簇中的運算而言的.要獲得Ax的結(jié)果，需要將不同的存儲簇中的結(jié)果整合起來.在這一過程中，需要全局地維護結(jié)果向量b.

1) 不同存儲簇結(jié)果的整合.遍歷所有的存儲簇進行運算，都會產(chǎn)生對應(yīng)的MVM運算結(jié)果bc，對所有存儲簇c整合到b中：

b[i+c.x]=b[i+c.x]+bc[i],
i=0,1,…,c.size-1.

(6)

2) 未分塊元素整合.計算的最后一步是將unblock_list中的元素整合到結(jié)果向量b中，對于unblock_list中的第i個元素(xi,yi,vali)進行計算(假設(shè)unblock_list的長度為len_u)：

b[xi]=b[xi]+x[yi] ×vali,
i=0,1,…,len_u-1.

(7)

4.3 預處理部分開銷分析

整個系統(tǒng)的開銷主要包括預處理開銷與運算開銷.其中，運算開銷主要包括模擬電路開銷(以運算陣列開銷為主)與數(shù)字電路開銷(以ADC轉(zhuǎn)換開銷為主).運算開銷的評估與分析將在第5節(jié)作詳細討論，本節(jié)主要對預處理部分的關(guān)鍵開銷進行定性分析與評估.

預處理部分的主要開銷來自于原始矩陣分塊開銷、確定塊內(nèi)指數(shù)范圍開銷、構(gòu)造葉子結(jié)點數(shù)可變的乘加縮減樹開銷這3個方面.其中，后2個開銷是原系統(tǒng)[11]中沒有涉及到的，即是由于應(yīng)用了尾數(shù)壓縮和對齊位優(yōu)化策略，系統(tǒng)在預處理部分多出的開銷.

1) 原始矩陣分塊開銷.假設(shè)一個矩陣中非0元素的個數(shù)為Nnz，按照本文提出的算法，將用4種不同大小的塊對原矩陣的非0元素密集部分進行捕捉，每個元素最多遍歷4次，即在最壞情況下需要進行4×Nnz次操作.已有的工作發(fā)現(xiàn)，在平均情況下，對于原始矩陣分塊的開銷大約為1.8×Nnz次操作[11].

2) 確定塊內(nèi)指數(shù)范圍開銷.為確定實際需要激活的運算對齊位個數(shù)，需求得每個矩陣塊元素中最大指數(shù)和最小指數(shù)的差值，這就要遍歷矩陣塊中所有非0元素.由于矩陣塊捕捉的都是原始矩陣中非0元素的密集區(qū)域，假設(shè)原矩陣中非0元素密度為K，邊長為l，則這一步驟總開銷(對整個原始矩陣)與Kl2成比例.

4.4 實現(xiàn)上的限制分析

本系統(tǒng)在實現(xiàn)時會遇到一些限制，本節(jié)將從硬件物理實現(xiàn)和實際測試負載應(yīng)用2個角度討論.

1) 硬件物理實現(xiàn).本文所提出系統(tǒng)的實現(xiàn)與評估都是通過編寫軟件程序模擬仿真實現(xiàn)的.這是由于目前的憶阻陣列制備工藝尚不成熟，不能夠支持如此大規(guī)模、擁有多陣列尺寸大小的憶阻系統(tǒng)的物理實現(xiàn).但使用軟件實現(xiàn)的模擬系統(tǒng)在評估測試中也有其獨特的優(yōu)勢，可以靈活調(diào)節(jié)軟件參數(shù)來改變分塊大小、壓縮位數(shù)等重要參數(shù)，快速得到各種情況下的理論結(jié)果.

2) 實際測試負載應(yīng)用.尚未找到一種方法，能對所有的測試負載都找到效果可觀的尾數(shù)壓縮和對齊位優(yōu)化方案.實際應(yīng)用該系統(tǒng)時，可能會由于部分參數(shù)選擇不當(如尾數(shù)壓縮位數(shù)、對齊位閾值)，達不到需要的優(yōu)化效果或精度要求.

5 評估與分析

本節(jié)將在線性代數(shù)求解框架中集成所實現(xiàn)的系統(tǒng)，對不同求解精度下的各種能耗節(jié)省比例進行評估.

5.1 評測數(shù)據(jù)與評估方法

下面將對模擬MVM內(nèi)核集成到算法框架中的方法、評估所用的數(shù)據(jù)以及評估方式作詳細論述.

5.1.1 MVM內(nèi)核的集成方法

1) 算法求解框架.常用的高性能線性代數(shù)迭代求解算法有共軛梯度法(conjugate gradient method, CG)[19]、雙共軛梯度穩(wěn)定法(biconjugate gradient stabilized method, BICGSTAB)[20].在本課題中，將使用Python的scipy.sparse.linalg庫中已經(jīng)實現(xiàn)的CG/BICGSTAB算法框架來進行測試評估.這種算法框架支持以線性運算符(linear operator)的形式將模擬MVM內(nèi)核應(yīng)用于迭代求解，算法框架(以CG為例)：

(8)

在輸出端元素中，x為線性系統(tǒng)的解，flag為收斂標志(0代表收斂).在輸入端參數(shù)中，A代表原始矩陣或者線性運算符，b代表線性方程組的右手向量，x0代表解向量的初始值，默認為零向量；tol代表解的容差，即|b-Ax|

2) 運用線性運算符集成MVM內(nèi)核.在式(8)中，輸入?yún)?shù)A可以是線性運算符，其作用是包裝一個函數(shù)，函數(shù)輸入為x，返回Ax的結(jié)果給算法求解框架.Python的scipy.sparse.linalg庫提供線性運算符的包裝方案，其框架：

A=LinearOperator(A.shape,matvec,rmatvec),

(9)

其中，matvec和rmatvec是2個關(guān)鍵函數(shù)，它們的輸入都是x，返回分別為Ax和ATx.為實現(xiàn)MVM內(nèi)核的集成，需編寫這2個函數(shù)，在函數(shù)體中調(diào)用實現(xiàn)的MVM模擬系統(tǒng)，通過輸入的向量x和固定的矩陣A，計算得到對應(yīng)的Ax或ATx，作為函數(shù)的返回.

3) 預條件子的使用.在式(8)中，輸入?yún)?shù)M代表預條件子矩陣.預條件子是在解決大型稀疏矩陣求解問題時的一種常用優(yōu)化手段[21-22].關(guān)于預條件子矩陣，通俗的解釋就是需要找到一個和稀疏矩陣A相似的矩陣M，對于對原等式Ax=b作出M-1Ax=M-1b的變換.其中M-1A和M-1b都可以直接算出，而且由于M和A為相似的矩陣，所以M-1A的值接近于單位矩陣，極大程度上減少了求解的步數(shù).有一種通用的預條件子方法為0級不完全LU分解(incomplete LU factorization with level-0 fill in, ILU(0))，該方法利用單位下三角矩陣L和上三角矩陣U構(gòu)造M=LU，使得M接近于A，并在矩陣中執(zhí)行0級填充，保持矩陣原有的稀疏模式.在Python的scipy.sparse.linalg庫中，有spilu這樣一個求解類，可以用于構(gòu)造ILU(0)預條件子，并集成到算法框架中.在后續(xù)測試中，存內(nèi)硬件求解器都將集成ILU(0)預條件子.

5.1.2 評測數(shù)據(jù)集

評測將采用SuiteSparse矩陣集合[2]中的矩陣作為數(shù)據(jù)集.SuiteSparse矩陣集合是科學與工程領(lǐng)域中的實際應(yīng)用產(chǎn)生的稀疏矩陣集合，該集合被廣泛地用于稀疏矩陣算法的開發(fā)和性能評估.

根據(jù)需要，選用了大小不一的6個稀疏矩陣，作為評測數(shù)據(jù).其中，部分數(shù)據(jù)文件只給出了矩陣A的值但是沒有給出右手向量b的值.對于這部分數(shù)據(jù)，在進行CG/BICGSTAB算法求解時參照文獻[21]的做法，構(gòu)造全是1的右手向量b進行測試.表3顯示了6個矩陣的名稱、矩陣行數(shù)(rows)、非0元素數(shù)量(Nnz)、平均每行非0元素個數(shù)(Nnz/rows)等基本信息.

為了后續(xù)評估不同求解精度下的能量消耗比例，需要一些有關(guān)于憶阻器件的參數(shù)來輔助進行能耗評定.表4列舉出了一些真實的憶阻器件參數(shù)[9,23].

Table 3 Evaluated Matrices Dataset表3 評測矩陣數(shù)據(jù)集

Table 4 Memristive Device Parameters表4 憶阻器件參數(shù)

5.1.3 評估方法

在下面的評估中，主要考慮只進行對齊位優(yōu)化、尾數(shù)保留35 b、尾數(shù)保留25 b、尾數(shù)保留15 b這4種優(yōu)化策略(3種尾數(shù)保留策略也是在對齊位優(yōu)化的基礎(chǔ)上進行).在憶阻MVM運算中，這4種策略所激活的陣列個數(shù)依次減少，會導致求解精度下降，但是會獲得極大的能耗節(jié)約.

為了評估所設(shè)計的自選尾數(shù)壓縮機制和動態(tài)對齊位優(yōu)化機制的有效性和效果，將在CG和BICGSTAB這2種算法下評估4種策略的3個指標：1)相對于軟件求解基線的相對殘差及對應(yīng)的迭代次數(shù).2)運算陣列能耗相對于能耗基線的優(yōu)化程度.3)ADC能耗相對于能耗基線的優(yōu)化程度.其中，能耗基線指在原系統(tǒng)[11]對應(yīng)的條件下(尾數(shù)53 b、對齊位64 b固定不變)測得的對應(yīng)種類的能耗(在本文的測試中，對齊位閾值也選擇為64 b).

通過評估這3個指標，可以分析出隨著激活運算陣列的減少，計算精度的下降(通過計算對于軟件求解基線的相對殘差來體現(xiàn)求解的精確程度)，運算陣列和ADC能耗的優(yōu)化程度.驗證所提出的系統(tǒng)是否可以在為高精度應(yīng)用提出無損浮點計算功能的同時，又能有效降低較低精度應(yīng)用的能耗開銷，以體現(xiàn)系統(tǒng)實現(xiàn)方案的有效性和良好效果.

5.2 不同策略下的求解精度與迭代次數(shù)

由于只進行對齊位優(yōu)化的情況下并不會影響MVM的計算精度，即此策略下的解向量與迭代次數(shù)與軟件求解完全一致，所以不再對其進行單獨測定.本節(jié)主要關(guān)注5.1節(jié)提到的3種尾數(shù)壓縮策略.

5.2.1 評估標準

用殘差來描述不同策略下求解精度相對于軟件求解基線的差距，求解殘差：

(10)

5.2.2 評估結(jié)果與分析

本節(jié)展示在各種尾數(shù)壓縮策略下，求解結(jié)果相對于軟件基線的相對殘差以及算法求解的迭代次數(shù).

1) 各種尾數(shù)壓縮策略下解的相對殘差.圖14分別顯示了在BICGSTAB(圖14(a))、CG(圖14(b))測試算法下，3種尾數(shù)壓縮策略相對于軟件基準的相對殘差.總體來說，在對6個數(shù)據(jù)集測試結(jié)果進行對數(shù)平均處理后，3種策略的相對殘差分別在10-9，10-7，10-3這一數(shù)量級.具體地，對于同一種尾數(shù)壓縮策略，2種算法解向量的相對殘差是幾乎相同的，可以認為壓縮策略影響的是映射矩陣本身的性質(zhì)，不會因為算法框架的改變而發(fā)生太大的結(jié)果變化.對于每一種算法的每個測試矩陣而言，隨著保留位數(shù)的減少，相對殘差都會成數(shù)量級地增大.

2) 各種尾數(shù)壓縮策略下算法的迭代次數(shù).圖15顯示了在2種測試算法下各策略對應(yīng)的算法迭代次數(shù).總體來說，進行適當范圍內(nèi)的尾數(shù)壓縮不會對迭代次數(shù)產(chǎn)生大的影響.但是，如果壓縮的位數(shù)太多，保留位數(shù)太少，可能會增大迭代次數(shù)(如圖15中保留15 b尾數(shù)時的測試矩陣msc00726的測試結(jié)果).為了進一步驗證該想法，分別繼續(xù)測試了測試矩陣msc00726在只保留15 b，12 b，9 b，6 b這4種尾數(shù)位情況下2種算法的迭代次數(shù).結(jié)果顯示，BISGTAB算法的迭代次數(shù)對應(yīng)分別為2，2，3，5，CG算法的迭代次數(shù)對應(yīng)分別為3，3，4，9.可以看出，如果尾數(shù)保留位數(shù)過少，會顯著增加求解迭代次數(shù)，使得迭代求解過程需要更多的時間和能耗.故在實際應(yīng)用中，需要謹慎選擇尾數(shù)壓縮位數(shù)，盡量不選擇過少的保留位數(shù)(如15 b以下)，以防止尾數(shù)壓縮策略給求解帶來的負面影響.

Fig. 15 Algorithm iteration times under various compaction strategies圖15 各種壓縮策略下的算法迭代次數(shù)

5.3 不同優(yōu)化策略下的運算陣列能耗評估

運算陣列能耗是模擬電路中最大的能耗來源.進行陣列運算時，電壓加載到電阻上，會產(chǎn)生瞬時功率.為了使得采樣—保持數(shù)組能夠獲得正確的電流值，通常需要使得電流維持一段時間，所以會產(chǎn)生能量消耗.本節(jié)將展示在4種優(yōu)化策略下，運算陣列能耗相對于能耗基線的減少比例.能耗基線的值是在已有未優(yōu)化系統(tǒng)[11]的條件下評估得到的，其尾數(shù)為固定的53 b，對齊位為固定的64 b.

5.3.1 評估指標

參數(shù)如表4所示，每一個存儲單元的瞬時功率：

(11)

采樣—保持數(shù)組需要電流維持的時間主要取決于其RC時間常數(shù)，它與分辨率，即lbN是成比例的，其中N是當前憶阻存儲陣列的大小(行數(shù))，即一個存儲單元的計算能耗與P×lbN是成比例的[11].

為了獲得進行一次Ax運算所對應(yīng)的運算陣列能耗Ecrossbar，需要依次進行3個具體步驟計算：

1) 對于加載1次電壓(1個向量位片)而言，需要對所有陣列中的所有存儲單元的能耗值進行求和，得到憶阻存儲結(jié)構(gòu)的整體能耗;

2) 對于每一次向量位片的加載，將上述整體能耗進行累加，得到Ecrossbar;

3) 為了求解迭代過程中的運算陣列總能耗，需要對所有Ax步驟的Ecrossbar求和，得到一個與迭代陣列消耗總能量成比例的值Ecrossbar_tol.

5.3.2 評估結(jié)果與分析

評估求解過程中，各策略Ecrossbar_tol的相對比例.

Fig. 16 Computation crossbar energy ratio under various strategies圖16 各種策略下的運算陣列能耗比例

1) 各種壓縮/優(yōu)化策略下的運算陣列能耗比例.以測試矩陣mesh3em5在BICGSTAB算法下的測試得到的Ecrossbar_tol為單位能量，畫出在2種算法下，所有矩陣在各種策略下的能耗值，如圖16所示.首先，對于同一個算法同一個測試矩陣而言，隨著每個存儲簇中激活運算的陣列越來越少(從只進行對齊位優(yōu)化，到尾數(shù)保留位數(shù)分別只有35 b，25 b，15 b)，運算陣列消耗的能量也逐漸變少，符合預期.同時，可以觀察到，矩陣中非0元素越多，消耗的陣列能量明顯越多.這是因為其非0元素多，導致其映射電阻中Ron較多，導致能耗較大.并且注意到，消耗的運算陣列能量與原始矩陣大小沒有直接關(guān)系，比如在6個測試集中，最大的矩陣(行數(shù)最多的)為測試矩陣Chem97ZtZ，但是其中的非0元素沒有測試矩陣bcsstk27多，所以消耗能量明顯比測試矩陣bcsstk27小很多.對于同一個測試矩陣而言，由于CG算法的迭代次數(shù)多于BICGSTAB，所以消耗的運算陣列能量也成比例的增加.

2) 各種壓縮/優(yōu)化策略下的運算陣列平均能耗優(yōu)化.表5展示了對于所有6個測試矩陣平均值而言，使用每一種策略時運算陣列能耗相對于能耗基線減少的比例.可以發(fā)現(xiàn)，BICGSTAB和CG這2種方法的能量優(yōu)化比例是幾乎相同的，并且隨著運算時激活陣列的減少，能耗優(yōu)化比例會越來越大.具體地，在相對于軟件方法有0～10-3范圍的求解殘差時，2種算法下的平均運算陣列能耗都相對已有的未優(yōu)化系統(tǒng)減少了5%～65%.

Table 5 Average Computation Crossbar Energy Optimization Ratio

5.4 不同優(yōu)化策略下的ADC能耗評估

ADC是將采集-保持數(shù)組中的電信號轉(zhuǎn)化為數(shù)字信號的部件，只有經(jīng)過這種轉(zhuǎn)換，乘加縮減樹才能對陣列中的計算結(jié)果求和.ADC能耗是外圍數(shù)字電路能耗的主要來源.本節(jié)中將展示4種優(yōu)化策略下的ADC的能耗相對于能耗基線(已有未優(yōu)化系統(tǒng)[11]條件下評估得到)的減少.

5.4.1 評估指標

假設(shè)一個陣列的大小是N，則其所需要的ADC分辨率為lbN.因為ADC的平均功耗隨分辨率呈指數(shù)增長[24-26]，故其功耗和N成比例.又因為ADC的轉(zhuǎn)換時間與分辨率lbN成正比，所以列電流的轉(zhuǎn)換能耗與N×lbN成正比[11].

假設(shè)某個存儲簇中有M個存儲陣列.即對于每一個存儲簇，每一個單獨的向量x位片要進行M×N次列電流的計算，則存儲簇進行一次計算ADC的總能耗與M×N×N×lbN成正比.

為了獲得一次Ax的完整ADC能耗EADC.需要在算法執(zhí)行過程中對所有存儲簇，所有的向量x位片的ADC能量求和，得到最終EADC.注意EADC的值只是與真實的ADC能耗成比例，而非真實ADC能耗值，其只會被用來計算相對于能耗基線的能量節(jié)省比例.

最后，為了計算求解迭代過程中的ADC總能耗，需要對所有Ax計算產(chǎn)生的EADC求和，得到求解整個算法ADC消耗的總能量值EADC_tol，該值與迭代求解過程中真實的ADC能量消耗成比例.

5.4.2 評估結(jié)果與分析

評估求解過程中，各策略下EADC_tol的相對比例.

1) 各種壓縮/優(yōu)化策略下的ADC能耗比例.以測試矩陣mesh3em5在BICGSTAB算法下的測試得到的EADC_tol為單位能量，畫出所有矩陣在各種算法、各種執(zhí)行策略下的能耗比例.圖17顯示出一些與運算陣列能耗優(yōu)化(圖16)相同的性質(zhì).比如，矩陣塊的非0元素越多，ADC的能耗明顯越多，因為需要更多的陣列進行運算，對應(yīng)的ADC轉(zhuǎn)換過程也就越多.再比如由于CG算法的迭代次數(shù)明顯比BICGSTAB算法要多，所以其ADC能耗也成比例的多.

對同一個算法同一個測試矩陣而言，隨著每個存儲簇中激活運算的陣列越來越少，ADC能耗基本上也逐漸變少，符合預期.但是需要注意，對于測試矩陣msc00726和測試矩陣nos5而言，應(yīng)用對齊位優(yōu)化策略相對于能耗基線并沒有一個明顯的ADC能耗下降，這是因為其矩陣塊內(nèi)指數(shù)范圍太大，超過了設(shè)定的閾值64，對齊位優(yōu)化效果不明顯.

對于大多數(shù)矩陣而言，ADC能耗隨著保留的尾數(shù)位數(shù)的減少而降低.但是對于測試矩陣msc00726而言，在保留15 b尾數(shù)的時候，其ADC能耗反常地比保留25 b時要多，這是由于尾數(shù)截斷過多使得算法迭代的次數(shù)增加，導致要用到更多計算步驟來使得算法收斂.所以須注意，由于壓縮位數(shù)過多導致迭代次數(shù)增加，最終導致計算能耗過多的情況也會發(fā)生，對于ADC能耗來說可能更明顯，需要謹慎選擇壓縮位數(shù).

Fig. 17 ADC energy ratio under various strategies圖17 各種策略下的ADC能耗比例

2) 各種壓縮/優(yōu)化策略下的ADC平均能耗優(yōu)化.表6展示了對于所有6個測試矩陣平均值而言，每種執(zhí)行策略的ADC能耗相對于能耗基線減少的比例.和運算陣列能耗一樣，各種策略在應(yīng)用于BICGSTAB和CG這2種算法的情況下，對ADC能耗減少的比例幾乎是一致的，且隨著運算激活陣列數(shù)的減少，能量優(yōu)化比例逐漸增大.具體地，在相對于軟件方法有0～10-3范圍的求解殘差時，2種算法下的平均ADC能耗都相對于已有的未優(yōu)化系統(tǒng)減少了30%～55%.

Table 6 Average ADC Energy Consumption Optimization Ratio

與運算陣列能耗優(yōu)化的情況不同的是，ADC能耗在只應(yīng)用對齊位優(yōu)化時，其能耗優(yōu)化程度約為30%，它大于同樣策略下的運算陣列優(yōu)化的數(shù)值(約5%).但是，隨著存儲簇中陣列激活數(shù)量(保留尾數(shù)位數(shù))的逐漸減少，ADC能量優(yōu)化的效果沒有運算陣列能量明顯.運算陣列的能耗在保留15 b的時候優(yōu)化比例在65%左右，而ADC能耗在尾數(shù)保留位數(shù)15 b時優(yōu)化比例只有55%左右.

6 結(jié)束語

在科學與工程領(lǐng)域中，浮點數(shù)線性系統(tǒng)的求解問題十分普遍.現(xiàn)有的研究提出將IEEE-754雙精度浮點數(shù)部署在憶阻陣列上，進行高能效的存內(nèi)計算.然而，還沒有研究提出為不同求解精度的浮點線性系統(tǒng)優(yōu)化運算能耗.為解決這個問題，本文創(chuàng)新性地提出了一種具有尾數(shù)壓縮與對齊位優(yōu)化策略的浮點MVM計算系統(tǒng)，并針對該策略的特征提出了一種葉子結(jié)點數(shù)可變的流水乘加縮減樹.評估結(jié)果表明：本系統(tǒng)既可進行高精度浮點數(shù)無損運算，又可對低精度求解應(yīng)用大幅度優(yōu)化運算陣列與ADC的能耗.

作者貢獻聲明：丁文隆提出了本文工作的思路，設(shè)計完成了實驗，并撰寫論文；汪承寧協(xié)助改進系統(tǒng)框架，在論文寫作過程中提供理論基礎(chǔ)的幫助；童薇對整個工作提出指導意見并修改論文.