因為計畫的關係,需要找出果蠅腦神經元 (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 block4: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 block19: {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 block4: 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_columnBlocks22: */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 secondsNeuron 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_columnBlocks3: *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 end18: 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
想問個題外問題,在cuda中建立struct with pointer 的data member。先看看我的程式碼:http://codepad.org/dUYFOKL6我現在在學習如何建立一個指標指向的struct,並正確的配置記憶體在device端。但是卡住在如何驗證我的結果是我要的…
我試了一下, 有試出一個解決問題的方法。我把改好的程式碼放在以下連結:http://codepad.org/UspyNSK6和你原本程式碼最大的不同, 是在於 device 那個 struct 變數的宣告。
非常感謝您的回覆。另外也附上我今天嘗試的做法..不過仍是不正確,請看看我嘗試的作法:http://codepad.org/C2bWIkBc不知道為什麼最後p的內容竟然沒有改變?
我才注意到我之前的 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], 我覺得用處並不大, 有興趣的話就做個參考。如果你找到解法, 請通知一下, 我也很有興趣想知道怎麼做。謝謝..
非常高興你對這個議題也有興趣。目前我也還沒試出一個比較好的方法,另外這是我的mail,youknowme09@hotmail.com如果有任何發現也可以寄信跟我討論!! 非常感謝!!
那是因為host 跟 device的memory的address是不一樣的, pointer紀錄的是address所以傳到另一個地方就會存取到錯誤的address了
不好意思請問一下:
1. 上方的FW_APSP()內的threadNum值是設為多少,因為我發現,不同的thread個數,所能接受執行個數的範圍好像會有不同的差異。
2. 對於APSP_kernel()內的dev_pre_matrix有什麼作用嗎?這個對於dev_matrix好像沒有資料相依性的作用,而且對於大型數據能確定是哪個block先執行嗎?會不會有資料相依性,會不會造成的race condition的問題?
謝謝!
1. threadNum 我是設16, 我當時有設不同的大小, 後來用16。
2. 那個 dev_pre_matrix 是要儲存路徑 predecessor 的matrix, 當算出 shortest path 後, 要重建 path 的話, 就會用到。如果純粹只要算出 shortest path 的長度, 就不用。
3. 有三個 phase, 只要三個 phase 是按照順序做, 每個 phase 裡 block 執行順序並不會影響。
感謝你的回覆!
對於單純的Folyd warshall apsp我已經測了好幾天了,我嘗試著套用你的結構,卻只能跑到最多250個node,不好意思請教你,除了上述的結構外,請問還有什麼其它的設定或者定義嗎?
謝謝!
Hi,
我查了一下, 並沒有什麼特殊設定。請問你用什麼卡呢?
真的很不好意思!
後來才發現我們實驗室的卡應該是太過於老舊且是一般的卡,可能裏面有一兩個單位有問題,所以才造成後續運算整個出錯,後來換了一張工作用的,數據才正確。
我們目前正在嘗試將此程式以shared memory的方式改得更快一些,不知是否可以跟您請教完整的原始碼?
謝謝!
freedom_chu@yahoo.com.tw