title: How to avoid un-coalesced accesses in matrix multiplication CUDA kernel?
body: <p>I am learning CUDA with the book '<strong>Programming Massively Parallel Processors</strong>'. A practice problem from chapter 5 confuses me: </p>
  
  <blockquote>
    <p>For tiled matrix multiplication out of possible range of values for
    BLOCK_SIZE, for what values of BLOCK_SIZE will the kernel completely
    avoid un-coalesced accesses to global memory? (you only need to consider square blocks)</p>
  </blockquote>
  
  <p>On my understanding, BLOCK_SIZE does little to memory-coalescing. As long as threads within single warp access consecutive elements, we will have a coalesced accesses. I could not figure out where the kernel has un-coalesced accesses to global memory. Any hints from you guys? </p>
  
  <p>Here is the kernel's source codes:</p>
  
  <pre><code>#define COMMON_WIDTH 512
  #define ROW_LEFT 500 
  #define COL_RIGHT 250
  #define K 1000
  #define TILE_WIDTH 32
  __device__ int D_ROW_LEFT = ROW_LEFT;
  __device__ int D_COL_RIGHT = COL_RIGHT;
  __device__ int D_K = K;
  .....
  __global__
  void MatrixMatrixMultTiled(float *matrixLeft, float *matrixRight, float *output){
      __shared__  float sMatrixLeft[TILE_WIDTH][TILE_WIDTH];
      __shared__  float sMatrixRight[TILE_WIDTH][TILE_WIDTH];  
     int bx = blockIdx.x; int by = blockIdx.y;
     int tx = threadIdx.x; int ty = threadIdx.y;
     int col = bx * TILE_WIDTH + tx;
     int row = by * TILE_WIDTH + ty;
     float value = 0;
     for (int i = 0; i &lt; ceil(D_K/(float)TILE_WIDTH); ++i){
         if (row &lt; D_ROW_LEFT &amp;&amp; row * D_K + i * TILE_WIDTH  +tx &lt; D_K){
          sMatrixLeft[ty][tx]  = matrixLeft[row * D_K + i * TILE_WIDTH  +tx];
         }
         if (col &lt; D_COL_RIGHT &amp;&amp; (ty + i * TILE_WIDTH) * D_COL_RIGHT  + col &lt; D_K ){
          sMatrixRight[ty][tx] = matrixRight[(ty + i * TILE_WIDTH) * D_COL_RIGHT  + col];
         }
         __syncthreads();
         for (int j = 0; j &lt; TILE_WIDTH; j++){
             value += sMatrixLeft[ty][j] * sMatrixRight[j][tx]; 
         }
         __syncthreads();
     }
     if (row &lt; D_ROW_LEFT &amp;&amp; col &lt; D_COL_RIGHT ){
          output[row * D_COL_RIGHT + col] = value;
         }
  }
  </code></pre>