張凌潔，趙 英

（北京化工大學 北京 100029）

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

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

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

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

1 Floyd-W arshall算法簡介

APSP問題就是求解給定圖中所有頂點對之間的全部最短路徑。此問題有兩種解法：

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

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

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

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

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

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

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

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

圖1 距離矩陣計算時元素對應關系Fig．1 Relationship of elements in distancematrix

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

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

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

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

圖2 距離矩陣計算時線程對應關系Fig．2 Relationship of threads on grid， block and computing distancematrix

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

3 Floyd-W arshall并行算法的改進

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

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

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

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

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

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

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

圖4 內存分配情況Fig．4 Memory allocation

4 實驗結果比較與分析

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

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

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

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

圖5 算法計算時間Fig．5 Computing time

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

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

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

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

5 結 論

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

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

[1]許楨．關于CPU+GPU異構計算的研究與分析[J]．科技信息，2010（17）：22-28．

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

[2]張舒，褚艷麗，趙開勇，等．GPU高性能運算之CUDA[M]．北京：中國水利出版社，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]程豪，張云泉，張先軼．CPU-GPU并行矩陣乘法的實現(xiàn)與性能分析[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]盧照，師軍．并行最短路徑搜索算法的設計與實現(xiàn)[J]．計算機工程與應用，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．