title: CUDA median calculation through reduction
body: <p>I'm probably doing something incredibly stupid, but I can't seem to make this reduction work (there is probably a library that does this already, but this is for self-learning, so please bear with me). I'm trying to find the median of an array of integer entries by taking the median of medians approach, which I've coded below:</p>
  
  <pre><code>__global__ void gpuMedOdd(int *entries, int *med) {
          extern __shared__ int sdata[];
  
          int tid = threadIdx.x;
          int i = blockIdx.x * blockDim.x + threadIdx.x;
  
          sdata[tid] = entries[i];
          __syncthreads();
  
          for(int s = blockDim.x / 3; s &gt; 0; s /= 3) {
                  if(tid &lt; s) {
                          int list[3];
                          list[0] = sdata[tid], list[1] = sdata[tid + s], list[2] = sdata[tid + 2 * s];
                          if(list[1] &lt; list[0])
                                  swapGpu(list[1], list[0]);
                          if(list[2] &lt; list[0])
                                  swapGpu(list[2], list[0]);
                          if(list[2] &lt; list[1])
                                  swapGpu(list[2], list[1]);
  
                          sdata[tid] = list[1];
                  }
  
                  __syncthreads();
          }
  
          *med = sdata[0];
  }
  </code></pre>
  
  <p>I invoke this kernel function as:</p>
  
  <pre><code>gpuMedOdd&lt;&lt;&lt;9, numEntries / 9&gt;&gt;&gt;(d_entries, d_med);
  </code></pre>
  
  <p>I then copy the value in d_med over into med and print out that value. Unfortunately, this value is always 0, regardless of input. What am I doing wrong?</p>
  
  <p>Edit: I forgot to mention, swapGpu(a, b) is defined as below:</p>
  
  <pre><code>__device__ inline void swapGpu(int a, int b) {
          int dum = a;
          a = b;
          b = dum;
  }
  </code></pre>
  
  <p>Edit2: As suggested below, here is the entirety of the code.</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;fstream&gt;
  #include &lt;cstdlib&gt;
  
  #define checkCudaErrors(err) __checkCudaErrors(err, __FILE__, __LINE__)
  #define getLastCudaError(msg) __getLastCudaError(msg, __FILE__, __LINE__)
  
  inline void __checkCudaErrors(cudaError err, const char *file, const int line) {
          if(cudaSuccess != err) {
                  std::cout &lt;&lt; file &lt;&lt; "(" &lt;&lt; line &lt;&lt; ") : CUDA Runtime API error " &lt;&lt; (int) err &lt;&lt; ": " &lt;&lt; cudaGetErrorString(err) &lt;&lt; std::endl;
                  exit(3);
          }
  }
  
  inline void __getLastCudaError(const char *errorMsg, const char *file, const int line) {
          cudaError_t err = cudaGetLastError();
          if(cudaSuccess != err) {
                  std::cout &lt;&lt; file &lt;&lt; "(" &lt;&lt; line &lt;&lt; ") : getLastCudaError() CUDA error : " &lt;&lt; errorMsg &lt;&lt; " : (" &lt;&lt; (int) err &lt;&lt; ") " &lt;&lt; cudaGetErrorString(err) &lt;&lt; std::endl;
                  exit(3);
          }
  }
  
  int cpuMin(int *entries, int numEntries) {
          int minVal = entries[0];
  
          for(int i = 1; i &lt; numEntries; i++)
                  if(entries[i] &lt; minVal)
                          minVal = entries[i];
  
          return minVal;
  }
  
  int cpuMax(int *entries, int numEntries) {
          int maxVal = entries[0];
  
          for(int i = 1; i &lt; numEntries; i++)
                  if(entries[i] &gt; maxVal)
                          maxVal = entries[i];
  
          return maxVal;
  }
  
  inline void swap(int a, int b) {
          int dum = a;
          a = b;
          b = dum;
  }
  
  __device__ inline void swapGpu(int a, int b) {
          int dum = a;
          a = b;
          b = dum;
  }
  
  __global__ void gpuMedOdd(int *entries, int *med, int numEntries) {
          extern __shared__ int sdata[];
  
          int tid = threadIdx.x;
          int i = blockIdx.x * (blockDim.x * 3) + threadIdx.x;
  
          if(i + 2 * blockDim.x &lt; numEntries) {
                  int list[3];
                  list[0] = entries[i], list[1] = entries[i + blockDim.x], list[2] = entries[i + 2 * blockDim.x];
                  if(list[1] &lt; list[0])
                          swapGpu(list[1], list[0]);
                  if(list[2] &lt; list[0])
                          swapGpu(list[2], list[0]);
                  if(list[2] &lt; list[1])
                          swapGpu(list[2], list[1]);
  
                  sdata[tid] = list[1];
          }
  
          __syncthreads();
  
          for(int s = blockDim.x / 3; s &gt; 0; s /= 3) {
                  if(tid &lt; s &amp;&amp; tid + 2 * s &lt; blockDim.x) {
                          int list[3];
                          list[0] = sdata[tid], list[1] = sdata[tid + s], list[2] = sdata[tid + 2 * s];
                          if(list[1] &lt; list[0])
                                  swapGpu(list[1], list[0]);
                          if(list[2] &lt; list[0])
                                  swapGpu(list[2], list[0]);
                          if(list[2] &lt; list[1])
                                  swapGpu(list[2], list[1]);
  
                          sdata[tid] = list[1];
                  }
  
                  __syncthreads();
          }
  
          *med = sdata[0];
  }
  
  __global__ void gpuMin(int *entries, int *min, int numEntries) {
          extern __shared__ int sdata[];
  
          int tid = threadIdx.x;
          int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
  
          if(i + blockDim.x &lt; numEntries)
                  sdata[tid] = (entries[i] &lt; entries[i + blockDim.x]) ? entries[i] : entries[i + blockDim.x];
          __syncthreads();
  
          for(int s = blockDim.x / 2; s &gt; 0; s &gt;&gt;= 1) {
                  if(tid &lt; s)
                          sdata[tid] = (sdata[tid] &lt; sdata[tid + s]) ? sdata[tid] : sdata[tid + s];
  
                  __syncthreads();
          }
  
          *min = sdata[0];
  }
  
  __global__ void gpuMax(int *entries, int *max, int numEntries) {
          extern __shared__ int sdata[];
  
          int tid = threadIdx.x;
          int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
  
          if(i + blockDim.x &lt; numEntries)
                  sdata[tid] = (entries[i] &gt; entries[i + blockDim.x]) ? entries[i] : entries[i + blockDim.x];
          __syncthreads();
  
          for(int s = blockDim.x / 2; s &gt; 0; s &gt;&gt;= 1) {
                  if(tid &lt; s)
                          sdata[tid] = (sdata[tid] &gt; sdata[tid + s]) ? sdata[tid] : sdata[tid + s];
  
              __syncthreads();
      }
  
      *max = sdata[0];
  }
  
  int partition(int *entries, int left, int right, int pivotIdx) {
          int i, storeIdx = left, pivot = entries[pivotIdx];
  
          swap(entries[pivotIdx], entries[right]);
  
          for(i = left; i &lt; right; i++)
                  if(entries[i] &lt; pivot) {
                          swap(entries[i], entries[storeIdx]);
                          storeIdx++;
                  }
  
          return storeIdx;
  }
  
  int cpuSelect(int *entries, int left, int right, int k) {
          if(left == right)
                  return entries[left];
  
          int pivotIdx = ((left + right) &gt;&gt; 2) + 1, pivotNewIdx, pivotDist;
  
          pivotNewIdx = partition(entries, left, right, pivotIdx);
  
          pivotDist = pivotNewIdx - left + 1;
  
          if(pivotDist == k)
                  return entries[pivotNewIdx];
  
          else if(k &lt; pivotDist)
                  return cpuSelect(entries, left, pivotNewIdx - 1, k);
  
          else
                  return cpuSelect(entries, pivotNewIdx + 1, right, k - pivotDist);
  }
  
  int main(int argc, char *argv[]) {
          if(argc != 3) {
                  std::cout &lt;&lt; "ERROR: Incorrect number of input arguments" &lt;&lt; std::endl;
                  std::cout &lt;&lt; "Proper usage: " &lt;&lt; argv[0] &lt;&lt; " fileName numEntries" &lt;&lt; std::endl;
                  exit(1);
          }
  
          std::ifstream inp(argv[1]);
  
          if(!inp.is_open()) {
                  std::cout &lt;&lt; "ERROR: File I/O error" &lt;&lt; std::endl;
                  std::cout &lt;&lt; "Could not find file " &lt;&lt; argv[1] &lt;&lt; std::endl;
                  exit(2);
          }
  
          int numEntries = atoi(argv[2]), i = 0;
  
          int *entries = new int[numEntries];
  
          while(inp &gt;&gt; entries[i] &amp;&amp; i &lt; numEntries)
                  i++;
  
          if(i &lt; numEntries) {
                  std::cout &lt;&lt; "ERROR: File I/O error" &lt;&lt; std::endl;
                  std::cout &lt;&lt; "Command-line input suggested " &lt;&lt; numEntries &lt;&lt; " entries, but only found " &lt;&lt; i &lt;&lt; " entries" &lt;&lt; std::endl;
                  exit(2);
          }
  
          if(inp &gt;&gt; i) {
                  std::cout &lt;&lt; "ERROR: File I/O error" &lt;&lt; std::endl;
                  std::cout &lt;&lt; "Command-line input suggested " &lt;&lt; numEntries &lt;&lt; " entries, but file contains more entries" &lt;&lt; std::endl;
                  exit(2);
          }
  
          int min, max;
          int *d_entries, *d_min, *d_max;
  
          checkCudaErrors(cudaMalloc(&amp;d_entries, sizeof(int) * numEntries));
          checkCudaErrors(cudaMalloc(&amp;d_min, sizeof(int)));
          checkCudaErrors(cudaMalloc(&amp;d_max, sizeof(int)));
  
          checkCudaErrors(cudaMemcpy(d_entries, entries, sizeof(int) * numEntries, cudaMemcpyHostToDevice));
  
          gpuMin&lt;&lt;&lt;16, numEntries / 16, numEntries / 16 * sizeof(int)&gt;&gt;&gt;(d_entries, d_min, numEntries);
          getLastCudaError("kernel launch failure");
          gpuMax&lt;&lt;&lt;16, numEntries / 16, numEntries / 16 * sizeof(int)&gt;&gt;&gt;(d_entries, d_max, numEntries);
          getLastCudaError("kernel launch failure");
  
          checkCudaErrors(cudaMemcpy(&amp;min, d_min, sizeof(int), cudaMemcpyDeviceToHost));
          checkCudaErrors(cudaMemcpy(&amp;max, d_max, sizeof(int), cudaMemcpyDeviceToHost));
  
          std::cout &lt;&lt; "The minimum value is: " &lt;&lt; min &lt;&lt; std::endl;
          std::cout &lt;&lt; "The maximum value is: " &lt;&lt; max &lt;&lt; std::endl;
  
          if(numEntries % 2) {
                  int med, *d_med;
                  checkCudaErrors(cudaMalloc(&amp;d_med, sizeof(int)));
  
                  gpuMedOdd&lt;&lt;&lt;16, numEntries / 16, 16 * sizeof(int)&gt;&gt;&gt;(d_entries, d_med, numEntries);
                  getLastCudaError("kernel launch failure");
  
                  checkCudaErrors(cudaMemcpy(&amp;med, d_med, sizeof(int), cudaMemcpyDeviceToHost));
  
                  std::cout &lt;&lt; "The median value is: " &lt;&lt; med &lt;&lt; std::endl;
          }
          else {
                  int *d_med;
                  cudaMalloc(&amp;d_med, sizeof(int));
                  gpuMedOdd&lt;&lt;&lt;16, numEntries / 16&gt;&gt;&gt;(d_entries, d_med, numEntries);
          }
  
          min = cpuMin(entries, numEntries);
  
          max = cpuMax(entries, numEntries);
  
          if(numEntries % 2) {
                  int median = cpuSelect(entries, 0, numEntries - 1, (numEntries - 1) / 2 + 1);
                  std::cout &lt;&lt; "The median value is: " &lt;&lt; median &lt;&lt; std::endl;
          }
  
          else {
                  int med2 = cpuSelect(entries, 0, numEntries - 1, numEntries / 2);
                  int med1 = cpuSelect(entries, 0, numEntries - 1, numEntries / 2 + 1);
  
                  float median = 0.5 * (med1 + med2);
                  std::cout &lt;&lt; "The median value is: " &lt;&lt; median &lt;&lt; std::endl;
          }
  
          std::cout &lt;&lt; "The minimum value is: " &lt;&lt; min &lt;&lt; std::endl;
          std::cout &lt;&lt; "The maximum value is: " &lt;&lt; max &lt;&lt; std::endl;
  
          exit(0);
  }
  </code></pre>