肖 漢，肖詩洋，李煥勤，周清雷

(1.鄭州師范學(xué)院 信息科學(xué)與技術(shù)學(xué)院，河南 鄭州 450044；2.東南大學(xué) 土木工程學(xué)院，江蘇 南京 211189；3.鄭州大學(xué) 計(jì)算機(jī)與人工智能學(xué)院，河南 鄭州 450001)

圖分析算法在許多領(lǐng)域具有重要應(yīng)用價(jià)值，全源對(duì)最短路徑算法是一種重要的圖最短路徑算法，在生物網(wǎng)絡(luò)比對(duì)和路徑推斷、交通運(yùn)輸、運(yùn)籌學(xué)、機(jī)器人動(dòng)作規(guī)劃、網(wǎng)絡(luò)通信等領(lǐng)域有廣泛的應(yīng)用.但是，隨著最短路徑問題應(yīng)用領(lǐng)域的數(shù)據(jù)規(guī)模迅速增長，數(shù)據(jù)結(jié)構(gòu)日益復(fù)雜，數(shù)據(jù)處理時(shí)間過長，限制了其在實(shí)際中的應(yīng)用[1].同時(shí)，全源對(duì)最短路徑處理算法的計(jì)算量大，如何提高最短路徑處理速度是當(dāng)前圖算法處理領(lǐng)域的研究熱點(diǎn)之一[2].

當(dāng)前，圖形處理器（Graphics Processing Unit,GPU）的浮點(diǎn)運(yùn)算處理能力是CPU 的10 倍以上，帶寬是CPU 的5 倍，但成本和能耗遠(yuǎn)低于CPU，在許多方面具有很大的優(yōu)勢[3].架構(gòu)方面，GPU 已經(jīng)采用了統(tǒng)一計(jì)算架構(gòu)單元，并實(shí)現(xiàn)了線程間通信和細(xì)粒度的任務(wù)分配，在通用計(jì)算領(lǐng)域GPU 已經(jīng)顯示出了強(qiáng)大的競爭力.統(tǒng)一計(jì)算設(shè)備架構(gòu)（Compute Unified Device Architecture,CUDA）是一種非開放標(biāo)準(zhǔn)，僅適用于在NVIDIA 的GPU 硬件架構(gòu)上運(yùn)行，從而阻礙了CUDA 技術(shù)的推廣.然而，隨著基于GPU 的通用計(jì)算（General-Purpose computing on GPU,GPGPU）被人們廣泛接受，業(yè)界已經(jīng)針對(duì)GPGPU 制定了如開放式計(jì)算語言（Open Computing Language，OpenCL）的開放標(biāo)準(zhǔn).設(shè)計(jì)適用在各種GPU 上且能夠高效執(zhí)行的基于OpenCL 的最短路徑并行算法是一項(xiàng)有挑戰(zhàn)性且非常有意義的工作[4-5].

本文利用OpenCL 架構(gòu)實(shí)現(xiàn)了一種高效的基于CPU+GPU 的全源對(duì)最短路徑并行算法（OCL_SP），解決了全源對(duì)最短路徑算法在不同異構(gòu)計(jì)算設(shè)備上快速處理大規(guī)模頂點(diǎn)集的問題.OCL_SP 并行算法在各種設(shè)備上表現(xiàn)出良好的可移植性.

1 相關(guān)研究

最短路徑算法是許多復(fù)雜圖算法的一個(gè)重要組成部分，普遍應(yīng)用于科學(xué)研究和工程分析中.Jamour 等[6]將圖表分解為雙連通組件，圖形中的所有點(diǎn)對(duì)最短路徑的計(jì)算性能提高3.7 倍.Zhang等[7]采用方向共享、并行請求和路由點(diǎn)優(yōu)化減少外部路由請求的數(shù)量，提出了一種能在外部路由請求中找到最佳路點(diǎn)集的貪心算法.Chandio 等[8]通過遵循批量同步并行參數(shù)預(yù)先計(jì)算路段的最短路徑距離和速度約束，提出了一種基于云環(huán)境的實(shí)時(shí)GPS 軌跡的完全自適應(yīng)地圖匹配策略.Panomruttanarug[9]利用經(jīng)典的路徑規(guī)劃方法找到最短的停車路徑，并利用經(jīng)驗(yàn)學(xué)習(xí)的能力跟蹤設(shè)計(jì)路徑.Hung 等[10]開發(fā)了一種分布式、并行且可擴(kuò)展的最短路徑算法消除網(wǎng)絡(luò)中的冗余間接邊緣.Gyongyosi[11]在量子中繼器網(wǎng)絡(luò)中以群體智能為基礎(chǔ)，使用并行線程執(zhí)行路由以通過糾纏梯度系數(shù)確定最短路徑.閆春望等[12]提出了基于自動(dòng)波理論的并行模糊神經(jīng)網(wǎng)絡(luò)最短路徑算法，其迭代次數(shù)和收斂速度得到優(yōu)化.李平等[13]采用雙線程進(jìn)行子網(wǎng)分割和拼接匯總得到最短路徑，獲得了2.07倍加速比.

Stpiczynski 等[14]通過基于Cilk 數(shù)組符號(hào)和內(nèi)聯(lián)函數(shù)的手動(dòng)矢量化技術(shù)，實(shí)現(xiàn)了Belman-Ford 單源最短路徑并行算法.Cabrera 等[15]提出了以稀疏矩陣-向量乘法的迭代為基礎(chǔ)的單源最短路徑算法，并行DBMS 系統(tǒng)獲得了3.8 倍加速比.Maleki 等[16]采用具備16 個(gè)核心節(jié)點(diǎn)的共享和分布式存儲(chǔ)器系統(tǒng)，設(shè)計(jì)了一種帶狀放松式Dijkstra 最短路徑并行算法，取得了1.66 倍加速比.Chakaravarthy 等[17]采用邊緣分類和方向優(yōu)化，并提出負(fù)載平衡策略處理更高度的頂點(diǎn)，實(shí)現(xiàn)了基于CPU 集群平臺(tái)的單源最短路徑并行算法.孫文彬等[18]采用多粒度通訊的策略，實(shí)現(xiàn)了基于MPI 的Dijkstra 最短路徑并行算法.

Liu 等[19]通過解決在稀疏矩陣-矩陣乘法中3個(gè)不規(guī)則性問題，最短路徑問題在GPU 平臺(tái)上獲得了大約兩倍加速比.李寅等[20]提出一種基于GPU架構(gòu)的采樣混合式全源對(duì)最短路徑并行算法，算法采用BFS 遍歷方法，能達(dá)到7.2 倍的加速比.葉穎詩等[21]設(shè)計(jì)實(shí)現(xiàn)了一種并行多標(biāo)號(hào)Dijkstra 算法，獲得了3.49 倍的加速比.Djidjev 等[22]提出了一種基于GPU 集群的Floyd-Warshall 最短路徑并行算法，取得了一個(gè)數(shù)量級(jí)的性能提升.

Aridhi 等[23]采用并行方式求解原始圖中每個(gè)子圖上的最短路徑，提出了一種基于MapReduce的解決道路網(wǎng)絡(luò)中的最短路徑問題.Gayathri 等[24]采用結(jié)合傳遞閉包屬性和Dijkstra 單源最短路徑算法中貪婪技術(shù)特征尋找到所有點(diǎn)對(duì)最短路徑，在MapReduce 上設(shè)計(jì)了ex-FTCD 算法.何亞茹等[25]在基于神威平臺(tái)的單個(gè)SW26010 處理器上，利用MPI 技術(shù)求解全源最短路徑的Floyd 并行算法，取得了106 倍的最高加速比.

綜上，最短路徑并行算法研究有些以新型理論為基礎(chǔ)提出新的最短路徑算法，有些側(cè)重對(duì)Dijkstra 等傳統(tǒng)算法本身進(jìn)行優(yōu)化，有些基于CPU集群計(jì)算、CUDA 和GPU 集群等并行計(jì)算方式實(shí)現(xiàn)最短路徑并行算法，有些則以大數(shù)據(jù)Hadoop 平臺(tái)為基礎(chǔ)設(shè)計(jì)最短路徑高性能算法.可是，全源對(duì)最短路徑問題采用多種算法多類平臺(tái)系統(tǒng)進(jìn)行性能評(píng)估報(bào)道目前較少.多數(shù)研究的加速效果不顯著且均受限于一種算法或平臺(tái).本文將根據(jù)OpenCL框架特征和最短路徑問題的特點(diǎn)，研究在GPU 加速下的最短路徑算法和在不同并行計(jì)算平臺(tái)下的算法移植性能.

2 最短路徑算法分析

2.1 GPU 加速處理通用計(jì)算OpenCL 是一種開放的、免版稅的標(biāo)準(zhǔn)，用于對(duì)個(gè)人計(jì)算機(jī)、服務(wù)器、移動(dòng)設(shè)備和嵌入式平臺(tái)中的各種處理器進(jìn)行跨平臺(tái)、并行編程.OpenCL 定義了具有模型層次結(jié)構(gòu)的計(jì)算體系結(jié)構(gòu)：平臺(tái)模型、內(nèi)存模型、執(zhí)行模型和編程模型.OpenCL 平臺(tái)由主機(jī)和連接到主機(jī)的計(jì)算設(shè)備組成[26].主機(jī)通常是CPU，OpenCL 設(shè)備可以是支持OpenCL 的設(shè)備，如CPU、GPU 和FPGA.OpenCL 設(shè)備執(zhí)行內(nèi)核，每個(gè)內(nèi)核都是程序代碼中的一個(gè)函數(shù).內(nèi)核以工作項(xiàng)為單位執(zhí)行，工作項(xiàng)是并行執(zhí)行的單位.整個(gè)工作項(xiàng)集合稱為一個(gè)N-Dimensional Range (NDRange).共享其結(jié)果的工作項(xiàng)被組織成工作組的一個(gè)單元[27].

支持OpenCL 的計(jì)算設(shè)備由OpenCL 計(jì)算架構(gòu)中的主機(jī)端系統(tǒng)統(tǒng)一管理和調(diào)度.當(dāng)kernel 被主機(jī)端提交到設(shè)備端時(shí)，OpenCL 為任務(wù)并行處理定義了索引空間NDRange.執(zhí)行內(nèi)核的各個(gè)工作項(xiàng)的標(biāo)識(shí)由其在NDRange 中的坐標(biāo)表示.多個(gè)工作項(xiàng)可組織為工作組，提供了對(duì)NDRange 更粗粒度的劃分，同一工作組中的多個(gè)工作項(xiàng)將在一個(gè)計(jì)算單元的不同處理單元上并發(fā)執(zhí)行[28].

2.2 算法原理矩陣乘法求最短路徑算法的核心思想是通過多次迭代計(jì)算求得點(diǎn)對(duì)間的最短距離.在該算法中是將原矩陣乘法中的“×”運(yùn)算以“+”運(yùn)算代替，“∑”運(yùn)算以“min”運(yùn)算代替[29].文中矩陣乘法運(yùn)算內(nèi)涵均按照此處說明進(jìn)行.而重復(fù)平方法是對(duì)矩陣乘法的變形應(yīng)用，重復(fù)平方法求最短路徑是對(duì)每次計(jì)算得到的結(jié)果矩陣本身進(jìn)行“×”運(yùn)算[30].

在重復(fù)平方法最短路徑算法處理的各個(gè)步驟中，建立有向圖的帶權(quán)鄰接矩陣步驟、對(duì)鄰接矩陣進(jìn)行初始化步驟和建立鄰接矩陣的鏡像矩陣步驟的時(shí)間復(fù)雜度均為O(n2).利用鄰接矩陣計(jì)算有向圖各頂點(diǎn)間的最短路徑步驟的時(shí)間復(fù)雜度為.因此，本文主要關(guān)注頂點(diǎn)間的最短路徑階段的并行性及優(yōu)化.

2.3 可并行性分析存在于算法中的數(shù)據(jù)關(guān)聯(lián)性與算法的可并行性好壞有關(guān).若數(shù)據(jù)運(yùn)算先后關(guān)聯(lián)性越強(qiáng)，算法可并行性越低.反之，若數(shù)據(jù)運(yùn)算先后關(guān)聯(lián)性越弱，算法可并行性越高，并行化后算法的性能會(huì)越好.對(duì)加權(quán)有向圖的重復(fù)平方法最短路徑算法的可并行性分析如下，如圖1 所示.

圖1 加權(quán)有向圖Fig.1 Weighted digraph

圖1 中有權(quán)鄰接矩陣L(1)如式(1)所示.按照重復(fù)平方法計(jì)算，則有權(quán)鄰接矩陣L(2)如式(2)所示（“ ⊕”符號(hào)代表兩個(gè)矩陣進(jìn)行重復(fù)平方法運(yùn)算）.在計(jì)算L(2)的過程中，發(fā)現(xiàn)L(2)中一個(gè)任意元素的計(jì)算和其余元素的計(jì)算過程相互無關(guān)，不存在依賴性.即在L(2)中計(jì)算某個(gè)元素值時(shí)，可對(duì)其余元素值同時(shí)進(jìn)行運(yùn)算.結(jié)合OpenCL 模型的特性，在工作項(xiàng)中進(jìn)行每個(gè)元素的計(jì)算，這樣，L(2)中每個(gè)元素的結(jié)果將由工作項(xiàng)計(jì)算得出.如果計(jì)算完矩陣中的每個(gè)元素，則本輪結(jié)束計(jì)算，如果需要迭代，則重復(fù)上面的過程.

3 最短路徑并行算法的設(shè)計(jì)和優(yōu)化

3.1 并行算法描述設(shè)有相同鄰接矩陣A、B，鄰接矩陣A、B中的子矩陣分別存儲(chǔ)于數(shù)組As和數(shù)組Bs中，每次計(jì)算之后的最小值由v保存，矩陣C存放的是完成計(jì)算后點(diǎn)對(duì)間最短路徑長度數(shù)據(jù).重復(fù)平方法最短路徑并行算法描述見算法1.

算法1 SIMT 模型上的重復(fù)平方法最短路徑并行算法

矩陣乘法是重復(fù)平方法最短路徑并行算法的核心運(yùn)算.每一個(gè)工作組負(fù)責(zé)矩陣C中v的計(jì)算，每一個(gè)工作項(xiàng)計(jì)算比較得出v的最小元素值.每個(gè)工作組中的工作項(xiàng)被劃分成為Bsize×Bsize的正方形計(jì)算工作項(xiàng)，與正方形計(jì)算工作項(xiàng)對(duì)應(yīng)的矩陣A、B中的數(shù)據(jù)就是需要進(jìn)行計(jì)算的源數(shù)據(jù).在每次進(jìn)行運(yùn)算前，將源數(shù)據(jù)存入本地存儲(chǔ)器As和Bs中.在每個(gè)工作組內(nèi)進(jìn)行計(jì)算時(shí)，工作組ID 不變，工作項(xiàng)ID 也不變，每循環(huán)一次需要更新一次本地?cái)?shù)組As和Bs中存放的數(shù)據(jù).當(dāng)As和Bs中數(shù)據(jù)更新完畢，對(duì)于As中的每一行元素、Bs中的每一列元素，由每個(gè)工作項(xiàng)負(fù)責(zé)數(shù)組As和數(shù)組Bs中對(duì)應(yīng)元素的求和、取最小值運(yùn)算，并將取得的最小值存入v中.當(dāng)該正方形中的元素全部讀取計(jì)算完畢，則讀取該計(jì)算條帶中的下一個(gè)正方形區(qū)域源數(shù)據(jù)，直到每個(gè)工作組中的元素全部讀取計(jì)算結(jié)束.工作組中的正方形計(jì)算區(qū)域每次移動(dòng)的步長為Bsize.當(dāng)工作組中所有元素讀取計(jì)算完畢時(shí)，v中保存的就是這些矩陣中對(duì)應(yīng)元素和的最小值.直到所有工作組中的所有線程全部計(jì)算結(jié)束，把最終得到的結(jié)果再從本地存儲(chǔ)器寫回全局存儲(chǔ)器.

3.2 算法并行化思路OCL_SP 算法具體執(zhí)行步驟如下.

步驟 1根據(jù)頂點(diǎn)數(shù)，在主機(jī)端初始化鄰接矩陣A并保存.

步驟 2初始化OpenCL 平臺(tái).

步驟 3命令序列對(duì)象通過上下文及指定設(shè)備創(chuàng)建.在設(shè)備和上下文之間通過clCreateCommandQueue建立邏輯連接，以協(xié)調(diào)內(nèi)核計(jì)算.

步驟 4輸入源文件并創(chuàng)建、編譯程序?qū)ο?依據(jù)在上下文中的設(shè)備特性，通過運(yùn)行時(shí)編譯構(gòu)建程序?qū)ο?，OCL_SP 并行算法具備了較強(qiáng)的可移植性.

步驟 5設(shè)置OpenCL 內(nèi)存對(duì)象和傳輸數(shù)據(jù).首先創(chuàng)建buffer 存儲(chǔ)器對(duì)象，接著將存儲(chǔ)器的訪問任務(wù)加載到命令隊(duì)列中，最后利用clCreateBuffer將鄰接矩陣A、B隱式的從主機(jī)寫入到設(shè)備的全局內(nèi)存中.

步驟 6構(gòu)建kernel 對(duì)象.利用clCreateKernel將內(nèi)核函數(shù)與參數(shù)封裝到指定的kernel 對(duì)象中.

步驟 7為kernel 對(duì)象設(shè)置參數(shù)傳遞.

步驟 8構(gòu)建頂點(diǎn)數(shù)目滿足的條件，創(chuàng)建內(nèi)核函數(shù)，調(diào)度內(nèi)核執(zhí)行.

步驟 9循環(huán)調(diào)用內(nèi)核函數(shù)對(duì)數(shù)據(jù)的處理，計(jì)算點(diǎn)對(duì)間的最短路徑距離值.

步驟 10在設(shè)備端中完成運(yùn)算任務(wù)后，將結(jié)果傳輸?shù)街鳈C(jī)端存儲(chǔ)器，釋放設(shè)備端內(nèi)存空間；并將最終計(jì)算結(jié)果保存至相應(yīng)文件.

3.3 并行算法的加速策略

3.3.1 理論設(shè)計(jì) 設(shè)計(jì)重復(fù)平方法最短路徑并行算法時(shí)，矩陣乘法采用了混合劃分法實(shí)現(xiàn).混合劃分法由帶狀劃分法和棋盤劃分法組成，在工作空間NDRange 中采取的是帶狀劃分法，在每個(gè)工作組內(nèi)采取的是棋盤劃分法.

工作空間中的帶狀劃分法如圖2 所示，其中一個(gè)工作組對(duì)應(yīng)矩陣A和 矩陣B中一對(duì)計(jì)算條帶的運(yùn)算，運(yùn)算完成后的結(jié)果放在矩陣C的對(duì)應(yīng)位置上.在矩陣A中陰影區(qū)域計(jì)算條帶的高度為Bsize，矩陣B中陰影區(qū)域計(jì)算條帶的寬度為Bsize.在每個(gè)計(jì)算條帶中，再對(duì)數(shù)據(jù)元素進(jìn)行棋盤式的劃分.每次計(jì)算時(shí)，塊坐標(biāo)不變，但是每次讀取數(shù)據(jù)的步長按照Bsize增加，以更新將要計(jì)算的數(shù)據(jù).計(jì)算條帶內(nèi)劃分方法如圖3 所示.

圖2 NDRange 帶狀劃分圖Fig.2 NDRange banded partition diagram

圖3 工作組內(nèi)棋盤劃分圖Fig.3 Chessboard division diagram in the working group

在進(jìn)行矩陣乘法運(yùn)算時(shí)，工作組每次加載一個(gè)計(jì)算條帶.由于計(jì)算條帶寬度為Bsize，所以每次進(jìn)行迭代的步長是Bsize.根據(jù)帶狀劃分的方法，矩陣A從全局存儲(chǔ)器完全被加載到本地存儲(chǔ)器，需要被工作組內(nèi)的本地存儲(chǔ)器讀取Bwidth/Bsize次 .而和矩陣A進(jìn)行矩陣乘法運(yùn)算的矩陣B需要被工作組內(nèi)的本地存儲(chǔ)器讀取的次數(shù)是Aheight/Bsize次，這樣的數(shù)據(jù)讀取方式大量減少了直接從全局存儲(chǔ)器讀取數(shù)據(jù)的延遲時(shí)間.工作組內(nèi)按照棋盤劃分法設(shè)計(jì)矩陣乘法.每個(gè)工作項(xiàng)負(fù)責(zé)讀取矩陣A中的一行元素和矩陣B中的一列元素，計(jì)算得出矩陣C中對(duì)應(yīng)的元素值，并儲(chǔ)存在全局存儲(chǔ)器中.

對(duì)于每一個(gè)工作組，都需要由A的一行元素和B的一列元素進(jìn)行比較計(jì)算之后才能得到最終結(jié)果，如圖4 所示.矩陣A中高度為Bsize長度為Awidth的矩形代表一個(gè)工作組，U和U′分別代表讀取到本地存儲(chǔ)器As和Bs的正方形區(qū)域源數(shù)據(jù).每個(gè)工作組內(nèi)的操作過程如下：

圖4 矩陣乘法的分塊計(jì)算Fig.4 Block calculation of matrix multiplication

步驟 1把U和U′從全局存儲(chǔ)器中讀出，放入本地存儲(chǔ)器As和Bs中 .在邊長為Bsize的正方形計(jì)算區(qū)域中，一個(gè)工作項(xiàng)對(duì)應(yīng)加載一對(duì)元素，然后利用矩陣乘法計(jì)算比較得出最小的元素值，并存入v.

步驟 2進(jìn)行for 循環(huán)，以邊長Bsize為步長進(jìn)行迭代，更新As和Bs中保存的數(shù)據(jù).重復(fù)步驟1 直到讀取完該計(jì)算條帶A的每行元素和相應(yīng)計(jì)算條帶B的每列元素.

步驟 3在最后一次內(nèi)層for 循環(huán)中，W和W′被放入As和Bs中，在最后一輪比較運(yùn)算中所有工作項(xiàng)將所有計(jì)算任務(wù)完成后，每個(gè)工作項(xiàng)計(jì)算比較得出的最終結(jié)果保存在v.相應(yīng)在圖4(c)中，計(jì)算得出了圖中的塊標(biāo)號(hào)為(0,0)的塊中的數(shù)據(jù).

步驟 4對(duì)于矩陣C中的每個(gè)小塊，都進(jìn)行上述相同的操作，直到全部計(jì)算結(jié)束，矩陣C中保存的就是該輪計(jì)算得出的結(jié)果.

3.3.2 kernel 函數(shù)的配置設(shè)計(jì) 整個(gè)鄰接矩陣A依據(jù)互不重疊的劃分原則，可以被分割為若干個(gè)計(jì)算條帶.計(jì)算條帶作為基本處理單位，存在于計(jì)算條帶中的大量矩陣乘法計(jì)算可交給OpenCL 工作組和工作項(xiàng)處理.文中采用二維索引空間來設(shè)計(jì)，從數(shù)據(jù)角度來講，二維索引空間在x、y方向上的維度均為D，工作組在x、y方向上的維度均為Bsize，工作空間在方向x、y方向上工作項(xiàng)總量均為base×Bsize.因此，內(nèi)核函數(shù)的網(wǎng)格配置如下：

3.4 優(yōu)化設(shè)計(jì)

3.4.1 本地內(nèi)存的設(shè)計(jì) GPU 計(jì)算性能常受制于存儲(chǔ)器速度.由于運(yùn)算速度與顯存帶寬間的不平衡，OpenCL 定義了多層存儲(chǔ)器體系結(jié)構(gòu).文中依據(jù)算法特征將存儲(chǔ)器的合理分配和優(yōu)化運(yùn)用于計(jì)算流程中的訪存方式和數(shù)據(jù)傳輸.在重復(fù)平方法最短路徑算法的矩陣乘法運(yùn)算中，鄰接矩陣An×n、Bn×n中的每一個(gè)數(shù)據(jù)都需要被重復(fù)讀取n次.如果將計(jì)算條帶放在具有400～800 個(gè)時(shí)鐘周期讀寫延遲的全局存儲(chǔ)器中，將大大降低并行數(shù)據(jù)訪問速率.通過引入僅具有1（無bank 沖突）～16 個(gè)時(shí)鐘周期（發(fā)生16 路bank 沖突）讀寫延遲的本地存儲(chǔ)器，將需要大量重用的計(jì)算條帶數(shù)據(jù)在計(jì)算前預(yù)置在本地內(nèi)存As和Bs中.這樣既能方便地在工作項(xiàng)間進(jìn)行通信，也使工作項(xiàng)的訪存速度更快從而大量的全局存儲(chǔ)器延遲被透明地隱藏起來.因此，并行算法中采用靜態(tài)分配方式在內(nèi)核函數(shù)中對(duì)As和Bs分配本地存儲(chǔ)器空間：

同時(shí)，在不同的本地存儲(chǔ)器 bank 中存儲(chǔ)了每條計(jì)算條帶中的數(shù)據(jù)，保證了在同一個(gè)warp 中相鄰的一組32 個(gè)訪問工作項(xiàng)與bank 是一一對(duì)應(yīng)，避免了bank 沖突，有效提高了運(yùn)算速度.

3.4.2 NDRange 優(yōu)化 計(jì)算單元應(yīng)該保持較高的占用率.因?yàn)楣ぷ黜?xiàng)指令在OpenCL 中是順序執(zhí)行，隱藏延遲的有效方法就是在一個(gè)warp 繁忙時(shí)執(zhí)行其他warp.若存在足夠多的工作組數(shù)量，就可避免因warp 塊填充不足導(dǎo)致OpenCL 系統(tǒng)的計(jì)算資源浪費(fèi).一個(gè)工作組中保持的工作項(xiàng)最好是32的倍數(shù).由于有向圖的頂點(diǎn)數(shù)是變化的，工作組和工作組中的工作項(xiàng)數(shù)應(yīng)該根據(jù)所占用的寄存器空間及本地存儲(chǔ)器空間動(dòng)態(tài)分配，以達(dá)到最佳效率.針對(duì)OpenCL 加速的全源對(duì)最短路徑圖算法，當(dāng)有向圖頂點(diǎn)數(shù)為5 000 時(shí)，不同工作組維度下，算法運(yùn)算時(shí)間比較如表1 所示.從表1 中可見，當(dāng)工作組維度為16×16 時(shí)，全源對(duì)最短路徑并行算法的運(yùn)算時(shí)間最短.

表1 算法運(yùn)算時(shí)間比較Tab.1 Compare the operation times of algorithm

3.5 其他優(yōu)化方案為了便于各種并行化方案的性能比較，本文分別在OpenMP 和CUDA 并行計(jì)算體系結(jié)構(gòu)上實(shí)現(xiàn)了一種符合2.2 節(jié)最短路徑算法原理的并行算法.

3.5.1 基于OpenMP 的最短路徑并行算法（OMP_SP）在Windows 系統(tǒng)中，使用OpenMP 內(nèi)置的Microsoft Visual Studio.Net 2017 編譯器，在項(xiàng)目的Microsoft.Cpp.x64.user 屬性頁中設(shè)置‘C/C++→Language→OpenMP Support’為“是(/OpenMP)”，在Visual Studio.net 2017 中打開OpenMP 編譯選項(xiàng).然后將選項(xiàng)“C/C++→Gode Generation→Runtime Library”更改為“多線程調(diào)試(/mtd)”.

OpenMP 是一種粗粒度并行.根據(jù)CPU 核心的數(shù)量，并行塊通過omp_set_num_thread (thread)啟用子線程進(jìn)行并行執(zhí)行.主線程執(zhí)行OMP_SP并行計(jì)算系統(tǒng)，當(dāng)遇到位于構(gòu)造初始鄰接矩陣L(1)、保存L(1)到A和B、矩陣乘法等函數(shù)之前的編譯引導(dǎo)指令#pragma omp for 時(shí)，將對(duì)后面的循環(huán)體進(jìn)行多線程分解.每個(gè)線程并行處理結(jié)果矩陣的一個(gè)元素.

3.5.2 基于CUDA 的最短路徑并行算法（CUDA_SP）CUDA 是一個(gè)細(xì)粒度的平行.將最短路徑算法可并行的矩陣乘法功能部分劃分為一個(gè)任務(wù).該任務(wù)由相應(yīng)的global 函數(shù)處理.在算法運(yùn)行過程中，global 函數(shù)與網(wǎng)格一一對(duì)應(yīng)，網(wǎng)格將任務(wù)分配給線程塊，線程塊將任務(wù)重新分配給線程進(jìn)行處理.主機(jī)端調(diào)用global 函數(shù)，在GPU 端的多個(gè)計(jì)算處理單元中進(jìn)行并行處理.

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

4.1 測試環(huán)境和數(shù)據(jù)記錄選擇單精度數(shù)據(jù)類型實(shí)現(xiàn)是因?yàn)樗臏?zhǔn)確性，而且GPU 針對(duì)現(xiàn)代計(jì)算機(jī)在單精度浮點(diǎn)運(yùn)算上進(jìn)行了高度優(yōu)化.

用OpenCL C 語言開發(fā)了全源對(duì)最短路徑并行算法的核心源代碼.該并行算法由兩部分構(gòu)成，在主機(jī)上執(zhí)行主程序部分，另外內(nèi)核程序部分是在OpenCL 設(shè)備上執(zhí)行.核心代碼實(shí)現(xiàn)了矩陣乘法運(yùn)算.

實(shí)驗(yàn)工作是在CPU/GPU 異構(gòu)混合并行計(jì)算平臺(tái)上進(jìn)行.硬件環(huán)境參數(shù)如下：平臺(tái)1：CPU 為AMD Ryzen5 1600X 3.6 GHz 六核；系統(tǒng)內(nèi)存容量為24 GB；顯卡采用NVIDIA Geforce GTX 1070 GPU.GTX 1070 擁有1 920 個(gè)CUDA 核，核心頻率是1 506 MHz.搭配8 GB GDDR5 顯存容量，顯存位寬256 位.平臺(tái)2：CPU 為AMD Ryzen5 1600X 3.6 GHz 六核；系統(tǒng)內(nèi)存容量為24 GB；顯卡采用AMD Radeon RX 570 GPU.Radeon RX 570 擁有2 048 個(gè)核心，核心頻率是1 168 MHz.搭配8 GB GDDR5 顯存容量，顯存位寬256 位.

軟件平臺(tái)：采用Microsoft Windows 10 64 位為操作系統(tǒng)；Microsoft Visual Studio.Net 2017 為集成開發(fā)環(huán)境；CUDA Toolkit 10.0 為系統(tǒng)編譯環(huán)境并支持OpenCL 1.2 標(biāo)準(zhǔn).

有向圖的頂點(diǎn)集大小n分別設(shè)為10、50、200、300、1 400、2 000、5 000、8 920 和10 000，從小規(guī)模頂點(diǎn)集到大規(guī)模頂點(diǎn)集都均勻覆蓋，有利于準(zhǔn)確反映全源對(duì)最短路徑并行算法的性能可擴(kuò)展性.rand()隨機(jī)數(shù)函數(shù)的隨機(jī)數(shù)種子采用n可生成一組隨機(jī)數(shù)，生成具有代表的節(jié)點(diǎn)數(shù)量為10、50、200、300、1 400、2 000、5 000、8 920 和10 000 的圖，并構(gòu)成相應(yīng)的鄰接矩陣A.鄰接矩陣是使用二維數(shù)組模擬線性代數(shù)中矩陣的結(jié)構(gòu).此種數(shù)據(jù)結(jié)構(gòu)所使用的儲(chǔ)存空間較大，但其使用較為便捷.根據(jù)重復(fù)平方法最短路徑算法的原理，分別基于OpenMP 和CUDA 并行計(jì)算架構(gòu)實(shí)現(xiàn)了最短路徑并行算法.

測試并記錄CPU_SP 系統(tǒng)（平臺(tái)1 測試），OMP_SP 系統(tǒng)（平臺(tái)1 測試）、CUDA_SP 系統(tǒng)、基于AMD GPU 的OCL_SP 系統(tǒng)和基于NVIDIA GPU的OCL_SP 系統(tǒng)的處理時(shí)間，如表2 所示.從表2可以看出，隨著頂點(diǎn)數(shù)的增大，最短路徑算法的執(zhí)行時(shí)間越來越長.基于GPU 的最短路徑算法運(yùn)行時(shí)間要比CPU_SP 和OMP_SP 算法的運(yùn)行時(shí)間少的多，而OMP_SP 算法運(yùn)行時(shí)間要比CPU_SP 算法的運(yùn)行時(shí)間更短.

表2 不同計(jì)算平臺(tái)下最短路徑算法執(zhí)行時(shí)間Tab.2 The shortest path algorithm execution time under different computing platforms

用加速比作為加速效果的衡量標(biāo)準(zhǔn)，可以驗(yàn)證各種架構(gòu)下并行算法的效率.加速比定義如下：

式中：Ts是指算法運(yùn)行在單線程CPU 平臺(tái)上的時(shí)間，Tp是指算法運(yùn)行在多線程CPU 或CPU+GPU協(xié)同計(jì)算平臺(tái)上的時(shí)間.

OMP_SP 并行算法執(zhí)行時(shí)間Tpo與基于NVIDIA GPU 的OCL_SP 并行算法執(zhí)行時(shí)間Tpnc的比值為相對(duì)加速比R1，定義為：

CUDA_SP 并行算法執(zhí)行時(shí)間Tpc與基于NVIDIA GPU 的OCL_SP 并行算法Tpnc執(zhí)行時(shí)間的比值為相對(duì)加速比R2，定義為：

式中：在CPU 和GPU 上OpenCL 內(nèi)核的計(jì)算時(shí)間之和為Tk，CPU 和GPU 之間數(shù)據(jù)傳輸時(shí)間之和為Th，系統(tǒng)初始化時(shí)間之和為To.

加速比可用于評(píng)估并行算法相比CPU 串行算法的性能提升.相對(duì)加速比R1可用于評(píng)估基于NVIDIA GPU 的OCL_SP 并行算法相比OMP_SP并行算法的性能提升.相對(duì)加速比R2可用于評(píng)估基于NVIDIA GPU 的OCL_SP 并行算法相比CUDA_SP 并行算法的性能提升.根據(jù)表2 可得算法加速比情況，如表3 所示.

表3 不同計(jì)算平臺(tái)下最短路徑并行算法性能對(duì)比Tab.3 Performance comparison of parallel algorithm of shortest path on different computing platforms

4.2 并行算法性能分析

4.2.1 系統(tǒng)性能瓶頸 OCL_SP 并行算法需要n×n×log(n-1)×2次右鄰接矩陣的存儲(chǔ)器讀取操作，以及n×n×log(n-1)次輸出矩陣的存儲(chǔ)器寫入操作.頂點(diǎn)集大小為n=10 000，每個(gè)數(shù)據(jù)占用2 B的存儲(chǔ)空間.因此，存儲(chǔ)器訪問的數(shù)據(jù)總量約為6 GB.除以設(shè)備上內(nèi)核執(zhí)行所用時(shí)間0.028 s.系統(tǒng)實(shí)際獲得的帶寬約為214.29 Gb/s，接近GeForce GTX1070 顯存的理論帶寬.由此可見，全局存儲(chǔ)器的帶寬限制了OCL_SP 并行算法效率的進(jìn)一步提高.因此，OCL_SP 并行算法的性能瓶頸是顯示內(nèi)存帶寬.

4.2.2 不同架構(gòu)下最短路徑算法處理時(shí)間分析 根據(jù)表2 顯示，對(duì)不同頂點(diǎn)數(shù)的加權(quán)有向圖進(jìn)行全源對(duì)最短路徑處理將在 4 種計(jì)算平臺(tái)上進(jìn)行.在頂點(diǎn)數(shù)一樣時(shí)，最短路徑算法在不同并行計(jì)算架構(gòu)下的處理時(shí)間相對(duì)串行處理時(shí)間均有不同程度縮減，即獲得了性能提升效果.如：對(duì)于頂點(diǎn)數(shù)為10 000的最短路徑串行運(yùn)算時(shí)間為407.10 s，OMP_SP 算法運(yùn)算時(shí)間縮短為76.09 s，CUDA_SP 算法運(yùn)算時(shí)間則大幅縮減為2.13 s，而基于OpenCL 并行計(jì)算平臺(tái)上則保持在2.07 s 左右的低運(yùn)算時(shí)間.

設(shè)n為加權(quán)有向圖G(V,E)的頂點(diǎn)數(shù)，c為多核CPU 核心數(shù)，OMP_SP 并行算法的時(shí)間復(fù)雜度為在相同頂點(diǎn)數(shù)下，存在OMP_SP 并行算法的時(shí)間復(fù)雜度因此，從表2 可以看出，OMP_SP 并行算法相比CPU_SP 串行算法的處理時(shí)間有明顯降低.同時(shí)可以看出，設(shè)置線程數(shù)c恒定時(shí)，隨著頂點(diǎn)數(shù)n的增大，OMP_SP并行算法的時(shí)間復(fù)雜度也隨之增加.每一個(gè)CPU線程負(fù)責(zé)處理的鄰接矩陣乘的數(shù)據(jù)規(guī)模也隨之增大，OMP_SP 并行算法處理時(shí)間也越來越長，基本符合線性增長的趨勢.

當(dāng)鄰接矩陣的計(jì)算規(guī)模較小時(shí)，在GPU 算法中啟動(dòng)的工作項(xiàng)計(jì)算負(fù)荷不足，系統(tǒng)在調(diào)度方面有較大時(shí)間開銷.此時(shí)GPU 算法與CPU 算法的處理時(shí)間接近，GPU 并行計(jì)算的優(yōu)勢沒有展現(xiàn).隨著鄰接矩陣的計(jì)算規(guī)模增大，每個(gè)工作項(xiàng)計(jì)算負(fù)荷也逐漸提高，而系統(tǒng)調(diào)度的時(shí)間占比降低，GPU 并行計(jì)算占比上升，速度提升明顯.

4.2.3 并行計(jì)算架構(gòu)下最短路徑算法加速效果對(duì)比分析 設(shè)n為加權(quán)有向圖G(V,E) 的頂點(diǎn)數(shù)，u為GPU 啟動(dòng)的工作項(xiàng)的數(shù)據(jù)量，CUDA_SP 和OCL_SP 并行算法的時(shí)間復(fù)雜度為由于在GPU 中可以維護(hù)大量的活動(dòng)線程和工作項(xiàng)，即u?c.因此，在相同頂點(diǎn)數(shù)下，存在CUDA_SP和OCL_SP 并行算法的時(shí)間復(fù)雜度從表3 可以看出，當(dāng)有向圖的頂點(diǎn)集合較小時(shí)，OMP_SP 并行算法與基于GPU 的最短路徑并行算法的加速效果相當(dāng).分析其原因，OpenMP的并行開銷要小于GPU.由于OMP_SP 并行算法是在CPU 端執(zhí)行，不同于GPU 并行算法需要在CPU 與GPU 之間進(jìn)行數(shù)據(jù)傳遞.而且OpenMP 啟動(dòng)的線程少，線程維護(hù)成本低.當(dāng)頂點(diǎn)集合逐步增大時(shí)，OMP_SP 并行算法的加速效果始終維持在一個(gè)較低的水平，而GPU 并行算法的加速效果更好.這主要是由于OMP_SP 并行算法需要操作系統(tǒng)將線程調(diào)度到對(duì)應(yīng)數(shù)目的CPU 核心上執(zhí)行，多核CPU 的核心數(shù)是有限的.反觀計(jì)算資源充裕的GPU 可以創(chuàng)建足夠的工作項(xiàng)以滿足大規(guī)模數(shù)據(jù)的并行處理.同時(shí)，隨著頂點(diǎn)數(shù)n的增大，CUDA_SP和OCL_SP 并行算法的時(shí)間復(fù)雜度也隨之增加，每一個(gè)線程或工作項(xiàng)負(fù)責(zé)處理的鄰接矩陣乘的數(shù)據(jù)量也隨之增大.從表2 可以看出，CUDA_SP 和OCL_SP 并行算法處理時(shí)間也適當(dāng)?shù)脑鲩L.

從表3 可以看出，當(dāng)頂點(diǎn)集合規(guī)模較小時(shí)，GPU 加速的最短路徑并行算法性能改善并不明顯.當(dāng)頂點(diǎn)集合規(guī)模較大時(shí)，CUDA_SP 和OCL_SP 并行算法的性能提升較為顯著，在3 種不同并行計(jì)算架構(gòu)下有效的性能收益明顯.當(dāng)頂點(diǎn)集合n=10 000時(shí)，基于NIVIDA GPU 的OCL_SP 并行算法的加速比更達(dá)到了196.35 倍.然而，當(dāng)頂點(diǎn)集合規(guī)模n＞5 000以后，GPU 加速的最短路徑算法加速比曲線呈現(xiàn)一種緩慢上升的過程.說明在頂點(diǎn)集規(guī)模較大的情況下，系統(tǒng)已不能獲得更快速的性能提高.主要原因是GPU 的計(jì)算資源接近飽和.主機(jī)和設(shè)備間的數(shù)據(jù)傳輸時(shí)間占比越來越大，并行計(jì)算時(shí)間的縮短已對(duì)加速比貢獻(xiàn)甚微.

采用離線編譯內(nèi)核讀寫數(shù)據(jù)文件方式的OCL_SP 并行算法，相比采用在線編譯內(nèi)核讀寫數(shù)據(jù)文件方式的CUDA_SP 并行算法減少了應(yīng)用初始化時(shí)間.如表3 所示，在同等數(shù)據(jù)集規(guī)模下，OCL_SP 并行算法的運(yùn)算耗時(shí)更少，與CUDA_SP并行算法性能相比略有提升，最大獲得了2.25 倍加速比.OCL_SP 并行算法性能與OMP_SP 并行算法性能相比則有很大的提高，最大獲得了36.76 倍加速比.

將本文提出的OCL_SP 算法的整體加速效果與文獻(xiàn)[25]進(jìn)行對(duì)比.由于其他文獻(xiàn)大多是利用最短路徑算法進(jìn)行各種應(yīng)用研究，很少專門針對(duì)全源對(duì)最短路徑算法進(jìn)行加速效果的研究.因此，無法直接進(jìn)行加速效果的對(duì)比.根據(jù)文獻(xiàn)[25]中提供的測試數(shù)據(jù)，當(dāng)每個(gè)核組開啟64 個(gè)從核時(shí)，經(jīng)過數(shù)組劃分優(yōu)化后的全源最短路徑并行算法獲得了106倍加速比.而根據(jù)表3 的測試結(jié)果可以看到，OCL_SP 并行算法最高獲得了196.35 倍加速比.因此，OCL_SP 并行算法相比文獻(xiàn)[25]中的算法取得了更好的加速性能.

4.2.4 OCL_SP 并行算法跨平臺(tái)性分析 源碼級(jí)別的算法可移植性不僅要求在不同的平臺(tái)上源碼可以成功地編譯和運(yùn)行，同時(shí)還需要算法保持相當(dāng)?shù)男阅?CUDA_SP 并行算法受到硬件平臺(tái)的限制，表3 中數(shù)據(jù)顯示，在多種硬件平臺(tái)上OCL_SP 并行算法具備了最大程度的兼容性和可移植性.

5 結(jié)束語

針對(duì)全源對(duì)最短路徑圖算法中數(shù)據(jù)量大、運(yùn)算量大和高計(jì)算密集度等問題，本文設(shè)計(jì)了采用GPU 加速的最短路徑并行算法.分別從算法可并行性分析、算法描述和框架的設(shè)計(jì)、矩陣乘法的加速策略以及存儲(chǔ)器優(yōu)化策略等方面進(jìn)行了闡述.實(shí)驗(yàn)結(jié)果表明，OCL_SP 并行算法分別與CPU_SP 串行算法、OMP_SP 并行算法和CUDA_SP 并行算法的性能相比，取得了196.35、36.76 和2.25 倍的加速比，不但實(shí)現(xiàn)了高性能，而且在不同GPU 計(jì)算平臺(tái)間實(shí)現(xiàn)了性能移植.擁有大量計(jì)算單元和強(qiáng)大計(jì)算能力的GPU，為圖論有關(guān)問題的快速求解提供了一種廉價(jià)而高效的計(jì)算平臺(tái)，為高效地大規(guī)模數(shù)值計(jì)算進(jìn)行了有益探索.在下面的工作中，我們將把多CPU 與多GPU 構(gòu)成的異構(gòu)計(jì)算平臺(tái)作為并行計(jì)算模式，采用并行I/O 技術(shù)，以達(dá)到更大規(guī)模數(shù)據(jù)集對(duì)計(jì)算能力的要求.