用 CUDA 來解 All Pair Shortest Path (APSP) 問題

2010-11-292010-11-29| cylcyl| 11 Comments| 15:17
Categories:

因為計畫的關係，需要找出果蠅腦神經元 (neuron) 之間的可能路徑。找出路徑，可以方便了解在果蠅腦內訊息的傳遞。

前置工作有很多，包含前處理，資料庫建立，這資料庫是先把三維腦神經影像裡的空間資料做整理而產生的（是另一個同事的辛苦成果），再由資料庫一筆一筆Query，建立出一個 connection matrix。這 connection matrix 就相當於是一個 weight 的 Matrix, 如果 i, j 這兩點有相連，m(i, j) = 1；不然就等於一個很大的數，代表沒有直接連結。

我的想法是直接找出 APSP，接下來用查表的就可以找出路徑。算出APSP最簡單的方式，應該是 Floyd Warshall 方法吧。

procedure FloydWarshall ()
for k := 1 to n
for i := 1 to n
for j := 1 to n
path[i][j] = min ( path[i][j], path[i][k] path[k][j] );

根據 Pawan Harish 和 P.J. Narayanan 在2007年發表的 paper, “Accelerating large graph algorithms on the GPU using CUDA” 裡所提到的演算法：

圖1

實際用 CUDA 來寫，程式碼如下：

  1: void FW_APSP()


  2: {


  3: 	int *dev_matrix;


  4: 	size_t msize = matrix_dim*matrix_dim*sizeof(int);


  5: 


  6: 	cudaMalloc((void**)&dev_matrix, msize);


  7: 	cudaMalloc((void**)&dev_pre_matrix, msize) ;


  8: 


  9: 	/* apsp_matrix is pre-loaded */


 10: 	cudaMemcpy(dev_matrix, apsp_matrix, msize, cudaMemcpyHostToDevice);


 11: 	cudaMemcpy(dev_pre_matrix, pre_matrix, msize, cudaMemcpyHostToDevice);


 12: 	dim3 block((matrix_dim threadNum-1)/threadNum, (matrix_dim threadNum-1)/threadNum);


 13: 	dim3 threads(threadNum, threadNum);


 14: 	


 15: 	for(int k=0; k<matrix_dim; k  )


 16: 	{


 17: 		APSP_kernel<<<block, threads>>>(k, matrix_dim, dev_matrix, dev_pre_matrix);


 18: 	}


 19: 


 20: 	cudaMemcpy(apsp_matrix, dev_matrix, msize, cudaMemcpyDeviceToHost);


 21: 	cudaMemcpy(pre_matrix, dev_pre_matrix, msize, cudaMemcpyDeviceToHost);


 22: 


 23: 	cudaFree(dev_matrix);


 24: }


Kernel 的地方也很容易：


  1: __global__ void APSP_kernel(int k, int matrix_dim, int *dev_matrix, int *dev_pre_matrix)


  2: {


  3: 	int ioffset, koffset;


  4: 	int Dij, Dik, Dkj;


  5: 	


  6: 	int i = blockIdx.x*blockDim.x threadIdx.x;


  7: 	int j = blockIdx.y*blockDim.y threadIdx.y;


  8: 


  9: 	if(i<matrix_dim && j<matrix_dim)


 10: 	{


 11: 		koffset = k*matrix_dim;


 12: 		ioffset = i*matrix_dim;


 13: 		Dij = dev_matrix[ioffset j];


 14: 		Dik = dev_matrix[ioffset k];


 15: 		Dkj = dev_matrix[koffset j];


 16: 		if(Dik Dkj<Dij)


 17: 		{


 18: 			dev_matrix[ioffset j] = Dik Dkj;


 19: 			dev_pre_matrix[ioffset j] = k;


 20: 		}


 21: 	}


 22: }


我去網路上找了一組 test data：rome99.txt（原始說明，請見此網頁）。這筆資料有 3353 個 nodes，以及 8870 個 edges。用上面的程式來跑，一點問題也沒有。改成用我實際的神經元資料，有 12529 個 neuron，所以相當於是一個 12529×12529 大小的 matrix。一跑，程式馬上就跳出來，檢查後才發現是在 cudaMalloc 那裡就死了。因為要顯卡上 malloc 一塊 12529x12529x4 bytes 的記憶體，我的顯卡是 GeForce 9800 GX2，才只有 512 MB 的記憶體。當然沒法執行。


上網搜尋相關文件，發現 2008 年有一篇論文，G. J. Katz 和 J. T. Kider Jr 所的 "All-Pairs Shortest-Paths for Large Graphs on the GPU” 裡，提到一種方法，以 block 的方式，來計算 APSP 的問題。這篇文章裡的方法是來自 G. Venkaaraman，Sartaj Sahni，和 S. Mukhopadhyaya 在2003 年的 paper："A Blocked All-Pairs Shortest Paths Algorithm”，2003 年的 paper，原本是要善用電腦的 cache 來加速。2008 年的作者就改成用 CUDA 來計算，效果不錯；最重要的，是可以一次只載入一部份，所以就算原始資料超過記憶體大小，還是可以把原始資料切成小 block 來做運算。下圖是此演算法的示意圖：




 圖2


先把整個 Matrix 切成比較小的 block。然後每一個在對角線的 block，都有三個步驟 (phases)。圖2 上面的二個圖，就是三個 phase 的示意圖。圖2上左，是第一個 phase，先計算對角線那個 block（稱為 primary block）。接下來第二個 phase(圖2上中)，是計算在和 primary block 同樣 column 和 row 的 blocks。第三個 phase，就是計算其餘的 blocks。把對角線的的 primary blocks 都重覆以上三步驟，就可以計算出 APSP matrix。三個 phases 的分別解說如下：


第一個 phase，只是一般的 APSP，可以直接用 FW 方法來做。


第二個 phase，需要用到第一個 phase 算出來的結果：


 圖3


第三個 phase，需要用到第二個 phase 所計算出來的 blocks：


 圖4


如圖5，白色粗線框起來的地方，是 Phase 3 要計算的一個 block，它需要那兩個黑線粗框的 block 裡的資料才能計算出結果。


 圖5




Blocked_APSP_GPU 程式碼如下：


  1: void Blocked_APSP_GPU()


  2: {


  3: 	int pStart, pEnd;	// primary block


  4: 	


  5: 	int block_size;


  6: 	clock_t startTime, endTime;


  7: 	double elapsedTime;


  8: 


  9: 	NUM_PBLOCK = (g_MATRIX_DIM g_BLOCK_DIM-1)/g_BLOCK_DIM;


 10: 


 11: 	RowBlocks = new int [g_BLOCK_DIM*g_MATRIX_DIM]; 


 12: 	HANDLE_ERROR(cudaMalloc((void**)&dev_lineBlocks, g_BLOCK_DIM*g_MATRIX_DIM*sizeof(int)));


 13: 	HANDLE_ERROR(cudaMalloc((void**)&dev_columnBlocks, g_BLOCK_DIM*g_MATRIX_DIM*sizeof(int)));


 14: 


 15: 	startTime = clock();


 16: 	t1=wallclock();


 17: 


 18: 	for(int pb =0; pb <NUM_PBLOCK; pb   )		// primary block


 19: 	{


 20: 		pStart = pb*g_BLOCK_DIM; 


 21: 		pEnd = pStart g_BLOCK_DIM-1;


 22: 		if(pEnd>=g_MATRIX_DIM)


 23: 		{


 24: 			pEnd = g_MATRIX_DIM-1;


 25: 		}


 26: 		


 27: 		PhaseI(pStart, pEnd);


 28: 		PhaseII_GPU(pb, pStart, pEnd);


 29: 		PhaseIII_GPU(pb, pStart, pEnd);


 30: 	}


 31: 	endTime = clock();


 32: 	t2 = wallclock();


 33: 


 34: 	HANDLE_ERROR(cudaFree(dev_lineBlocks));


 35: 	HANDLE_ERROR(cudaFree(dev_columnBlocks));


 36: 


 37: 	elapsedTime = (double(endTime-startTime)/CLOCKS_PER_SEC);


 38: 	printf("Blocked APSP: %lf seconds", elapsedTime);


 39: 	printf("time: %.3lf", t2-t1);


 40: 	


 41: }


因為顯示卡記憶體有限，所以原始 matrix 先切成小的 blocks。在 PhaseII 裡，每次在顯卡上，我先載入並計算位於和 primary block 相同 column 的blocks，再載入和計算位於和 primary block 相同 row 的 blocks。在上面列表 Line 12、Line 13 裡的 g_MATRIX_DIM 就是指整個 matrix 的 dimension (以本例來說，是12529)。g_BLOCK_DIM 是自訂的參數，以本例來說，我是設 32。


Phase II 


在Phase II 裡，顯卡上宣告的 memory 大小為 32x12529x4 bytes。當程式在載入和 primary block 相同 column (或 row) 的 blocks 時，連 primary block 也會被載入(下面程式碼 line 12)，所以可以算出結果，沒有問題。在 PhaseII 裡資料是被放在 dev_columnBlocks 這個空間。Line 18 – Line 21 是計算和 primary block 相同 column 的 blocks； Line 30-35 是計算和 primary block 相同 row 的 blocks。Line 27 是一個函式，先把整個 row 的 blocks 複製到一塊連續記憶體 (RowBlocks)，Line 28 再從 RowBlocks 複製到顯示卡記憶體裡。


  1: void PhaseII_GPU(int bid, int pStart, int pEnd)


  2: {


  3: 	int cStart, cEnd;	// current block


  4: 	int dim_i;


  5: 	int Dij, Dik, Dkj;


  6: 	int ioffset;


  7: 


  8: 


  9: 	int pdim = pEnd-pStart 1;


 10: 


 11: 	// copy blocks of column bid to to device


 12: 	HANDLE_ERROR(cudaMemcpy(dev_columnBlocks, &(apsp_matrix[pStart*g_MATRIX_DIM]),  


 13: 					sizeof(int) * pdim * g_MATRIX_DIM, cudaMemcpyHostToDevice));


 14: 


 15: 	dim3 block(NUM_PBLOCK);   //(i, *) and (*, j)


 16: 	dim3 threads(g_BLOCK_DIM);


 17: 


 18: 	for(int k=0; k<pdim; k  )		


 19: 	{


 20: 		PhaseII_Column_kernel<<<block, threads>>>(k, bid, pStart, pEnd, g_MATRIX_DIM, dev_columnBlocks);	


 21: 	}


 22: 


 23: 	// copy values from device memory back to matrix


 24: 	HANDLE_ERROR(cudaMemcpy(&(apsp_matrix[pStart*g_MATRIX_DIM]), dev_columnBlocks, 


 25: 			sizeof(int) * pdim * g_MATRIX_DIM, cudaMemcpyDeviceToHost));


 26: 


 27: 	CopyToRowBlocks(pStart, pEnd);


 28: 	HANDLE_ERROR(cudaMemcpy(dev_columnBlocks, RowBlocks,  


 29: 					sizeof(int) * pdim * g_MATRIX_DIM, cudaMemcpyHostToDevice));


 30: 	for(int k=0; k<pdim; k  )		


 31: 	{


 32: 		PhaseII_Row_kernel<<<block, threads>>>(k, bid, pStart, pEnd, g_MATRIX_DIM, 


 33: 													dev_columnBlocks);	


 34: 	}


 35: 	HANDLE_ERROR(cudaMemcpy(RowBlocks, dev_columnBlocks,  


 36: 					sizeof(int) * pdim * g_MATRIX_DIM, cudaMemcpyDeviceToHost));


 37: 	CopyFromRowBlocks(pStart, pEnd);


 38: 


 39: 


 40: }


Phase III


但在 PhaseIII 裡，需要兩個 input block 才能算出一個 output block。我就在 Phase III 一開始時，先把和 primary block 相同 column 的 blocks 先載入 dev_lineBlocks (下表 line 10 ，line 11)；然後再其他blocks，一次 copy 一整個  column 的 blocks 到 dev_columnBlocks，計算結果，再複製回來(下圖 line 16- line 38)。


  1: void PhaseIII_GPU(int pbid, int pStart, int pEnd)


  2: {


  3: 


  4: 	int iStart;


  5: 	int iEnd;


  6: 	int pdim = pEnd-pStart 1;


  7: 	int idim;


  8: 


  9: 	// copy blocks of column pbid to dev_lineBlocks


 10: 	HANDLE_ERROR(cudaMemcpy(dev_lineBlocks, &(apsp_matrix[pStart*g_MATRIX_DIM]),  


 11: 							sizeof(int) * pdim * g_MATRIX_DIM, cudaMemcpyHostToDevice));


 12: 


 13: 	dim3 block(NUM_PBLOCK);   


 14: 	dim3 threads(g_BLOCK_DIM);


 15: 


 16: 	for(int I=0; I<NUM_PBLOCK; I  )


 17: 	{


 18: 		if(I==pbid) continue;


 19: 


 20: 		/* 


 21: 		    copy blocks of Column I to dev_columnBlocks 


 22: 		*/


 23: 		iStart = g_BLOCK_DIM*I;


 24: 		iEnd = iStart g_BLOCK_DIM-1; 	


 25: 		if(iEnd>=g_MATRIX_DIM) iEnd = g_MATRIX_DIM - 1;


 26: 


 27: 		idim = iEnd-iStart 1;


 28: 		HANDLE_ERROR(cudaMemcpy(dev_columnBlocks, &(apsp_matrix[iStart*g_MATRIX_DIM]),  


 29: 							sizeof(int) * idim * g_MATRIX_DIM, cudaMemcpyHostToDevice));


 30: 		// envoke  kernel (pbid, j),  0<=j<NUM_PBLOCK, j!=bid)


 31: 		for(int k=0; k<pdim; k  )


 32: 			PhaseIII_kernel<<<NUM_PBLOCK, threads>>>(pbid, k, iStart, iEnd, pStart, pEnd,


 33: 									g_MATRIX_DIM, dev_lineBlocks, dev_columnBlocks);


 34: 		


 35: 		// copy the values in dev_columnBlocks back to matrix


 36: 		HANDLE_ERROR(cudaMemcpy(&(apsp_matrix[iStart*g_MATRIX_DIM]), dev_columnBlocks, 


 37: 			sizeof(int) * idim * g_MATRIX_DIM, cudaMemcpyDeviceToHost));


 38: 	}


 39: }


結果：


後來我有申請 GPU Cluster 帳號，並在上面做測試。GPU Cluster 的每台機器，都是裝 4G 的 Tesla T10 卡，一次可以把整個 12529×12529 的 matrix 載入也沒問題。


rome99.txt:




Sequential FW         61.699 seconds

FW_GPU                 60.607 seconds

block_APSP_GPU   11.875 seconds




Neuron connection table:




Sequential FW          ~3100 seconds

FW_GPU                  3452.3 seconds

block_APSP_GPU    500.1 seconds




以GPU 來計算 block_APSP，時間只要約1/6。成效很顯著。


後記：


這程式麻煩的地方之一，是在於 matrix 的 dimension 無法剛好是 block dimension 的倍數。解決最後一個 primary block 的問題，讓我腦筋迷糊了好一會。Kernel 的程式碼，我附上其中之一，其餘的也類似。


  1: /**


  2:  *	a. copy blocks (I, *) to dev_columnBlocks


  3:  *


  4:  *  b. dev_lineBlocks:


  5:  *     ------


  6:  *     |I, *|


  7:  *     ------


  8:  *     |I, *|


  9:  *     ------


 10:  *     |I, *|


 11:  *     ------


 12:  *     |I, *|


 13:  *		....


 14:  */


 15: __global__ void PhaseII_Column_kernel(int k, int bid, int pStart, int pEnd, int matrix_dim, int *dev_Blocks)


 16: {


 17: 	int cStart, cEnd;		// current block start, current block end


 18: 	int pdim, dim, i, j;


 19: 	int ioffset, koffset;


 20: 	int Dij, Dik, Dkj;


 21: 


 22: 


 23: 	if(blockIdx.x==bid) return;


 24: 


 25: 	pdim = pEnd-pStart 1;


 26: 


 27: 	cStart = blockDim.x*blockIdx.x;


 28: 	cEnd = cStart blockDim.x-1;


 29: 	if(cEnd>=matrix_dim)


 30: 	{


 31: 		cEnd = matrix_dim-1;


 32: 	}


 33: 	dim = cEnd-cStart  1;


 34: 


 35: 	j = cStart threadIdx.x;


 36: 	koffset = k*matrix_dim;


 37: 


 38: 	if(j<matrix_dim)


 39: 	{


 40: 		//j = threadIdx.x;


 41: 		for(i=0; i<pdim; i  )


 42: 		{


 43: 			ioffset = i*matrix_dim; 


 44: 			Dij = dev_Blocks[ioffset j]; 


 45: 			Dik = dev_Blocks[ioffset pStart k];


 46: 			Dkj = dev_Blocks[koffset j];


 47: 			if(Dik Dkj<Dij)


 48: 			{


 49: 				dev_Blocks[ioffset j] = Dik Dkj;


 50: 			}


 51: 		}


 52: 		__syncthreads();


 53: 	}


 54: }


 55: 




本文同步發佈於：http://a-small-place.blogspot.com/2010/11/cuda-all-pair-shortest-path-apsp.html

11 thoughts on “用 CUDA 來解 All Pair Shortest Path (APSP) 問題”

  1. 想問個題外問題，在cuda中建立struct with pointer 的data member。先看看我的程式碼：http://codepad.org/dUYFOKL6我現在在學習如何建立一個指標指向的struct，並正確的配置記憶體在device端。但是卡住在如何驗證我的結果是我要的…

    回覆

  4. 我才注意到我之前的 code 有個錯誤, 我之前的 code 雖然可以成功的把 hp->data 和 dp->data互相交換, 但沒法把 hp 複製到 dp (反之亦然). 所以 hp->size 這欄並沒有被更新。我想這問題應是在於用 cudaMalloc 一塊記憶體時, 應就只能在 device kernel 中被 access。如下兩行程式： My_data *dp; cudaMalloc( (void**) &dp, sizeof(My_data));雖沒問題, 但是那樣一宣告, dp->data 反而沒法在 CPU 下 access。想在 CPU 下 用 cudaMalloc 來指定 dp->data 就會有錯誤。我們可以在 CPU 下用 cudaMalloc 宣告一塊記憶體, 但沒法在CPU裡把那塊記憶體指定給 dp->data, 因為 dp->data 只能在 GPU 裡才可以存取。我想到的作法是可以把 p 傳入 kernel, 然後在 kernel 裡指定 ip->data = p。我試了是沒問題, 但是這樣一來, 宣告 struct 就沒什麼意思了。 struct 裡每個變數都得另外傳入的話, 那就不需要用 struct 了。 我上網搜尋一下, 也找不到其他解決方法。我新改的程式碼在這裡 [urlhttp://codepad.org/MG272POs[/url], 我覺得用處並不大, 有興趣的話就做個參考。如果你找到解法, 請通知一下, 我也很有興趣想知道怎麼做。謝謝..

    回覆

  5. 非常高興你對這個議題也有興趣。目前我也還沒試出一個比較好的方法，另外這是我的mail,youknowme09@hotmail.com如果有任何發現也可以寄信跟我討論!! 非常感謝!!

    回覆

  6. 那是因為host 跟 device的memory的address是不一樣的, pointer紀錄的是address所以傳到另一個地方就會存取到錯誤的address了

    回覆

  7. 不好意思請問一下：
    1. 上方的FW_APSP()內的threadNum值是設為多少，因為我發現，不同的thread個數，所能接受執行個數的範圍好像會有不同的差異。
    2. 對於APSP_kernel()內的dev_pre_matrix有什麼作用嗎？這個對於dev_matrix好像沒有資料相依性的作用，而且對於大型數據能確定是哪個block先執行嗎？會不會有資料相依性，會不會造成的race condition的問題？
    謝謝！

    回覆

  8. 1. threadNum 我是設16, 我當時有設不同的大小, 後來用16。
    2. 那個 dev_pre_matrix 是要儲存路徑 predecessor 的matrix, 當算出 shortest path 後, 要重建 path 的話, 就會用到。如果純粹只要算出 shortest path 的長度, 就不用。
    3. 有三個 phase, 只要三個 phase 是按照順序做, 每個 phase 裡 block 執行順序並不會影響。

    回覆

  9. 感謝你的回覆！
    對於單純的Folyd warshall apsp我已經測了好幾天了，我嘗試著套用你的結構，卻只能跑到最多250個node，不好意思請教你，除了上述的結構外，請問還有什麼其它的設定或者定義嗎？
    謝謝！

    回覆

  10. Hi,

    我查了一下, 並沒有什麼特殊設定。請問你用什麼卡呢?

    回覆

  11. 真的很不好意思！
    後來才發現我們實驗室的卡應該是太過於老舊且是一般的卡，可能裏面有一兩個單位有問題，所以才造成後續運算整個出錯，後來換了一張工作用的，數據才正確。
    我們目前正在嘗試將此程式以shared memory的方式改得更快一些，不知是否可以跟您請教完整的原始碼？
    謝謝！
    freedom_chu@yahoo.com.tw

    回覆

Leave a Reply

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *

Related Posts

GitLab 簡單的 C++ 專案腳本範例

之前寫過《Gitlab CI/CD 簡單介紹》，大概介紹過 GitLab CI/CD 的架構了，而 Heres […] ...

Read MoreRead More

使用 QGraphicsScene 繪製 widget、產生 OpenGL Texture

從去年開始，Heresy 就有在試著用 Qt 搭配 OpenVR，開發支援 SteamVR 的虛擬實境的程式。 […] ...

Read MoreRead More

使用 Boost 的Locale 進行字碼轉換

這一篇會比較簡短一點，來大概提一下怎麼用 Boost C Libraries 的 Locale（官方文件）這個 […] ...

Read MoreRead More