title: Using CUDA atomicInc to get unique indices
body: <p>I have CUDA kernel where basically each thread holds a value, and it needs to add that value to one or more lists in shared memory. So for each of those lists, it needs to get an index value (unique for that list) to put the value.</p>
  <p>The real code is different, but there are lists like:</p>
  <pre><code>typedef struct {
      unsigned int numItems;
      float items[MAX_NUM_ITEMS];
  } List;
  __shared__ List lists[NUM_LISTS];
  </code></pre>
  <p>The values <code>numItems</code> are initially all set to 0, and then a <code>__syncthreads()</code> is done.</p>
  <p>To add its value to the lists, each thread does:</p>
  <pre><code>for(int list = 0; list &lt; NUM_LISTS; ++list) {
      if(should_add_to_list(threadIdx, list)) {
          unsigned int index = atomicInc(&amp;lists[list].numItems, 0xffffffff);
          assert(index &lt; MAX_NUM_ITEMS); // always true
          lists[list].items[index] = my_value;
      }
  }
  </code></pre>
  <p>This works most of the time, but it seems that when making some unrelated changes in other parts of the kernel (such as not checking asserts that always succeed), sometimes two threads get the same index for one list, or indices are skipped.
  The final value of <code>numSamples</code> always becomes correct, however.</p>
  <p>However, when using the following custom implementation for <code>atomicInc_</code> instead, it seems to work correctly:</p>
  <pre><code>__device__ static inline uint32_t atomicInc_(uint32_t* ptr) {
      uint32_t value;
      do {
          value = *ptr;
      } while(atomicCAS(ptr, value, value + 1) != value);
      return value;
  }
  </code></pre>
  <p>Are the two <code>atomicInc</code> functions equivalent, and is it valid to use <code>atomicInc</code> that way to get unique indices?</p>
  <p>According the the <a href="https://docs.nvidia.com/cuda/cuda-c-programming-guide/#atomic-functions" rel="nofollow noreferrer">CUDA programming guide</a>, the atomic functions do not imply memory ordering constraints, and different threads can access the <code>numSamples</code> of different lists at the same time: could this cause it to fail?</p>
  <p><strong>Edit:</strong></p>
  <p>The real kernel looks like this:</p>
  <p>Basically there is a list of <em>spot blocks</em>, containing <em>spots</em>. Each <em>spot</em> has XY coordinates (<em>col</em>, <em>row</em>). The kernel needs to find, for each spot, the spots that are in a certain window (col/row difference) around it, and put them into a list in shared memory.</p>
  <p>The kernel is called with a fixed number of warps. A CUDA block corresponds to a group of <em>spot blocks</em>. (here 3) These are called the <em>local</em> spot blocks.</p>
  <p>First it takes the spots from the block's 3 spot blocks, and copies them into shared memory (<code>localSpots[]</code>).
  For this it uses one warp for each spot block, so that the spots can be read coalesced. Each thread in the warp is a spot in the local spot block.
  The spot block indices are here hardcoded (<code>blocks[]</code>).</p>
  <p>Then it goes through the <em>surrounding</em> spot blocks: These are all the spot blocks that may contain spots that are close enough to a spot in the <em>local spot blocks</em>. The surrounding spot block's indices are also hardcoded here (<code>sblock[]</code>).</p>
  <p>In this example it only uses the first warp for this, and traverses <code>sblocks[]</code> iteratively. Each thread in the warp is a spot in the surrounding spot block.
  It also iterates through the list of all the local spots. If the thread's spot is close enough to the local spot: It inserts it into the local spot's list, using <code>atomicInc</code> to get an index.</p>
  <p>When executed, the printf shows that for a given local spot (here the one with row=37, col=977), indices are sometimes repeated or skipped.</p>
  <p>The real code is more complex/optimized, but this code already has the problem. Here it also only runs one CUDA block.</p>
  <pre><code>#include &lt;assert.h&gt;
  #include &lt;stdio.h&gt;
  
  #define MAX_NUM_SPOTS_IN_WINDOW 80
  
  __global__ void Kernel(
      const uint16_t* blockNumSpotsBuffer,
      XGPU_SpotProcessingBlockSpotDataBuffers blockSpotsBuffers,
      size_t blockSpotsBuffersElementPitch,
      int2 unused1,
      int2 unused2,
      int unused3 ) {
      typedef unsigned int uint;
  
      if(blockIdx.x!=30 || blockIdx.y!=1) return;
  
      int window = 5;
  
      ASSERT(blockDim.x % WARP_SIZE == 0);
      ASSERT(blockDim.y == 1);
  
      uint numWarps = blockDim.x / WARP_SIZE;
      uint idxWarp = threadIdx.x / WARP_SIZE;
      int idxThreadInWarp = threadIdx.x % WARP_SIZE;
  
      struct Spot {
          int16_t row;
          int16_t col;
          volatile unsigned int numSamples;
          float signalSamples[MAX_NUM_SPOTS_IN_WINDOW];
      };
  
      __shared__ uint numLocalSpots;
      __shared__ Spot localSpots[3 * 32];
  
      numLocalSpots = 0;
  
      __syncthreads();
  
      ASSERT(numWarps &gt;= 3);
      int blocks[3] = {174, 222, 270};
      if(idxWarp &lt; 3) {
          uint spotBlockIdx = blocks[idxWarp];
          ASSERT(spotBlockIdx &lt; numSpotBlocks.x * numSpotBlocks.y);
  
          uint numSpots = blockNumSpotsBuffer[spotBlockIdx];
          ASSERT(numSpots &lt; WARP_SIZE);
  
          size_t inOffset = (spotBlockIdx * blockSpotsBuffersElementPitch) + idxThreadInWarp;
  
          uint outOffset;
          if(idxThreadInWarp == 0) outOffset = atomicAdd(&amp;numLocalSpots, numSpots);
          outOffset = __shfl_sync(0xffffffff, outOffset, 0, 32);
  
          if(idxThreadInWarp &lt; numSpots) {
              Spot* outSpot = &amp;localSpots[outOffset + idxThreadInWarp];
              outSpot-&gt;numSamples = 0;
  
              uint32_t coord = blockSpotsBuffers.coord[inOffset];
              UnpackCoordinates(coord, &amp;outSpot-&gt;row, &amp;outSpot-&gt;col);
          }
      }
  
  
  
      __syncthreads();
  
  
      int sblocks[] = { 29,30,31,77,78,79,125,126,127,173,174,175,221,222,223,269,270,271,317,318,319,365,366,367,413,414,415 };
      if(idxWarp == 0) for(int block = 0; block &lt; sizeof(sblocks)/sizeof(int); ++block) {
          uint spotBlockIdx = sblocks[block];
          ASSERT(spotBlockIdx &lt; numSpotBlocks.x * numSpotBlocks.y);
  
          uint numSpots = blockNumSpotsBuffer[spotBlockIdx];
          uint idxThreadInWarp = threadIdx.x % WARP_SIZE;
          if(idxThreadInWarp &gt;= numSpots) continue;
  
          size_t inOffset = (spotBlockIdx * blockSpotsBuffersElementPitch) + idxThreadInWarp;
  
          uint32_t coord = blockSpotsBuffers.coord[inOffset];
          if(coord == 0) return; // invalid surrounding spot
  
          int16_t row, col;
          UnpackCoordinates(coord, &amp;row, &amp;col);
  
          for(int idxLocalSpot = 0; idxLocalSpot &lt; numLocalSpots; ++idxLocalSpot) {
              Spot* localSpot = &amp;localSpots[idxLocalSpot];
  
              if(localSpot-&gt;row == 0 &amp;&amp; localSpot-&gt;col == 0) continue;
              if((abs(localSpot-&gt;row - row) &gt;= window) &amp;&amp; (abs(localSpot-&gt;col - col) &gt;= window)) continue;
  
              int index = atomicInc_block((unsigned int*)&amp;localSpot-&gt;numSamples, 0xffffffff);
              if(localSpot-&gt;row == 37 &amp;&amp; localSpot-&gt;col == 977) printf(&quot;%02d  &quot;, index); // &lt;-- sometimes indices are skipped or duplicated
  
              if(index &gt;= MAX_NUM_SPOTS_IN_WINDOW) continue; // index out of bounds, discard value for median calculation
              localSpot-&gt;signalSamples[index] = blockSpotsBuffers.signal[inOffset];
          }
      } }
  </code></pre>
  <p>Output looks like this:</p>
  <pre><code>00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  23                                                                                                                   
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24                                                                                                                 
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24        
  00  01  02  02  03  03  04  05  06  07  08  09  10  11  12  06  13  14  15  16  17  18  19  20  21        
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24        
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24        
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  23        
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24        
  00  01  02  03  04  05  06  07  08  09  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24    
  </code></pre>
  <p>Each line is the output of one execution (the kernel is run multiple times). It is expected that indices appear in different orders. But for example on the third-last line, index 23 is repeated.</p>
  <p>Using <code>atomicCAS</code> seems to fix it. Also using <code>__syncwarp()</code> between executions on the outer for-loop seems to fix it. But it is not clear why, and if that always fixes it.</p>
  <hr />
  <p><strong>Edit 2:</strong>
  This is a full program (main.cu) that shows the problem:</p>
  <p><a href="https://pastebin.com/cDqYmjGb" rel="nofollow noreferrer">https://pastebin.com/cDqYmjGb</a></p>
  <p>The CMakeLists.txt:</p>
  <p><a href="https://pastebin.com/iB9mbUJw" rel="nofollow noreferrer">https://pastebin.com/iB9mbUJw</a></p>
  <p>Must be compiled with -DCMAKE_BUILD_TYPE=Release.</p>
  <p>It produces this output:</p>
  <pre><code>00(0:00000221E40003E0)
  01(2:00000221E40003E0)
  02(7:00000221E40003E0)
  03(1:00000221E40003E0)
  03(2:00000221E40003E0)
  04(3:00000221E40003E0)
  04(1:00000221E40003E0)
  05(4:00000221E40003E0)
  06(6:00000221E40003E0)
  07(2:00000221E40003E0)
  08(3:00000221E40003E0)
  09(6:00000221E40003E0)
  10(3:00000221E40003E0)
  11(5:00000221E40003E0)
  12(0:00000221E40003E0)
  13(1:00000221E40003E0)
  14(3:00000221E40003E0)
  15(1:00000221E40003E0)
  16(0:00000221E40003E0)
  17(3:00000221E40003E0)
  18(0:00000221E40003E0)
  19(2:00000221E40003E0)
  20(4:00000221E40003E0)
  21(4:00000221E40003E0)
  22(1:00000221E40003E0)
  </code></pre>
  <p>For example the lines with 03 show that two threads (1 and 2), get the same result (3), after calling <code>atomicInc_block</code> on the same counter (at <code>0x00000221E40003E0</code>).</p>