title: Techniques for increasing arithmetic intensity of elementwise CUDA kernel
body: <p>I have written a CUDA kernel for calculating pairwise interactions between a set of source particles and target particles.</p>
  <p>My M target particles look like,</p>
  <pre><code>[[x_1,y_1,z_1],...,[x_m, y_m, z_m]]
  </code></pre>
  <p>And my N source particles look like</p>
  <pre><code>[[x_1,y_1,z_1],...,[x_n, y_n, z_n]]
  </code></pre>
  <p>with M &lt; N.</p>
  <p>I begin by transferring all source/target data to the GPU, and looping through batches of the source particles and evaluating the pairwise interactions with <em>all</em> targets.</p>
  <p>Like so (in cupy syntax)</p>
  <pre class="lang-py prettyprint-override"><code>for i in range(n_blocks):
      left_idx = i*width
      right_idx = (i+1)*width
      gpu_func(
         grid_dimensions, block_dimensions, targets, 
         sources[left_idx:right_idx,:], width
      )
  </code></pre>
  <p>where <code>sources</code> and <code>targets</code> are arrays on the GPU containing the source and target data.</p>
  <p>My question is what do I have to do to avoid this loop? I am a novice to CUDA. My thinking is to check whether the global thread indices from each thread satisfy the 'left_idx' and 'right_index' conditions I've got above in my host device loop, is this correct? Is there a better way of doing this? I feel like I'm massively underutilising the GPU currently as the number of sources/targets is significantly less than the number of CUDA cores on my machine for my problem.</p>