title: radix select using cuda
body: <p>I have been working to develop a radix select using CUDA which utilizes k smallest element to sort given number of elements. The main idea behind this radix select is that is scans through 32 bit integer starting from its MSB to LSB. It partitions all 0 bit on left side and all 1 bit on the right side. The side with contains k smallest elements is solved recursively. My partition process works just fine but I am having problem dealing with recursive function calls. I am unable to stop the recursion. Please help me on that!
  My kernel function looks like this: This is kernel.h</p>
  
  <pre><code>  #include "header.h"
    #define WARP_SIZE 32
    #define BLOCK_SIZE 32
  
  __device__ int Partition(int *d_DataIn, int firstidx, int lastidx, int k, int N, int bit)
  {
  int threadID = threadIdx.x + BLOCK_SIZE * blockIdx.x;
  int WarpID = threadID &gt;&gt; 5;
  int LocWarpID = threadID - 32 * WarpID;
  int NumWarps = N / WARP_SIZE;
  int pivot;
  
  __shared__ int DataPartition[BLOCK_SIZE];
  __shared__ int DataBinary[WARP_SIZE];
  
  for(int i = 0; i &lt; NumWarps; i++)
  {
      if(LocWarpID &gt;= firstidx &amp;&amp; LocWarpID &lt;=lastidx)
      {
          int r = d_DataIn[i * WARP_SIZE + LocWarpID];
          int p = (r&gt;&gt;(31-bit))&amp;1;
          unsigned int B = __ballot(p);
          unsigned int B_flip = ~B;
          if(p==1)
          {
              int b = B &lt;&lt; (32-LocWarpID);
              int RightLoc = __popc(b);
              DataPartition[lastidx - RightLoc] = r;
          }
          else
          {
              int b_flip = B_flip &lt;&lt; (32 - LocWarpID);
              int LeftLoc = __popc(b_flip);
              DataPartition[LeftLoc] = r;
          }
  
      if(LocWarpID &lt;= lastidx - __popc(B))
          {
  
              d_DataIn[LocWarpID] = DataPartition[LocWarpID];
          }
      else
          {
  
              d_DataIn[LocWarpID] = DataPartition[LocWarpID];
          }
          pivot = lastidx - __popc(B); 
          return pivot+1;
      }
  
  }   
  }
  
  
  __device__ int RadixSelect(int *d_DataIn, int firstidx, int lastidx, int k, int N, int bit)
  {
  
  if(firstidx == lastidx)
      return *d_DataIn;
  int q = Partition(d_DataIn, firstidx, lastidx, k, N, bit);
  int length = q - firstidx;
  if(k == length)
      return *d_DataIn; 
  else if(k &lt; length)
      return RadixSelect(d_DataIn, firstidx, q-1, k, N, bit+1);
  else
      return RadixSelect(d_DataIn, q, lastidx, k-length, N, bit+1);
  }
  
   __global__ void radix(int *d_DataIn, int firstidx, int lastidx, int k, int N, int bit)
  {
  RadixSelect(d_DataIn, firstidx, lastidx, k, N, bit);
  }
  </code></pre>
  
  <p>Host code is main.cu and it looks like:</p>
  
  <pre><code>#include "header.h"
  #include &lt;iostream&gt;
  #include &lt;fstream&gt;
  #include "kernel.h"
  
  #define BLOCK_SIZE 32
  using namespace std;
  
  int main()
  {
  int N = 32;
  thrust::host_vector&lt;float&gt;h_HostFloat(N);
  thrust::counting_iterator &lt;unsigned int&gt; Numbers(0);
  thrust::transform(Numbers, Numbers + N, h_HostFloat.begin(),         RandomFloatNumbers(1.f, 100.f));
  thrust::host_vector&lt;int&gt;h_HostInt(N);
  thrust::transform(h_HostFloat.begin(), h_HostFloat.end(), h_HostInt.begin(), FloatToInt());
  thrust::device_vector&lt;float&gt;d_DeviceFloat = h_HostFloat;
  thrust::device_vector&lt;int&gt;d_DeviceInt(N);   
  thrust::transform(d_DeviceFloat.begin(), d_DeviceFloat.end(), d_DeviceInt.begin(), FloatToInt());
  int *d_DataIn = thrust::raw_pointer_cast(d_DeviceInt.data());
  int *h_DataOut;
  float *h_DataOut1;
  int fsize = N * sizeof(float);
  int size = N * sizeof(int);
  h_DataOut = new int[size];
  h_DataOut1 = new float[fsize];
  
  int firstidx = 0;
  int lastidx = BLOCK_SIZE-1;
  int k = 20;
  int bit = 1;
  int NUM_BLOCKS = N / BLOCK_SIZE;
  radix &lt;&lt;&lt; NUM_BLOCKS, BLOCK_SIZE &gt;&gt;&gt; (d_DataIn, firstidx, lastidx, k, N, bit);
  cudaMemcpy(h_DataOut, d_DataIn, size, cudaMemcpyDeviceToHost);  
  WriteData(h_DataOut1, h_DataOut, 10, N);
  
  
  return 0;
     }
  </code></pre>
  
  <p>List of headers that I used: </p>
  
  <pre><code>   #include "cuda.h"
     #include "cuda_runtime_api.h"
     #include "device_launch_parameters.h"
     #include &lt;thrust/host_vector.h&gt;
     #include &lt;thrust/device_vector.h&gt;
     #include &lt;thrust/transform.h&gt;
     #include &lt;thrust/generate.h&gt;
     #include "functor.h"
     #include &lt;thrust/iterator/counting_iterator.h&gt;
     #include &lt;thrust/copy.h&gt;
     #include &lt;thrust/device_ptr.h&gt;
  </code></pre>
  
  <p>Another header file "functor.h" to convert floating point numbers to int type and to generate random floating numbers. </p>
  
  <pre><code>  #include &lt;thrust/random.h&gt;
    #include &lt;sstream&gt;
    #include &lt;fstream&gt;
    #include &lt;iomanip&gt;
    struct RandomFloatNumbers
     {
  float a, b;
  __host__ __device__
      RandomFloatNumbers(float _a, float _b) : a(_a), b(_b) {};  
  __host__ __device__
      float operator() (const unsigned int n) const{
          thrust::default_random_engine rng;
          thrust::uniform_real_distribution&lt;float&gt; dist(a,b);
          rng.discard(n);
          return dist(rng);
  }
  
      };
  
     struct FloatToInt
      {
  __host__ __device__
      int operator() (const float &amp;x)
  const {
      union {
          float f_value;
          int i_value;
      } value;
  
      value.f_value = x;
      return value.i_value;
  }
      };
  
     float IntToFloat(int &amp;x)
       {
      union{
            float f_value;
            int i_value;
          }value;
  
          value.i_value = x;
          return value.f_value;
       }
  
     bool WriteData(float *h_DataOut1, int *h_DataOut, int bit, int N)
      {
  std::ofstream data;
  std::stringstream file;
  file &lt;&lt; "out\\Partition_";
  file &lt;&lt; std::setfill('0') &lt;&lt;std::setw(2) &lt;&lt; bit;
  file &lt;&lt; ".txt";
  data.open((file.str()).c_str());
  if(data.is_open() ==  false)
  {
      std::cout &lt;&lt; "File is not open" &lt;&lt; std::endl;
      return false;
  }
  
      for(int i = 0; i &lt; N; i++)
      {
  
          h_DataOut1[i] = IntToFloat(h_DataOut[i]);
          //cout &lt;&lt; h_HostFloat[i] &lt;&lt; " \t" &lt;&lt;  h_DataOut1[i] &lt;&lt; endl;
          //std::bitset&lt;32&gt;bitshift(h_DataOut[i]&amp;1&lt;&lt;31-bit);
          //data &lt;&lt; bitshift[31-bit] &lt;&lt; "\t" &lt;&lt;h_DataOut1[i] &lt;&lt;std::endl;
          data &lt;&lt; h_DataOut1[i] &lt;&lt; std::endl;
      }
      data &lt;&lt; std::endl;
  data.close();
  std::cout &lt;&lt; "Partition=" &lt;&lt;bit &lt;&lt;"\n"; 
  return true;
      }
  </code></pre>