Shared-memory matmul

The tile

__shared__ float As[32][32], Bs[32][32];
int row = blockIdx.y * 32 + threadIdx.y;
int col = blockIdx.x * 32 + threadIdx.x;
float acc = 0.f;
for (int k0 = 0; k0 < N; k0 += 32) {
  As[threadIdx.y][threadIdx.x] = A[row * N + k0 + threadIdx.x];
  Bs[threadIdx.y][threadIdx.x] = B[(k0 + threadIdx.y) * N + col];
  __syncthreads();
  for (int k = 0; k < 32; ++k) acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
  __syncthreads();
}
C[row * N + col] = acc;

A 32 by 32 block is 1024 threads, the maximum, and each tile is 4 KiB, so the pair is 8 KiB per block. Both global loads are coalesced: threadIdx.x varies fastest and walks contiguous elements in each case.

What the tile buys

Each k-step loads 2 * 32 * 32 floats and performs 32 * 32 * 32 multiply-adds, so each element crosses the global boundary once and is read from shared 32 times. In the one-element-per-thread kernel every FMA reads two values from global; here traffic per FLOP falls by the tile width. The saving is the reuse, not the shared memory.

Two barriers, not one

The first barrier is between the load and the compute. The second is between the compute and the next load, and it is the one that gets dropped: a warp that finishes its inner loop can start overwriting As while another warp is still reading it. That version passes its first run and fails its third, which is the worst failure mode available.

Things that break

Missing the second barrier gives results that change with scheduling. Shared-memory bank conflicts depend on the pattern: the reads above are free, but a warp that reads down a column of a 32-float row puts every lane in one bank, so pad to [32][33] when the kernel needs that access. A static __shared__ array above 48 KiB is rejected at compile time; dynamic shared memory plus an opt-in through cudaFuncSetAttribute is the way past it. The default limit is architecture and toolkit dependent; check sharedMemPerBlock from your device query and the guide for your release.