title: CUDA Reduction: Warp Unrolling (School)
body: <p>I am currently working on a project in which I am unrolling the last warp of a reduction. I have finished the code above; however, some modifications were done by guessing and I'd like an explanation why. The code I have written is only the function kernel4</p>
  
  <pre><code>// in is input array, out is where to store result, n is number of elements from in
  // T is a float (32bit)
  __global__ void kernel4(T *in, T *out, unsigned int n)
  </code></pre>
  
  <p>which is a reduction algorithm, the rest of the code was already provided. </p>
  
  <p>Code:</p>
  
  <pre><code>#include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  
  #include "timer.h"
  #include "cuda_utils.h"
  
  typedef float T;
  
  #define N_ (8 * 1024 * 1024)
  #define MAX_THREADS 256
  #define MAX_BLOCKS 64
  
  #define MIN(x,y) ((x &lt; y) ? x : y)
  #define tid threadIdx.x 
  #define bid blockIdx.x 
  #define bdim blockDim.x
  #define warp_size 32
  
  unsigned int nextPow2( unsigned int x ) {
      --x;
      x |= x &gt;&gt; 1;
      x |= x &gt;&gt; 2;
      x |= x &gt;&gt; 4;
      x |= x &gt;&gt; 8;
      x |= x &gt;&gt; 16;
      return ++x;
  }
  
  void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &amp;blocks, int &amp;threads)
  {
      if (whichKernel &lt; 3) {
          threads = (n &lt; maxThreads) ? nextPow2(n) : maxThreads;
          blocks = (n + threads - 1) / threads;
      } else {
          threads = (n &lt; maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
          blocks = (n + (threads * 2 - 1)) / (threads * 2);
      }
      if (whichKernel == 5)
          blocks = MIN(maxBlocks, blocks);
  }
  
  T reduce_cpu(T *data, int n) {
      T sum = data[0];
      T c = (T) 0.0;
      for (int i = 1; i &lt; n; i++)
      {
          T y = data[i] - c;
          T t = sum + y;
          c = (t - sum) - y;
          sum = t;
      } 
      return sum;
  }
  
  __global__ void
  kernel4(T *in, T *out, unsigned int n)
  {
      __shared__ volatile T d[MAX_THREADS];
  
      unsigned int i = bid * bdim + tid;
  
      n &gt;&gt;= 1;
      d[tid] = (i &lt; n) ? in[i] + in[i+n] : 0;
      __syncthreads ();
  
      for(unsigned int s = bdim &gt;&gt; 1; s &gt; warp_size; s &gt;&gt;= 1) {
          if(tid &lt; s)
              d[tid] += d[tid + s];
          __syncthreads ();
      }
  
      if (tid &lt; warp_size) {
          if (n &gt; 64) d[tid] += d[tid + 32];
          if (n &gt; 32) d[tid] += d[tid + 16];
          d[tid] += d[tid + 8];
          d[tid] += d[tid + 4];
          d[tid] += d[tid + 2];
          d[tid] += d[tid + 1];
      }
  
      if(tid == 0)
          out[bid] = d[0];
  }
  
  int main(int argc, char** argv)
  {
      T *h_idata, h_odata, h_cpu;
      T *d_idata, *d_odata;   
  
      struct stopwatch_t* timer = NULL;
      long double t_kernel_4, t_cpu;
  
      int whichKernel = 4, threads, blocks, N, i;
      if(argc &gt; 1) {
          N = atoi (argv[1]);
          printf("N: %d\n", N);
      } else {
          N = N_;
          printf("N: %d\n", N);
      }
  
      getNumBlocksAndThreads (whichKernel, N, MAX_BLOCKS, MAX_THREADS, blocks, threads);
  
      stopwatch_init ();
      timer = stopwatch_create ();
  
      h_idata = (T*) malloc (N * sizeof (T));
      CUDA_CHECK_ERROR (cudaMalloc (&amp;d_idata, N * sizeof (T)));
      CUDA_CHECK_ERROR (cudaMalloc (&amp;d_odata, blocks * sizeof (T)));
  
      srand48(time(NULL));
      for(i = 0; i &lt; N; i++)
          h_idata[i] = drand48() / 100000;
      CUDA_CHECK_ERROR (cudaMemcpy (d_idata, h_idata, N * sizeof (T), cudaMemcpyHostToDevice));
  
      dim3 gb(blocks, 1, 1);
      dim3 tb(threads, 1, 1);
  
      kernel4 &lt;&lt;&lt;gb, tb&gt;&gt;&gt; (d_idata, d_odata, N);
      cudaThreadSynchronize ();
  
      stopwatch_start (timer);
  
      kernel4 &lt;&lt;&lt;gb, tb&gt;&gt;&gt; (d_idata, d_odata, N);
      int s = blocks;
      while(s &gt; 1) {
          threads = 0;
          blocks = 0;
          getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
  
          dim3 gb(blocks, 1, 1);
          dim3 tb(threads, 1, 1);
  
          kernel4 &lt;&lt;&lt;gb, tb&gt;&gt;&gt; (d_odata, d_odata, s);
          s = (s + threads * 2 - 1) / (threads * 2);
      }
      cudaThreadSynchronize ();
  
      t_kernel_4 = stopwatch_stop (timer);
      fprintf (stdout, "Time to execute unrolled GPU reduction kernel: %Lg secs\n", t_kernel_4);
  
      double bw = (N * sizeof(T)) / (t_kernel_4 * 1e9);   // total bits / time
      fprintf (stdout, "Effective bandwidth: %.2lf GB/s\n", bw);
      CUDA_CHECK_ERROR (cudaMemcpy (&amp;h_odata, d_odata, sizeof (T), cudaMemcpyDeviceToHost));
  
      stopwatch_start (timer);
      h_cpu = reduce_cpu (h_idata, N);
      t_cpu = stopwatch_stop (timer);
      fprintf (stdout, "Time to execute naive CPU reduction: %Lg secs\n", t_cpu);
  
      if(abs (h_odata - h_cpu) &gt; 1e-5)
          fprintf(stderr, "FAILURE: GPU: %f  CPU: %f\n", h_odata, h_cpu);
      else
          printf("SUCCESS: GPU: %f  CPU: %f\n", h_odata, h_cpu);
      return 0;
  }
  </code></pre>
  
  <p>My first question is: when declaring</p>
  
  <pre><code>__shared__ volatile T d[MAX_THREADS];
  </code></pre>
  
  <p>I would like to verify my understanding of volatile. Volatile prevents compilers from incorrectly optimizing my code and promises that load/stores are completed through the cache and not just registers (please correct me if wrong). For reduction, if partial reduction sums are still stored in registers, why is this a problem?</p>
  
  <p>My second question is: when doing the actual warp reduction</p>
  
  <pre><code>    if (tid &lt; warp_size) { // Final log2(32) = 5 strides
          if (n &gt; 64) d[tid] += d[tid + 32];
          if (n &gt; 32) d[tid] += d[tid + 16];
          d[tid] += d[tid + 8];
          d[tid] += d[tid + 4];
          d[tid] += d[tid + 2];
          d[tid] += d[tid + 1];
      }
  </code></pre>
  
  <p>The reduction sum will yield incorrect results without (n > 64) and (n > 32) conditions. The results I get are:</p>
  
  <pre><code>FAILURE: GPU: 41.966557  CPU: 41.946209
  </code></pre>
  
  <p>With 5 trials, the GPU reduction consistently yields an error of 0.0204. I am wary to think this is a floating point operation error.</p>
  
  <p>To be honest as well, my teacher's assistant suggested this change to add the (n > 64) and (n > 32) conditions but did not explain why it would fix the code.</p>
  
  <p>Since n in my trials are over 64, why does this conditional change the results. I am having difficulty tracing back the problem because I cannot use print functions like I would in a CPU.</p>