title: How to write an efficient CUDA program with two stages
body: <p>I want to write a lightweight PIC (Particle-in-cell) program. By "lightweight" I mean it doesn't need to scale up: just assume all data can fit into both the memory of a single GPU device and the memory of the host system. However I want it to be as fast as possible.</p>
  
  <p>The problem is, the typical structure of PIC is the interation of two stages: field solver and particle pusher. The workflow is like such:
  Initialize system -> push particle -> solve field -> push particle -> solve field... -> output</p>
  
  <p>The next push particle or solve field must wait until the previous solve field or push particle to finish. It could take millions of iteration to get the final output.</p>
  
  <p>As a test, omitting the field solver, the particle pusher can be written as:</p>
  
  <pre><code> __device__
   void push(Particle &amp;par) {
     // some routines to move a particle. same excecutiong time for every particle.
   }
  </code></pre>
  
  <p>And use kernel_1 like this to excecute it:</p>
  
  <pre><code>__global__
  void kernel_1(int n, Particle* parlist)
  {
    int i = blockIdx.x*blockDim.x + threadIdx.x;
    if (i &lt; n) {
        push(parlist[i]);
    }
  }
  </code></pre>
  
  <p>In the main loop</p>
  
  <pre><code>for (int i=0;i&lt;M;i++) {
      kernel_1&lt;&lt;&lt;(n+255)/256, 256&gt;&gt;&gt;(n, parlist);
  }
  </code></pre>
  
  <p>M is the required number of iteration. However, the performance is prohibitively slow: On my system with octa-core Intel E5-2640 v3 and Nvidia Quadro m4000, CUDA gives similiar performance as a pure CPU version using openmp. For a particle number of 10,000,000 and M=1000, it takes around 10s.</p>
  
  <p>However if I move the loop into the kernel:</p>
  
  <pre><code>void kernel_2(int n, Particle* parlist)
  {
    int i = blockIdx.x*blockDim.x + threadIdx.x;
    if (i &lt; n) {
        for (int i=0;i&lt;M;i++) {
            push(parlist[i]);
        }
    }
  }
  </code></pre>
  
  <p>and</p>
  
  <pre><code>kernel_2&lt;&lt;&lt;(n+255)/256, 256&gt;&gt;&gt;(n, parlist);
  </code></pre>
  
  <p>For the same M=1000, it takes only 100ms, which is a 10000% speedup. I have verified the results are the same and correct for both cases. Maybe the calling cost of running a kernel for M times is too high.</p>
  
  <p>The performance improvement of moving the loop into the kernel is so unbelieveable yet true. For the first case, it's easy to add the field solver: just write a new kernel and excecute two kernels sequentially in the main loop. However the performance should be medicore.</p>
  
  <p>I find it difficult to add the field solver routine into the second case: there doesn't seem to be a sychronization mechanism between blocks without calling kernels for multiple times, however the field solver must wait until ALL particles are pushed, which has to be distributed to different blocks (because the number of particles is very high).</p>
  
  <p>So is it possible to implement a two stage iteration inside one kernel? The performance gain is too much to be neglected.</p>
  
  <p><strong>EDIT:</strong>
  I found out the performance difference very confusing: the difference of 100ms and 10s is just a single line of code or even the sequence of loops. I have modified the push() to be a little more complex (a 2d Boris pusher):</p>
  
  <pre><code>class Particle
  {
  public:
      float x, y;        //m
      float vx, vy;      //m/s
      float m;           //kg
      float q;           //ee
  };
  
  __device__
  void run(Particle&amp; 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;
  }
  </code></pre>
  
  <p>I created 1 n-element array for Particle and 1 n-element float array for B. They are created on host and cudaMemcpy to device. I then examined the performance of the following three kernels:</p>
  
  <pre><code>__global__
  void kernel_A(int n, int m, Particle* parlist, float* Blist)
  {
      int i = blockIdx.x*blockDim.x + threadIdx.x;
      int j;
      if (i&lt;n) {
          for (j=0;j&lt;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&lt;n) {
          B = Blist[i];
          for (j=0;j&lt;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&lt;n) {
          B = Blist[i];
          for (j=0;j&lt;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&lt;n) {
          B = Blist[i];
      }
      for (j=0;j&lt;m;j++) {
          if (i&lt;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&lt;n) {
          for (j=0;j&lt;m;j++) {
              run(parlist[i], Blist[i]);
              __syncthreads();
          }
      }
  }
  </code></pre>
  
  <p>And the run time is quite different. For n=10,000,000 and m=1000:</p>
  
  <ul>
  <li>Kernel_A : 7.6s  </li>
  <li>Kernel_B : 66ms</li>
  <li>Kernel_C : 9.9s</li>
  <li>Kernel_D : 10.0s</li>
  <li>Kernel_E : 10.0s</li>
  </ul>
  
  <p>The results of the three kernels are all the same and correct (checked against the CPU version).</p>
  
  <p>I understand from official CUDA Programming guide that branching is expensive so kernel_C should be slower than kernel_B, though I suspect the difference is of two order of magnitude. What I don't understand is why kernel_B performs so much better than kernel_A. Kernel_B doesn't have to access Blist for 1000 times while kernel_A does, but however they both need to access parlist for 1000 times right? And why is accessing Blist so slow? </p>
  
  <p>Kernel_A, kernel_D and kernel_E have similiar performance, which confuse me further: so the extra time compared to kernel_B is spent on accessing Blist or sync?</p>
  
  <p>I want to achieve the performance of kernel_B in my PIC program.</p>