如何編寫一個有兩個階段的高效CUDA程序

Question

我想編寫一個輕量級PIC（顆粒單元）程序。通過「輕量級」，我的意思是不需要擴大規模：假設所有數據都可以同時適用於單個GPU設備的內存和主機系統的內存。不過，我希望它儘可能快。

問題是，PIC的典型結構是兩個階段的集成：場解算器和粒子推進器。工作流程如下： 初始化系統 - >推粒子 - >求解域 - >推粒子 - >求解域... - >輸出

下一個推粒子或求解域必須等到前一個求解域或推粒子完成。可能需要數百萬次迭代才能獲得最終輸出。

作爲測試，省略場解算器，顆粒推動器可被寫爲：

__device__ 
void push(Particle &par) { 
    // some routines to move a particle. same excecutiong time for every particle. 
} 

並使用kernel_1像這樣excecute它：

__global__ 
void kernel_1(int n, Particle* parlist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    if (i < n) { 
     push(parlist[i]); 
    } 
} 

在主循環

for (int i=0;i<M;i++) { 
    kernel_1<<<(n+255)/256, 256>>>(n, parlist); 
} 

M是需要的迭代次數。但是，性能非常慢：在採用八核英特爾E5-2640 v3和Nvidia Quadro m4000的系統上，CUDA使用openmp提供了與純CPU版本類似的性能。對於10,000,000的粒子數和M = 1000，大約需要10秒。

但是如果我移動循環到內核：

void kernel_2(int n, Particle* parlist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    if (i < n) { 
     for (int i=0;i<M;i++) { 
      push(parlist[i]); 
     } 
    } 
} 

和

kernel_2<<<(n+255)/256, 256>>>(n, parlist); 

出於同樣的M = 1000，只需要100毫秒，這是一個10000％的加速。我已經驗證兩種情況下的結果是相同的和正確的。也許M次運行內核的調用成本太高。

將循環移入內核的性能改進非常令人難以置信。對於第一種情況，可以很容易地添加字段求解器：只需編寫一個新的內核並在主循環中順序執行兩個內核。然而，表現應該是藥物。

我發現很難在第二種情況下添加字段解算器例程：在沒有多次調用內核的情況下，塊之間似乎沒有同步機制，但字段解析器必須等待，直到所有粒子被推入，必須將其分配到不同的塊（因爲粒子的數量非常高）。

那麼是否有可能在一個內核中實現兩階段迭代？性能增益太多，不容忽視。

編輯： 我發現的性能差異非常混亂：100毫秒和10秒的差別僅僅是一個單一的代碼行或循環的甚至序列。我已修改推（）是一個小更復雜的（二維鮑里斯推杆）：他們創建上

class Particle 
{ 
public: 
    float x, y;  //m 
    float vx, vy;  //m/s 
    float m;   //kg 
    float q;   //ee 
}; 

__device__ 
void run(Particle& par, float B) 
{ 
    float t, s, vpx, vpy; 
    t = (par.q*ee*B/par.m)*dt/2; 
    s = 2*t/(1+t*t); 
    vpx = par.vx+t*par.vy; 
    vpy = par.vy-t*par.vx; 
    par.vx += s*vpy; 
    par.vy -= s*vpx; 
    par.x += par.vx*dt; 
    par.y += par.vy*dt; 
} 

我創建1 n元素陣列的粒子和1 n元素浮子陣列B.主機和cudaMemcpy到設備。然後我檢查了以下三個內核的性能：

__global__ 
void kernel_A(int n, int m, Particle* parlist, float* Blist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    int j; 
    if (i<n) { 
     for (j=0;j<m;j++) { 
      run(parlist[i], Blist[i]); 
     } 
    } 
} 

__global__ 
void kernel_B(int n, int m, Particle* parlist, float* Blist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    int j; 
    float B; 
    if (i<n) { 
     B = Blist[i]; 
     for (j=0;j<m;j++) { 
      run(parlist[i], B); 
     } 
    } 
} 

__global__ 
void kernel_C(int n, int m, Particle* parlist, float* Blist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    int j; 
    float B; 
    if (i<n) { 
     B = Blist[i]; 
     for (j=0;j<m;j++) { 
      run(parlist[i], B); 
      __syncthreads(); 
     } 
    } 
} 

__global__  
void kernel_D(int n, int m, Particle* parlist, float* Blist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    int j; 
    float B; 
    if (i<n) { 
     B = Blist[i]; 
    } 
    for (j=0;j<m;j++) { 
     if (i<n) { 
      run(parlist[i], B); 
     } 
    } 
} 

__global__ 
void kernel_E(int n, int m, Particle* parlist, float* Blist) 
{ 
    int i = blockIdx.x*blockDim.x + threadIdx.x; 
    int j; 
    float B; 
    if (i<n) { 
     for (j=0;j<m;j++) { 
      run(parlist[i], Blist[i]); 
      __syncthreads(); 
     } 
    } 
} 

而且運行時間有很大不同。對於n = 10,000,000且m = 1000：

• Kernel_A：7.6s
• Kernel_B：66ms
• Kernel_C：9.9S
• Kernel_D：10.0S
• Kernel_E：10.0S

這三個內核的結果完全相同並且正確（根據CPU版本進行檢查）。

我從官方的CUDA編程指南瞭解到，分支比較昂貴，所以kernel_C應該比kernel_B慢，但我懷疑差異是兩個數量級。我不明白的是爲什麼kernel_B比kernel_A執行得更好。 Kernel_B在kernel_A不需要訪問Blist 1000次，但是他們都需要訪問1000次的Parlist？爲什麼訪問Blist太慢？

內核_A，內核_D和內核_E具有相似的性能，這使我更加困惑：因此與kernel_B相比，額外的時間花在訪問Blist或sync上？

我想在我的PIC程序中實現kernel_B的性能。

Answer 1

不，不可能在塊之間進行同步。通常，內核調用帶來的開銷並不重要。我可以想象，你的內核不夠大，無法在很大程度上利用你的設備。如果你想檢查這個，你可以使用nvprof來分析你的程序並尋找瓶頸。

實現快速PIC代碼並不容易。你有沒有考慮過使用像PIConGPU這樣的庫？您可以在此鏈接下找到它：https://github.com/ComputationalRadiationPhysics/picongpu