title: CUDA - why is warp based parallel reduction slower?
body: <p>I had the idea about a warp based parallel reduction since all threads of a warp are in sync by definition.   </p>
  
  <p>So the idea was that the input data can be reduced by factor 64 (each thread reduces two elements) without any synchronization need.</p>
  
  <p>Same as the original implementation by Mark Harris the reduction is applied on block-level and data is on shared memory.  <a href="http://gpgpu.org/static/sc2007/SC07_CUDA_5_Optimization_Harris.pdf">http://gpgpu.org/static/sc2007/SC07_CUDA_5_Optimization_Harris.pdf</a></p>
  
  <p>I created a kernel to test his version and my warp based version.<br>
  The kernel itself is completely identically storing BLOCK_SIZE elements in shared memory and outputting its result at its unique block index in an output array.</p>
  
  <p>The algorithm itself works fine. Tested with full array of one's to test the "counting".</p>
  
  <p><strong>Function body of the implementations:</strong></p>
  
  <pre><code>/**
   * Performs a parallel reduction with operator add 
   * on the given array and writes the result with the thread 0
   * to the given target value
   *
   * @param inValues T* Input float array, length must be a multiple of 2 and equal to blockDim.x
   * @param targetValue float 
   */
  __device__ void reductionAddBlockThread_f(float* inValues,
      float &amp;outTargetVar)
  {
      // code of the below functions
  }
  </code></pre>
  
  <p><strong>1. Implementation of his version:</strong></p>
  
  <pre><code>if (blockDim.x &gt;= 1024 &amp;&amp; threadIdx.x &lt; 512)
      inValues[threadIdx.x] += inValues[threadIdx.x + 512];
  __syncthreads();
  if (blockDim.x &gt;= 512 &amp;&amp; threadIdx.x &lt; 256)
      inValues[threadIdx.x] += inValues[threadIdx.x + 256];
  __syncthreads();
  if (blockDim.x &gt;= 256 &amp;&amp; threadIdx.x &lt; 128)
      inValues[threadIdx.x] += inValues[threadIdx.x + 128];
  __syncthreads();
  if (blockDim.x &gt;= 128 &amp;&amp; threadIdx.x &lt; 64)
      inValues[threadIdx.x] += inValues[threadIdx.x + 64];
  __syncthreads();
  
  //unroll last warp no sync needed
  if (threadIdx.x &lt; 32)
  {
      if (blockDim.x &gt;= 64) inValues[threadIdx.x] += inValues[threadIdx.x + 32];
      if (blockDim.x &gt;= 32) inValues[threadIdx.x] += inValues[threadIdx.x + 16];
      if (blockDim.x &gt;= 16) inValues[threadIdx.x] += inValues[threadIdx.x + 8];
      if (blockDim.x &gt;= 8) inValues[threadIdx.x] += inValues[threadIdx.x + 4];
      if (blockDim.x &gt;= 4) inValues[threadIdx.x] += inValues[threadIdx.x + 2];
      if (blockDim.x &gt;= 2) inValues[threadIdx.x] += inValues[threadIdx.x + 1];
  
      //set final value
      if (threadIdx.x == 0)
          outTargetVar = inValues[0];
  }
  </code></pre>
  
  <p><strong>Ressources:</strong></p>
  
  <p>4 syncthreads used<br>
  12 if statements used<br>
  11 read + add + write operations<br>
  1 final write operation<br>
  5 register usage  </p>
  
  <p><strong>Performance:</strong></p>
  
  <p>five test runs average: ~ 19.54 ms</p>
  
  <p><strong>2. Warp based approach:</strong> (Same function body as above)</p>
  
  <pre><code>/*
   * Perform first warp based reduction by factor of 64
   *
   * 32 Threads per Warp -&gt; LOG2(32) = 5
   *
   * 1024 Threads / 32 Threads per Warp = 32 warps
   * 2 elements compared per thread -&gt; 32 * 2 = 64 elements per warp
   *
   * 1024 Threads/elements divided by 64 = 16
   * 
   * Only half the warps/threads are active
   */
  if (threadIdx.x &lt; blockDim.x &gt;&gt; 1)
  {
      const unsigned int warpId = threadIdx.x &gt;&gt; 5;
      // alternative threadIdx.x &amp; 31
      const unsigned int threadWarpId = threadIdx.x - (warpId &lt;&lt; 5);
      const unsigned int threadWarpOffset = (warpId &lt;&lt; 6) + threadWarpId;
  
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 32];
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 16];
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 8];
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 4];
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 2];
      inValues[threadWarpOffset] += inValues[threadWarpOffset + 1];
  }
  
  // synchronize all warps - the local warp result is stored
  // at the index of the warp equals the first thread of the warp
  __syncthreads();
  
  // use first warp to reduce the 16 warp results to the final one
  if (threadIdx.x &lt; 8)
  {
      // get first element of a warp
      const unsigned int warpIdx = threadIdx.x &lt;&lt; 6;
  
      if (blockDim.x &gt;= 1024) inValues[warpIdx] += inValues[warpIdx + 512];
      if (blockDim.x &gt;= 512) inValues[warpIdx] += inValues[warpIdx + 256];
      if (blockDim.x &gt;= 256) inValues[warpIdx] += inValues[warpIdx + 128];
      if (blockDim.x &gt;= 128) inValues[warpIdx] += inValues[warpIdx + 64];
  
      //set final value
      if (threadIdx.x == 0)
          outTargetVar = inValues[0];
  }
  </code></pre>
  
  <p><strong>Ressources:</strong></p>
  
  <p>1 syncthread used<br>
  7 if statements<br>
  10 read add write operations<br>
  1 final write operation<br>
  5 register usage  </p>
  
  <p>5 bit shifts<br>
  1 add<br>
  1 sub  </p>
  
  <p><strong>Performance:</strong></p>
  
  <p>five test runs average: ~ 20.82 ms</p>
  
  <p>Testing both kernels multiple times on a <strong>Geforce 8800 GT 512 mb</strong> with <strong>256</strong> mb of float values. 
  And running kernel with <strong>256 threads per block</strong> (100 % occupancy).</p>
  
  <p>The warp based version is ~ <strong>1.28</strong> milliseconds slower.</p>
  
  <p><em>If future card's allow larger block sizes the warp based approach would still need no further sync statement since the max is 4096 which get reduced to 64 which get reduced by final warp to 1</em></p>
  
  <p><strong>Why is it not faster?, or where is the flaw in the idea, kernel?</strong> </p>
  
  <p>From ressources usage the warp approach should be ahead?</p>
  
  <p><strong>Edit1: Corrected the kernel that only half the threads are active not resulting in out of bound reads, added new performance data</strong></p>