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

因為計畫的關係，需要找出果蠅腦神經元 (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，所以相當於是一個 12529x12529 大小的 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\n", elapsedTime);
 39:  printf("time: %.3lf\n", 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 卡，一次可以把整個 12529x12529 的 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: 

第一支 CUDA 程式

CUDA 已經熾手可熱一段時間了，好像程式經過CUDA加持，效能都加了不少。可惜CUDA不能加在人身上，不然原本一百公尺跑 個15秒，CUDA上身也許可以跑個9秒9。

我也 觀望一段時間了，想看看CUDA和OpenCL誰勝出，再來學。但等著等著 ，火都快熄了。上網查了一下，看到有一本新書，"CUDA By Example"，在台灣買要新台幣1800元，去 Amazon.com 查了一下，含運費（ 最一般的運費，號稱要40天才能到）才美金不到五十元。心想也 不急，就訂了。十天左右，就拿到書。翻了前幾章，有點基本概念， 剛好有一筆資料要處理，雖然已經先用傳統方式解決了，心裡覺得可 以用CUDA試試，就動手來寫第一支CUDA程式了。

資料來源是模擬的結果，想要看 Iso-surface。但我一開 始只得到一筆200x300的網格資料，如下圖：

那是對一個圓柱形流場做模擬，所以雖然只給一片資料，但因 為是對稱資料，得先由那筆資料算出整個圓柱體的資料，示意圖如下 ：











想法很簡單，就是建立一個三維的volume data 把這圓柱包 住。資料的格式，是200x300，所以每讀進200的點，就可以產生出一 個 slice 的資料(如下圖）

如上圖右，這 slice 的每個點，都可以由和原點的距離來決定。Pseudo Code 如下：

---

CreateOneSlice(float *line, int r)
{
  for (i=0 to 2*r-1)
 for (j=0 to 2*r-1)
 {
      dist = Dist(O, Pij)
      Vij = Interpolation(dist, line)
    }
}

---

如此重覆300次，就完成了整個Volume data了。

CUDA Implementation

因為一面在看 CUDA 的書，就想到，在單一 Slice 上的每個點，都是獨立的，可以分別用一個 thread 來計算，很適合平行計算，就決定用 CUDA 來試試。

---

void CreateOneSliceGPU(int attr_index, FILE *outf, float *dev_Slice)
{
  int sliceDim = g_DIMY*2-1;

  if(g_Slice==0)    g_Slice = new float[sliceDim*sliceDim];

  double t1, t2;
t1 = wallclock();

  cudaMemcpy(dev_fYData, g_fYData, g_DIMY*sizeof(float), cudaMemcpyHostToDevice);

  dim3 grid((sliceDim+block_dim-1)/block_dim,
      (sliceDim+block_dim-1)/block_dim);
  dim3 threads(block_dim, block_dim);

  slice_kernel<<<grid, threads>>>(g_DIMY, dev_fYData, dev_Slice);

  cudaMemcpy(g_Slice, dev_Slice, sliceDim*sliceDim*sizeof(float), cudaMemcpyDeviceToHost);

t2 = wallclock();

TIME_USED += (t2-t1);

fwrite(g_Slice, sizeof(float), sliceDim*sliceDim, outf)
}

---

dev_fYData 就儲存200個點的值，做為 linear-interpolation用。每讀入200點資料就可以輸出一張Slice。dev_Slice 已事先用 cudaMalloc 宣告一塊記憶體了。

Kernel code 如下：

---

__global__ void slice_kernel(int dimy, float *dev_fYData, float* dev_Slice)
{
int x = threadIdx.x + blockIdx.x*blockDim.x;
int y = threadIdx.y + blockIdx.y*blockDim.y;
int offset;
float dx2, dy2;
float r = (float)(dimy)-1;
int sliceDim = 2*dimy-1;
float dist;
int lindex, uindex;

float vl, vu;

offset = x+y*sliceDim;
if(x < sliceDim && y < sliceDim)
{
 dx2 = (r-x) * (r-x);
 dy2 = (r-y) * (r-y);
 dist = sqrtf(dx2+dy2);
 if(dist>r)
 {
     dev_Slice[offset] = 0.f;
 }
 else
 {
     lindex = (int)dist;
     if(lindex==r)
     {
         dev_Slice[offset] = dev_fYData[dimy-1];
     }
     else
     {
         uindex = lindex+1;;
         vl = dev_fYData[lindex];
         vu = dev_fYData[uindex];
         dev_Slice[offset] = vl+(dist-(float)lindex)*(vu-vl);
     }
 }
}
else
{
 dev_Slice[offset] = 0.f;
}
}

---

結果：

後來真的實做，是用兩筆資料，150x300 和 300x600 兩組。block_dim 我有做不同的測試，以下是其中一組結果(block_dim=10)：

CPU vs. GPU

150x300 : 1086 ms vs. 195 ms

300x600 : 9030 ms vs. 1276 ms

運算時間的確縮短不少。

後記：

1. 我後來還有再把 g_fYData 放在 constant memory 裡，時間還會再縮短一些。

2. 這支程式，其實是一個多月前寫的。這幾個星期，又再試了 APSP (All Pair Shortest Path)的 CUDA code，有空再放上來(FW 方法，和 blocked 加速方法)。

3. 本次結果的幾張 snapshots:

圓柱的三個方向切面，加上 Contour Lines.

Iso-Surfaces

另一組 Iso-Surfaces