title: Nested loops modulo permutation in cuda
body: <p>I need to perform a function on triplets taken from an array and add the result to a Histogram, but I want to avoid permutations since the function is invariant under those [F(i,j,k) = F(j,i,k) and so on].</p>
  
  <p>Normally I would code something like this:</p>
  
  <pre><code>def F(int i, int j, int k){
      int temp_index;
      /* Do something */
      return temp_index;
  }
  
  for(int i=0;i&lt;N;i++){
      for(int j=i+1;j&lt;N;j++){
          for(int k=j+1;k&lt;N;k++){
              hist[F(i,j,k)]++;
          }
      }
  }
  </code></pre>
  
  <p>As N is quite big (approx. 10^5), I would like to call perform this on a GPU using cuda.</p>
  
  <p>I have written a code to call this function on the GPU, but I have no idea how to prevent multiple calls of the same triple of indices. So far I call cuda with a 3-dimensional grid, like:</p>
  
  <pre><code>__global__ void compute_3pcf(float *theta, float *hist) {
    int i,j,k;
    i = blockIdx.x*blockDim.x + threadIdx.x;
    j = blockIdx.y*blockDim.y + threadIdx.y;
    k = blockIdx.z*blockDim.z + threadIdx.z;
    if(i&gt;=j || j&gt;=k) return;
    atomicAdd(&amp;hist[F(i,j,k)],1);
  }
  
  int main(){
    /*
    Allocation of memory and cudaMemcpy
    */
    dim3 grid((N+15)/16,(N+7)/8,(N+7)/8);
    dim3 block(16,8,8);
  
    //Launch on GPU
    compute_3pcf&lt;&lt;&lt;grid,block&gt;&gt;&gt;(d_theta, d_hist);
  }
  </code></pre>
  
  <p>However, now for each combination (i,j,k) a new thread is launched and then aborted, which seems very inefficient to me, as then only 1/6 of the threads perform the actual computation. What I would like to have is something like this:</p>
  
  <pre><code>__global__ void compute_3pcf(float *theta, float *hist) {
    int i,j,k,idx;
    idx = blockIdx.x*blockDim.x + threadIdx.x;
    i = H_i(idx);
    j = H_j(idx,i);
    k = H_k(idx,j);
    atomicAdd(&amp;hist[F(i,j,k)],1);
  }
  
  int main(){
    /*
    Allocation of memory and cudaMemcpy
    */
    long long int N_combinations = N*(N-1)*(N-2)/6;
    long int grid = (N_combinations+1023)/1024;
    int block = 1024;
    //Launch on GPU
    compute_3pcf&lt;&lt;&lt;grid,block&gt;&gt;&gt;(d_theta, d_hist);
  }
  </code></pre>
  
  <p>However, I am unable to find the functions H_i, H_j, H_k. If anyone can tell me how I could solve or avoid this problem, I would be very thankful.</p>
  
  <p>Edit: The histogram contains about 10^6 bins, so that I can not have one histogram per block in a shared memory, like in the example code for cuda. Instead, it lies in the global memory of the GPU.</p>