曾弘揚(yáng)

（北京工業(yè)大學(xué) 軟件學(xué)院， 北京 100124）

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

1 矩陣乘法的串行實現(xiàn)

矩陣相乘串行實現(xiàn)的i－j－k 形式， 其中最內(nèi)層運(yùn)算是點積，數(shù)據(jù)用到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]；

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

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

2.1 串行計算不適應(yīng)于大矩陣乘法

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

2.2 并行計算時數(shù)據(jù)在多核間的傳輸速率慢

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

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

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

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

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

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

圖1 中劃分矩陣A 的下標(biāo)k 是按行連續(xù)劃分的下標(biāo)，下標(biāo)j 是與處理器個數(shù)相對應(yīng)的列上劃分， 它們都與處理器個數(shù)p 相關(guān)，因此范圍都是[0，…，p－1]。

n 階矩陣A 的子塊分成了n/p 階， 同理后面所要乘的向量x，與所得新向量y 的子塊，也都要分成行為n/p 的子塊，因此也都要按 行來進(jìn)行 劃分成：x=[x0x1…xp－1]T，y=[y0y1… yp－1]T，它們對應(yīng)成的第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 上進(jìn)行的第k 塊乘積yk的計算公式為：

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

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

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

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

4 矩陣乘法的并行Cannon 算法實現(xiàn)

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

4.1 二維連續(xù)塊狀劃分

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

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

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

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

4.2 對 準(zhǔn)

相乘的關(guān)鍵是相乘的兩個元素下標(biāo)要滿足對準(zhǔn)的要求。

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

對準(zhǔn)完成后，就可以從k=0 起始，從0 到p-1 按步進(jìn)行計算了。

4.3 反對準(zhǔn)

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

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

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

5.1 程序在核間的存儲

圖3 單核中的內(nèi)存分配Fig. 3 The mononuclear memory allocation

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

5.2 子塊在多核間的傳遞

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

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

send 和recv 分別表示在Pmyid中的發(fā)送和接收操作。send（x，j）表示將運(yùn)行當(dāng)前進(jìn)程的Pmyid處理器中的數(shù)據(jù)x 發(fā)送給Pj處理器。recv（x，j）表示將Pj處理器中的數(shù)據(jù)x 接收到運(yùn)行當(dāng)前進(jìn)程的Pmyid處理器中。

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

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

6 實驗與結(jié)果分析

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

通過表格中的數(shù)據(jù)可以發(fā)現(xiàn)：

與在標(biāo)準(zhǔn)臺式電腦上得到的結(jié)果相比較，我們看到基于Epiphany 系統(tǒng)的巨大優(yōu)勢——在一個負(fù)載平衡的系統(tǒng)中和采用更低的時鐘頻率，運(yùn)算速度卻可以提高5 倍以上。

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

進(jìn)一步研究， 在基于16 核的Epiphany 芯片在上運(yùn)行的并行程序的各階段的時間花費(fèi)，測量時間如下：

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

通過表格中的數(shù)據(jù)可以發(fā)現(xiàn)：

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

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

7 結(jié)束語

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

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

[2] 廖曉鐘，賴汝. 科學(xué)與工程計算[M]. 北京：電子工業(yè)出版社，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] 姜弘道，張健飛. 工程科學(xué)中的高性能計算[D]. 北京:科學(xué)出版社，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ī)群上的實現(xiàn)[J]. 計算機(jī)工程，2002（6）:102-130.