基于Epiphany 計算并行矩陣乘法的研究

曾弘揚

（北京工業大學 軟件學院， 北京 100124）

對于工程科學中的許多計算問題中，數值計算問題都是最基本的內容。 矩陣乘法作為基本的線性代數運算，被廣泛地應用于工程學領域[1]。 矩陣乘積在實際應用中是經常要用到的，例如，在解決有限元素的方法中，模型元素之間的關系被表示為狀態矩陣元素。 對狀態的改變，被表示成輸入一個矩陣或向量，通過進行矩陣乘法來計算出新的狀態。 所以如何高效的計算矩陣乘法是十分重要的。

1 矩陣乘法的串行實現

矩陣相乘串行實現的i－j－k 形式， 其中最內層運算是點積，數據用到A 的行和B 的列。

for（i=0；i＜M；i++）

for（j=0；j＜N；j++）

for（k=0；k＜P；k++）

C[i][j]+=A[i][k]*B[K][j]；

矩陣相乘的前提條件，前一個矩陣列數，須與后一個矩陣行數相等。 矩陣A 為M 行P 列，矩陣B 為P 行N 列，則矩陣乘法，繼承前一個矩陣的行數M，和后一個矩陣的列數N，相乘矩陣中每一項，都要經過P 次加和乘。 算法分析：該算法需要進行n3個乘法運算和n3個加法運算， 順序代碼的時間復雜度為O（n3）。

2 串行計算大矩陣乘法的問題

2.1 串行計算不適應于大矩陣乘法

許多先進的計算機都配有高效的用于解矩陣乘法的的串行程序庫，比如常用的SGEMM（單精度普通矩陣乘法）函數就是用于實現矩陣乘法運算一個常規形式的標準API，但它不適用于很大的矩陣。 因此，為了在并行計算環境下實現矩陣乘積，研究并行矩陣乘法是非常必要的[2]。

2.2 并行計算時數據在多核間的傳輸速率慢

在多核處理器上，當操作的矩陣足夠大時，數據就必須從外部內存中取出，例如系統的DRAM。 而Epiphany 芯片的外部內存DRAM，通過eMesh 網絡在STARTIX FPGA 上來連接。 這個擴展的eMesh 時鐘頻率是50 MHz，相較于芯片的時鐘來說是很慢。 因此，數據在DRAM 的傳輸非常緩慢。 通過eMesh 連接的Epiphany 能夠在每個時鐘周期讀寫一個字節的數據。 在一個負載均衡的系統上，400 MHz 的Epiphany E16G3 芯片， 從芯片到DRAM 的理論寫速率是400 MB/s，但實際測量到的寫速率是82 MB/s[3]。

可見，為了生成優化的代碼，應盡可能減少對速度相對緩慢的DRAM 的訪問。因此，在設計算法時，我們必須試圖盡可能多的去重復利用這些讀取出來的數據。

3 矩陣向量乘法的并行計算

矩陣與向量相乘是線性方程組迭代求解的核心問題，其并行計算的效率直接影響迭代求解的整體效率。n 階矩陣A，與n 維向量x=[x1x2… xn]T，并行計算A 與x 的乘積y=[y1y2…yn]T。

并行矩陣向量乘法，對矩陣采用一維塊狀劃分（帶狀劃分），即把矩陣按整行（或整列）劃分成組，將每組的數據指派給一個處理器存儲。 劃分的這些行列可以是連續的（連續帶狀劃分），也可以是等距相間的（循環帶狀劃分）。 這里，我們對矩陣A 進行按列連續劃分， 將矩陣A 按行連續劃分成p塊，則每塊所占行數為n/p 行（其中n 能被p 整除）。同時再將每個行塊按列也相應地劃分成p 個子塊， 則第k 個行塊Ak又可 進一步劃 分為[Ak，0Ak，1… Ak，p－1]。 則Ak，j可以表 示如下：

圖1 矩陣的二維網格劃分Fig. 1 The matrix of the two-dimensional grid

圖1 中劃分矩陣A 的下標k 是按行連續劃分的下標，下標j 是與處理器個數相對應的列上劃分， 它們都與處理器個數p 相關，因此范圍都是[0，…，p－1]。

n 階矩陣A 的子塊分成了n/p 階， 同理后面所要乘的向量x，與所得新向量y 的子塊，也都要分成行為n/p 的子塊，因此也都要按 行來進行 劃分成：x=[x0x1…xp－1]T，y=[y0y1… yp－1]T，它們對應成的第k 塊分別是：xk=[xk*n/p+1xk*n/p+2… xk*n/p+n/p]T，yk=[yk*n/p+1yk*n/p+2… yk*n/p+n/p]T。 劃分到處理器上，在處理器k 上進行的第k 塊乘積yk的計算公式為：

列下標j 按[0，1， …，p－1]移動時，下標（k+j） mod xk的模依次為[k，k+1，…，k－1]，即先做本地處理器上的Ak，kxk，再做Ak，k+1xk+1，…，就這樣不斷地循環向下移動，直至做到Ak，k+1xk－1完成整個循環， 即按順序逐次遍歷完整個x 向量上的每一個子塊為止。

矩陣A 的子塊行坐標始終為k 不變，即按上式進行計算時，所用到的矩陣A 中的子塊，就存儲在當前運行的處理器上。而向量x 的子向量下標為（k+j） mod p，會隨著j 的變化而變化， 即在計算的過程中需要用到整個x 向量上的所有子向量，因此可以通過在每次計算完成后，就將處理器中所存儲的x 子向量，循環向上移動到上一列的處理器中。

圖2 并行計算矩陣向量乘積Fig. 2 Parallel computing matrix vector product

如圖2 所示， 將各處理器中所存儲的矩陣A 中的子塊Ak，與處理器中存儲的向量x 的當前子塊xk相乘，進行p 次，每次x 向量中的子向量都向上循環移動一個位置， 對同一處理器中p 次的乘積進行累加求和， 即可得出矩陣向量的乘積向量y。

4 矩陣乘法的并行Cannon 算法實現

根據矩陣A 與矩陣B 在處理器中的不同存儲方式可得到多種并行計算矩陣乘積的劃分方法。 比如：行列劃分算法，將矩陣A 按行連續劃分成p 個一維塊狀子矩陣（行塊），將矩陣B 按列連續劃分成p 個一維塊狀子矩陣（列塊）。 但無論是按行列、行行、列列，還是按列行來對矩陣進行劃分，其本質都是基于矩陣的一維塊狀劃分。而Cannon 算法是基于矩陣的二維塊狀劃分的。 其特點是矩陣A 中的子塊在指定行中循環移動，矩陣B 中的子塊在指定列中循環移動，對每個處理器來說，計算量都是相同的，具有很好的負載均衡。 當p＞=4 時，Cannon 算法具有優越性[4]。

4.1 二維連續塊狀劃分

Cannon 算法中并行計算矩陣的乘法，與并行計算矩陣向量的乘法原理類似。但不同于矩陣向量的乘法中，只讓一邊的矩陣的各子矩陣循環移動。Cannon 算法是通過矩陣A 中的各行塊，與矩陣B 中的各列塊同時進行循環移位，來完成對矩陣C 中子塊的計算。

一個二維網絡由p*p 個處理器組成，將矩陣A、B、C 各劃分成p*p 塊，按二維連續塊狀進行劃分，如上面的矩陣A 圖的形式。 則處理器Pi，j上面的子矩陣Ci，j的計算公式如下：

其中，i 是行下標，j 是列下標，范圍都是[0，p－1]。 而k 是點積次數，即在m*l*n 中所乘的相同維度的l，因此k 的范圍同其他維度一樣也是[0，p－1]。

k 從0 變化到p－1 時，（i+j+k） mod p 也從0 取到p－1，k每增加1， 相應的行列坐標也循環遞增1。 因此，i，（i+j+k）mod p 是i 行的所有塊，（i+j+k） mod p，j 是j 列的所有塊。 所以，式子中的Ci，j就是A 的i 行上與B 的j 列上對應的所有塊乘積之和。

4.2 對 準

相乘的關鍵是相乘的兩個元素下標要滿足對準的要求。

在存儲數據時， 處理器Pi，j上存儲的數據應當是矩陣A的子塊Ai，j，與矩陣B 的子塊Bi，j。 而在計算開始時，從初始狀態k=0 開 始，在處理器Pi，j上處 理 的 是Ai，（i+j）modp* B（i+j）modp，j的乘積。 可見，需要處理的兩個子塊均不在當前的處理器上。 因此，在計算k=0 狀態之前，須要先進行對準操作，將處理器Pi，（i+j）modp上的子塊Ai，（i+j）modp及處理器P（i+j）modp，j上的子塊B（i+j）modp，j都先傳遞到Pi，j上后，才能開始計算。 即對p*p 二維網格上的每一個處理器Pi，j來說，都要進行的對準操作是，先將其上A的子塊向左循環移動i 個位置，傳遞給位于i 行（i+j） mod p列上的處理器； 再將其上B 的子塊向上循環移動j 個位置，傳遞給位于（i+j）mod p 行j 列上的處理器。

對準完成后，就可以從k=0 起始，從0 到p-1 按步進行計算了。

4.3 反對準

當計算p 次，即第p-1 次計算結束后，在Pi，j上參與計算的子塊實際是Ai，（i+j+k）modp和B（i+j+k）modp，j。要想將之 前執行的 對準進行還原，需要進行反對準，使Pi，j上存儲的操作子塊恢復成Ai，j與Bi，j。 因此，要將多核 上的每一個處理 器Pi，j上存儲 的A矩陣子塊向右循環移動i 個位置，B 矩陣子塊向下循環移動j個位置。

5 基于Epiphany 的矩陣乘法的Cannon 算法

為了實現一個任意大小矩陣乘法的高效的多核并行運算，需要使用單線程的C 代碼來構建運算中的子塊。其中，內存的分配被分為兩個主要部分。 第一部分，是在芯片上每一個單核上的內存。 里面存放的是子塊點積計算的結果。 第二部分，是芯片外的DRAM。 在這里，存放的是要操作的A、B、C矩陣。 以及，應用于系統層的主從設備間通信的Mailbox，和用于控制主從設備交互的結構體。

5.1 程序在核間的存儲

圖3 單核中的內存分配Fig. 3 The mononuclear memory allocation

Epiphany 芯片的每個核都具有32 kb 的SRAM， 將其分配成4 個8 kb 的bank。每個核都可并發訪問多個bank，且不會帶來性能損失。 bank 0 用來存儲程序代碼。 bank 1 用來存放操作矩陣A 中子塊的兩個緩存。 同理，bank 2 用來存放操作矩陣B 中子塊的兩個緩存。 將bank 3 中的內存等分為兩部分。 bank 3 中低位的部分，用于存放矩陣C 中子塊的計算結果，因為計算結束后，也不需寫回計算后的結果，因此不需要用兩個緩存來存放矩陣C 中的子塊。 將bank 3 中高位的部分， 用于存放對單核進行控制的結構體， 和核間交互的mailbox，以及運行bank 0 中代碼時，程序所需要的棧。內存分塊情形如圖3 所示。

5.2 子塊在多核間的傳遞

為了討論方便，先給出一些記號的約定：

假設有p 個處理器，每個處理器上運行一個進程，則Pj表示第j 個處理器或進程，Pmyid表示當前運行程序的處理器或進程。 算法都是在Pmyid上運算的，處理器排序均以0 為起始，因此分成的塊也從0 起始，矩陣的行列均從1 開始，如n階矩陣A 為[a1，1… an，n]。

send 和recv 分別表示在Pmyid中的發送和接收操作。send（x，j）表示將運行當前進程的Pmyid處理器中的數據x 發送給Pj處理器。recv（x，j）表示將Pj處理器中的數據x 接收到運行當前進程的Pmyid處理器中。

因此在子塊 移 動 的實現為順序 為Send （Ai，j，Pi，（p+i-j）modp）；Recv（Ai，j，P（i+j）modp）；即 先 把 本 處 理 器 中 計 算 出 的 數 據 循 環 向前移動到第i 個位置處理器的內存PING 中， 再接收循環向后第i 個位置處理器的數據放到本地內存的PONG 中[5]。 在Epiphany 芯片中的實際操作如圖4 所示。

圖4 子塊在多核間移動Fig. 4 Sub-block mobile between multicore

6 實驗與結果分析

通過計算兩個512*512 大小的矩陣乘積來測試該算法。分別使用在雙核AMD 處理器的臺式機， 和單核的Epiphany芯片上，使用串行算法；以及在16 核的Epiphany 芯片在上運行優化后的并行程序， 來計算相同的512*512 矩陣乘法，每組進行10 次實驗得到數據如下：

通過表格中的數據可以發現：

與在標準臺式電腦上得到的結果相比較，我們看到基于Epiphany 系統的巨大優勢——在一個負載平衡的系統中和采用更低的時鐘頻率，運算速度卻可以提高5 倍以上。

表1 計算矩陣乘法耗時時間Tab. 1 Calculate matrix multiplication takes time

進一步研究， 在基于16 核的Epiphany 芯片在上運行的并行程序的各階段的時間花費，測量時間如下：

表2 并行算法各階段的耗時Tab. 2 Every stage of the parallel algorithm’s time-consuming

通過表格中的數據可以發現：

1）總時間并不等于各項時間的和，在并行運算在多核上的性能會更低些。

2） 只有大約10%的時間是用來做真正的矩陣乘法運算的，其余時間都用來進行對DRAM 的讀寫。

7 結束語

本文在分析矩陣乘法的串行算法的基礎上，找出了影響其計算效率的問題，給出了通過并行計算來實現可擴展的大矩陣乘法的Cannon 算法， 以此來提高大矩陣乘法的計算效率。Cannon 算法在多核并行計算矩陣乘法中具有廣泛的應用性[6]。

[1] 李曉梅，吳建平. 數值并行算法與軟件[M]. 北京：科學出版社，2007.

[2] 廖曉鐘，賴汝. 科學與工程計算[M]. 北京：電子工業出版社，2003.

[3] Yaniv Sapir. Scalable Parallel Multiplication of Big Matrices[EB/OL].（2013-06）[2015-03].http://www.adapteva.com/wpcontent/uploads/2012/10/Scalable -Parallel -Multiplication -of-Big-Matrices.pdf.

[4] 姜弘道，張健飛. 工程科學中的高性能計算[D]. 北京:科學出版社，2013.

[5] Adapteva， Inc. Approaching Peak Theoretical Performance with Standard C [EB/OL].（2013-06）[2015-03].http://www.adapteva.com/white-papers/approaching-peak-theoreticalperformance-with-standard-c/.

[6] 李小洲， 李慶華. 矩陣相乘cannon并行算法在工作站機群上的實現[J]. 計算機工程，2002（6）:102-130.