基于GPU的并行APSP問(wèn)題的研究

張凌潔，趙 英

（北京化工大學(xué) 北京 100029）

最短路徑問(wèn)題是圖論中的一個(gè)典型問(wèn)題，已經(jīng)得到廣泛的研究和應(yīng)用。此問(wèn)題包含兩類：1）單源最短路徑（SSSP，Single Source Shortest Paths），求給定始末頂點(diǎn)之間的最短路徑，如旅行商問(wèn)題；2）所有頂點(diǎn)間最短路徑（APSP，All Pair Shortest Paths），求圖中所有頂點(diǎn)對(duì)之間的全部最短路徑。

目前針對(duì)此問(wèn)題的相關(guān)研究成果已被廣泛應(yīng)用于眾多領(lǐng)域，如城市物流規(guī)劃、圖形學(xué)圖像分割等。由此可見(jiàn)最短路徑問(wèn)題不僅具有較高的科研意義，還具有相當(dāng)高的實(shí)際應(yīng)用價(jià)值。

另一方面，當(dāng)今社會(huì)已步入信息化時(shí)代。海量數(shù)據(jù)以及復(fù)雜的科學(xué)計(jì)算都對(duì)計(jì)算機(jī)處理能力提出了更高的要求。HPC已經(jīng)成為熱點(diǎn)研究問(wèn)題。其中GPU作為低成本低功耗的新型高性能技術(shù)得到廣泛關(guān)注和長(zhǎng)足的發(fā)展。尤其是在CUDA平臺(tái)出現(xiàn)后，打破了圖形處理編程與普通編程之間的壁壘，對(duì)于GPU的發(fā)展可謂是如虎添翼，使得更多的程序員跨入到GPU通用計(jì)算的研究當(dāng)中[1-2]。

在GPU成為研究熱點(diǎn)的背景下，利用其高性能和低成本的優(yōu)勢(shì)，通過(guò)對(duì)APSP最短路徑并行算法的研究，解決Floyd算法在GPU并行化中產(chǎn)生的問(wèn)題。并根據(jù)GPU特點(diǎn)，提出相應(yīng)改進(jìn)策略，從而進(jìn)一步加快計(jì)算速度。希望將來(lái)可以把研究成果運(yùn)用到如物流配送中心的選址，城市規(guī)劃等實(shí)踐中。

1 Floyd-W arshall算法簡(jiǎn)介

APSP問(wèn)題就是求解給定圖中所有頂點(diǎn)對(duì)之間的全部最短路徑。此問(wèn)題有兩種解法：

第1種解法是，每次以一個(gè)頂點(diǎn)作為起始頂點(diǎn)，依次將其他頂點(diǎn)作為終止頂點(diǎn)求它們之間的最短路徑，然后更換起始頂點(diǎn)，依次再求其他頂點(diǎn)間的最短路徑。如果對(duì)每對(duì)頂點(diǎn)之間都使用Dijkstra算法，則總的時(shí)間復(fù)雜度將是O（n4）。

第2種解法是Floyd-Warshall算法，也稱傳遞閉包算法。假定一張帶權(quán)的有向圖 G，包含 n個(gè)頂點(diǎn){Vi，1≤i≤n}，及 m條有向邊，有鄰接矩陣Arcs和初始距離矩陣 D（-1）。串行的Floyd Warshall算法的計(jì)算步驟如下：

1）初始化：初始化鄰接矩陣D，如果Vi和Vj之間存在直達(dá)路徑，則將D（-1）[i][j]的值設(shè)為路徑長(zhǎng)度；如果不存在直接路徑，則 D（-1）[i][j]的值設(shè)為 ∞。 注意，默認(rèn)情況下，頂點(diǎn)到頂點(diǎn)自身的路徑長(zhǎng)度為0；

2）依次計(jì)算 D（k）[i][j]（k=0，1，…n-1）中 D[i][j]的值，如果經(jīng)過(guò)中間頂點(diǎn) Vk，D[i][k]+D[k][j]的值小于D[i][j]，則表明經(jīng)過(guò)頂點(diǎn)k存在更短路徑，更新D[i][j]值，直至圖中的所有頂點(diǎn)都做過(guò)中間頂點(diǎn)得到最后的D矩陣。

使用Floyd-Warshall算法的串行版本去處理ASAP問(wèn)題，在數(shù)據(jù)規(guī)模比較小的時(shí)候，由于邏輯簡(jiǎn)單，還是非常實(shí)用的。然而，算法主體是三重循環(huán)，使得當(dāng)數(shù)據(jù)規(guī)模隨著圖中頂點(diǎn)數(shù)量增大的時(shí)候，算法的執(zhí)行時(shí)間也將呈現(xiàn)指數(shù)級(jí)的增加，在處理現(xiàn)實(shí)問(wèn)題時(shí)效果并不理想。

2 基于GPU的Floyd算法的并行化

Floyd-Warshall串行算法的核心是一個(gè)三重循環(huán)的迭代計(jì)算，算法的時(shí)間復(fù)雜度是O（n3）。如果將算法并行化，將極大的縮短計(jì)算時(shí)間[5]。對(duì)于三重循環(huán)體的處理，類似于GPU芯片對(duì)矩陣乘法的計(jì)算過(guò)程。所以，可以將循環(huán)體內(nèi)的比較權(quán)值的操作改為CUDA多線程并行執(zhí)行。實(shí)現(xiàn)多線程并行后可以大大提高計(jì)算速度，發(fā)揮GPU并行計(jì)算的優(yōu)勢(shì)。

該并行化算法的原理如下：首先選取距離矩陣中主對(duì)角線位置的元素作為Key，在n×n距離矩陣中共有n個(gè)Key，其坐標(biāo)值為（0，0），（1，1）…（n，n）。 每次迭代選取一個(gè) Key 作為基準(zhǔn)，選出與Key共行或共列的元素作為參考值，計(jì)算余下（n-1）×（n-1）個(gè)元素。每個(gè)元素 c計(jì)算時(shí)，從參考值中找到和本元素共行的元素a，和參考值中與本元素共列的元素b，判斷a+b與c的大小關(guān)系。如a+b更小，則更新該元素值為a+b；否則，不更新。計(jì)算完本輪的（n-1）×（n-1）個(gè)元素后，Key元素沿對(duì)角線向下移動(dòng)，按上述步驟重新計(jì)算。圖1展示了計(jì)算元素與參考元素的位置關(guān)系。

圖1 距離矩陣計(jì)算時(shí)元素對(duì)應(yīng)關(guān)系Fig．1 Relationship of elements in distancematrix

GPU并行程序的重點(diǎn)就是線程的分配和核函數(shù)（Kernel Function）的設(shè)計(jì)[3]。首先介紹如何通過(guò)核函數(shù)完成多線程的并行計(jì)算。并行化的Floyd-Warshall算法中使用到了以下幾個(gè)核函數(shù)：

1）初始化核函數(shù)：根據(jù)鄰接矩陣的值初始化距離矩陣中的權(quán)值。這里需要注意線程的對(duì)應(yīng)關(guān)系，將鄰接矩陣等價(jià)為一維數(shù)組，依次復(fù)制其中的值到grid中，并賦值給距離矩陣。所以每個(gè)線程在循環(huán)處理中的增量大小為整個(gè)grid中線程的個(gè)數(shù)。

2）計(jì)算核函數(shù)：計(jì)算距離矩陣中每個(gè)元素位置對(duì)應(yīng)的最短路徑。距離矩陣中的每個(gè)元素都與主軸上對(duì)應(yīng)的同行和同列的兩個(gè)元素的權(quán)值之和做比較。如，Dij表示頂點(diǎn)i到頂點(diǎn)j的路徑權(quán)值，需要和Dik+Dkj的權(quán)值作比較，其中Dik與Dij同行，Dkj與Dij同列，而k表示當(dāng)前迭代時(shí)是以頂點(diǎn)k作為中間頂點(diǎn)。

下面介紹線程的分配。這時(shí)會(huì)遇到下面的兩種情況：當(dāng)要處理的矩陣小于或等于grid的大小，那么在線程在計(jì)算相應(yīng)的頂點(diǎn)權(quán)值的時(shí)候，需要首先對(duì)是否在邊界當(dāng)中進(jìn)行判斷。如果在邊界當(dāng)中，則繼續(xù)進(jìn)行處理和計(jì)算，反之，則放棄計(jì)算；當(dāng)要處理的矩陣規(guī)模大于grid中線程個(gè)數(shù)時(shí)，需要計(jì)算其相對(duì)位置，分批次計(jì)算，如圖2所示。

圖2 距離矩陣計(jì)算時(shí)線程對(duì)應(yīng)關(guān)系Fig．2 Relationship of threads on grid， block and computing distancematrix

當(dāng)距離矩陣規(guī)模較大時(shí)，根據(jù)線程數(shù)分塊和相應(yīng)的映射原則，將數(shù)據(jù)交由GPU的grid計(jì)算。GPU屬于兩層并行結(jié)構(gòu)，其中每個(gè)grid中含有多個(gè)block，每個(gè)block又含有多個(gè)線程，上圖中左側(cè)的兩個(gè)圖就表示了這個(gè)關(guān)系[4]。當(dāng)grid中的線程數(shù)不足以為每個(gè)矩陣元素分配線程時(shí)，需要按照grid的大小將原距離矩陣分塊，依次分配給grid計(jì)算。在上面最右側(cè)的圖片中，1號(hào)子矩陣分配給grid計(jì)算完成后，再將2號(hào)子矩陣分配給grid計(jì)算，然后是3號(hào)子矩陣。注意，這里距離矩陣的大小并不能正好整除grid的大?。?%3！=0），所以3號(hào)子矩陣的元素個(gè)數(shù)比較少，可以分配給grid計(jì)算，剩余部分不做計(jì)算。以此類推，向左向下繼續(xù)計(jì)算各個(gè)子矩陣。這里使用相應(yīng)的行偏移和列偏移變量去表明GPU上grid應(yīng)該對(duì)應(yīng)的距離子矩陣的方位。

3 Floyd-W arshall并行算法的改進(jìn)

在上一節(jié)中，已經(jīng)實(shí)現(xiàn)了基于GPU的并行的Floyd-Warshall算法，但是為了進(jìn)一步提高程序的計(jì)算效率，可以利用GPU當(dāng)中的共享內(nèi)存對(duì)程序進(jìn)行改進(jìn)。這里借用了Venkataraman提出了一個(gè)更復(fù)雜的、基于 “塊”的Floyd的APSP算法[7]。對(duì)于大規(guī)模的圖而言，這個(gè)算法可以更好的利用緩存，并提供了更加具體的性能分析。在課題研究過(guò)程中，對(duì)這個(gè)的應(yīng)用進(jìn)行了擴(kuò)展，使其通過(guò)CUDA平臺(tái)運(yùn)行在GPU上，并獲得更高的效率。原始的F-W算法并不能夠拆分成相互之間無(wú)關(guān)聯(lián)的子矩陣，然后分給GPU上的每一個(gè)處理器來(lái)分開(kāi)計(jì)算。其原因就在于每一個(gè)子矩陣需要和整個(gè)數(shù)據(jù)集產(chǎn)生對(duì)應(yīng)關(guān)系。查看F-W算法串行，可以發(fā)現(xiàn)D矩陣中的每個(gè)值的更新都需要檢查矩陣中的其他的每一個(gè)數(shù)據(jù)。所以想要簡(jiǎn)單地拆分矩陣并不是可行的。但是，通過(guò)仔細(xì)選擇的一些子矩陣序列，就可以解決這個(gè)問(wèn)題。而這個(gè)分塊的F-W算法就是使用之前處理過(guò)的某些子矩陣去決定現(xiàn)在正需要處理的子矩陣，這就使得問(wèn)題并行化成為可能。

算法首先需要對(duì)距離矩陣進(jìn)行分塊，而劃分完成后，就可以開(kāi)始分階段進(jìn)行計(jì)算了。矩陣主模塊個(gè)數(shù)等于迭代的次數(shù)[7]。以下以6×6的矩陣為例，對(duì)這個(gè)算法進(jìn)行簡(jiǎn)要地介紹。首先如下圖將這個(gè)矩陣分出3個(gè)主模塊（主對(duì)角線上的3個(gè)子矩陣），所以就產(chǎn)生了3次迭代。每個(gè)階段根據(jù)自己的主模塊來(lái)運(yùn)算，算法的每一階段都分為同樣的3步計(jì)算。

圖3 分塊F-W算法過(guò)程演示Fig．3 Example of GPU-Tiled-FW parallel computing

步驟1：只計(jì)算主模塊中的最短路徑，觀察是否僅僅經(jīng)由主模塊所包含的頂點(diǎn)能夠獲得更短的路徑。此時(shí)將主模塊分配給CUDA中的一個(gè)block去計(jì)算，所以此步驟只需要用到了一個(gè)多處理器。主模塊就是當(dāng)前需要計(jì)算的模塊，所以除了主模塊本身，不在需要加載任何其他數(shù)據(jù)到共享存儲(chǔ)器（Shared Memory）。每一個(gè)線程可以直接加載它們的線程編號(hào)對(duì)應(yīng)的數(shù)據(jù)塊，并在這一階段的計(jì)算結(jié)束后，將結(jié)果保存到全局存儲(chǔ)器（Global Memory）中。如圖4所示，為當(dāng)?shù)?個(gè)主模塊為當(dāng)前主模塊時(shí)進(jìn)行步驟一時(shí)的情形。

步驟2：計(jì)算那些只依靠自己和相應(yīng)主模塊來(lái)運(yùn)算的子模塊，即與主模塊在同一行和同一列的模塊。此時(shí)，主模塊和當(dāng)前計(jì)算模塊都需要加載到共享存儲(chǔ)器中，這就使得GPU中每一個(gè)線程都能從block的共享存儲(chǔ)器中加載到一個(gè)對(duì)應(yīng)的單元。同樣，在此階段完成后，需要每個(gè)線程將數(shù)據(jù)傳回到全局存儲(chǔ)器中。如圖4（b）所示，為當(dāng)?shù)?個(gè)主模塊為當(dāng)前主模塊時(shí)進(jìn)行步驟2時(shí)的情形。

步驟3：計(jì)算剩余的子模塊，如圖4所示，當(dāng)以7號(hào)模塊為主模塊的計(jì)算到達(dá)步驟3時(shí)，只有1、2、3、4號(hào)模塊需繼續(xù)計(jì)算。也就是說(shuō)，當(dāng)完成步驟1和步驟2時(shí)，還要處理剩下的（n-1）×（n-1）個(gè)子模塊。此時(shí)需放入共享存儲(chǔ)器中的模塊除了自己之外，還需要將以當(dāng)前主模塊為軸的同行和同列對(duì)應(yīng)的模塊一并載入到共享存儲(chǔ)器當(dāng)中。如圖4（c）所示，為當(dāng)?shù)?個(gè)主模塊為當(dāng)前主模塊時(shí)進(jìn)行步驟3時(shí)的內(nèi)存分配情形。

完成上面3個(gè)步驟，沿對(duì)角線向下移動(dòng)主模塊，重復(fù)上面的計(jì)算步驟。直至全部主模塊都計(jì)算后，才算完成依次迭代。再進(jìn)行第二次迭代，比較兩次迭代后的距離矩陣：如果經(jīng)過(guò)后一次迭代，距離矩陣并未發(fā)生變化，則迭代結(jié)束；否則，進(jìn)行新一輪的迭代，直至距離矩陣不再發(fā)生改變。

圖4 內(nèi)存分配情況Fig．4 Memory allocation

4 實(shí)驗(yàn)結(jié)果比較與分析

文中所有的程序運(yùn)行的硬件平臺(tái)為惠普公司的ProLiant系列SL390s G7服務(wù)器。其中每個(gè)節(jié)點(diǎn)服務(wù)器都由兩個(gè)2．4-GHZ的英特爾Xeon系列E5620處理器以及24 GB大小的內(nèi)存和兩塊NVIDIA的Tesla系列M2050的GPU構(gòu)成。本課題的軟件平臺(tái)為 64位 Linux RedHat 6．0操作系統(tǒng)，NVIDIA Tesla 驅(qū)動(dòng)器，CUDA SDK 3．2，CUDA Toolkit 4．0。 所以的程序運(yùn)行數(shù)據(jù)都是在上述平臺(tái)上運(yùn)行15次取平均值獲得的。

從表1的數(shù)據(jù)和性能比較圖中可以看出，無(wú)論是原始的GPU并行F-W算法還是之后改進(jìn)的分塊的并行F-W算法，當(dāng)數(shù)據(jù)規(guī)模較小時(shí)（≤100）時(shí)，算法執(zhí)行時(shí)間維持在一個(gè)較低的水平，大概在3 s左右。而且從圖中也無(wú)法看出改進(jìn)的GPU并行程序優(yōu)越于原始的并行程序，可以認(rèn)為它們是相當(dāng)?shù)?。而串行的CPU版本的F-W算法卻極其出色的完成了任務(wù)，基本在0．5 s內(nèi)就完成計(jì)算。這可能是由于在GPU并行算法當(dāng)中，需要進(jìn)行GPU設(shè)備與主機(jī)內(nèi)存之間的數(shù)據(jù)傳輸，而這個(gè)傳輸是非常耗費(fèi)時(shí)間的，且由于計(jì)算規(guī)模太小，使得這一瓶頸突顯出來(lái)。

對(duì)于改進(jìn)的GPU并行F-W算法來(lái)說(shuō)，也是由于較小的計(jì)算量，使得算法中的分塊以及加載使用共享存儲(chǔ)器并沒(méi)有顯現(xiàn)出優(yōu)勢(shì)，所以使用兩種并行算法在效率上沒(méi)有什么分別。另外，值得注意的是，可以從圖5（a）中看到，串行版本的F-W算法在數(shù)據(jù)增加的后期，已經(jīng)開(kāi)始發(fā)生明顯地時(shí)間增加，盡管增加的絕對(duì)值并不大。

表1 光頂點(diǎn)數(shù)≤100的Floyd算法串行程序和GPU并行程序耗時(shí)比較Tab.1 Cost time comparison with less than 100 vertices on CPU，GPU，and Tiled-GPU

圖5 算法計(jì)算時(shí)間Fig．5 Computing time

表2 光頂點(diǎn)數(shù)≤1000的Floyd算法串行程序和GPU并行程序耗時(shí)比較Tab.2 Cost time comparison with less than 1000 vertices on CPU，GPU，and Tiled-GPU

表2的數(shù)據(jù)表示了在頂點(diǎn)數(shù)為100～1 000，間隔為100的輸入數(shù)據(jù)下，各種F-W算法計(jì)算時(shí)間的比較。首先從表2中，可以看出兩種GPU并行F-W算法在計(jì)算時(shí)間上只有相對(duì)較小幅度的增加，即從3 s左右的計(jì)算時(shí)間增加到了5 s左右。而且在頂點(diǎn)數(shù)大于等于700時(shí)，改進(jìn)后的分塊GPU并行F-W算法才穩(wěn)定地表現(xiàn)出了一些優(yōu)勢(shì)，盡管這個(gè)優(yōu)勢(shì)并不大。而CPU版本的串行程序的運(yùn)行時(shí)間在這個(gè)階段卻產(chǎn)生了指數(shù)級(jí)的增長(zhǎng)，這使得在圖5（b）中，GPU并行算法和GPU分塊并行算法的計(jì)算時(shí)間幾乎和坐標(biāo)軸平行。

表3的數(shù)據(jù)可以看出，在這個(gè)階段CPU對(duì)F-W算法的計(jì)算時(shí)間已經(jīng)超出需求范圍，這使得在圖5（c）中，并沒(méi)有出現(xiàn)串行版本的執(zhí)行時(shí)間曲線。很顯然，由于算法中的一個(gè)二維循環(huán)和一個(gè)三維循環(huán)，使得在這個(gè)階段，每一次的數(shù)據(jù)增加都會(huì)帶來(lái)迭代次數(shù)的“暴增”。而原始的GPU并行F-W算法也因?yàn)轫旤c(diǎn)數(shù)的增加，使其自身的迭代次數(shù)和循環(huán)處理次數(shù)增加。而且這一算法是直接讀取全局內(nèi)存中的數(shù)據(jù)，這種延遲也會(huì)累積到一個(gè)可觀的數(shù)值的。相對(duì)而言，改進(jìn)的分塊GPU并行F-W算法的計(jì)算時(shí)間增加就小得多。盡管這種算法當(dāng)然也會(huì)因?yàn)轫旤c(diǎn)數(shù)的增加，使其迭代次數(shù)增加，但是由于使用了共享存儲(chǔ)器的優(yōu)化，且訪問(wèn)共享存儲(chǔ)器是訪問(wèn)全局內(nèi)存的若干倍，這使得迭代次數(shù)的增加被訪問(wèn)的加速而掩蓋，進(jìn)而達(dá)到了令人滿意的并行效果。

表3 光頂點(diǎn)數(shù)≤8000的Floyd算法串行程序和GPU并行程序耗時(shí)比較Tab.3 Cost time comparison with less than 8000 vertices on CPU，GPU，and Tiled-GPU

5 結(jié) 論

文中針對(duì)APSP問(wèn)題中的Floyd-Warshall算法的并行化進(jìn)行了較為全面的研究和分析。實(shí)現(xiàn)了基于GPU的并行的Floyd算法，在數(shù)據(jù)規(guī)模超過(guò)100個(gè)頂點(diǎn)后，加速比平均達(dá)到100倍以上。然后又利用GPU架構(gòu)中的共享儲(chǔ)存器，根據(jù)Venkataraman提出的分塊思想，將上述GPU并行算法進(jìn)行優(yōu)化。進(jìn)一步優(yōu)化之后，加速比達(dá)到8倍以上，計(jì)算速度明顯加快。

在后續(xù)的研究中，希望可以將GPU版本的并行F-W算法應(yīng)用在多GPU上，得到更高的性能加速比。并將這一并行算法應(yīng)用到實(shí)際中，如物流中心的選址、城市規(guī)劃等問(wèn)題當(dāng)中。

[1]許楨．關(guān)于CPU+GPU異構(gòu)計(jì)算的研究與分析[J]．科技信息，2010（17）：22-28．

XU Zhen．Research and analysis of CPU+GPU heterogeneous computing[J]．Technology Information，2010（17）：22-28．

[2]張舒，褚艷麗，趙開(kāi)勇，等．GPU高性能運(yùn)算之CUDA[M]．北京：中國(guó)水利出版社，2009．

[3]Harish P，Narayanan P J．Accelerating large graph algorithms on the GPU using CUDA[J]．High Performance Computing-HiPC，2007：197-208．

[4]程豪，張?jiān)迫?，張先軼．CPU-GPU并行矩陣乘法的實(shí)現(xiàn)與性能分析[J]．計(jì)算機(jī)工程，2010（7）：24-25．

CHENG Hao，ZHANG Yun-quan，ZHANG Xian-yi．Implementation and performance analysis of CPU-GPU parallelmatrix multiplication[J]．Computer Engineering，2010（7）：24-25．

[5]Nvidia：Cuda 3．0．Technical report，Nvidia Corp．（2010） Available at[EB/OL]．http：//developer．nvidia．com/object/cuda_3_0_downloads．htm l．

[6]盧照，師軍．并行最短路徑搜索算法的設(shè)計(jì)與實(shí)現(xiàn)[J]．計(jì)算機(jī)工程與應(yīng)用，2010（46）：69-71．

LU Zhao，SHIJun．Design and implement of parallel computing of shortest paths[J]．Computer Engineering and Applications，2010（46）：69-71．

[7]Venkataraman G，Sahni S，Mukhopad H S．A blocked all-pairs shortest-pathsalgorithm[M]．J．Exp．Algorithmics，2003．