The most obvious kernel assigns one thread to each element of C. The thread walks the k dimension, reading one row of A and one column of B, and accumulates the dot product. Written as a grid-stride loop so that any launch configuration covers an arbitrary matrix:
__global__ void sgemm_naive( int M, int N, int K, const float *A, const float *B, float *C ){ for ( int row = blockIdx.y*blockDim.y + threadIdx.y; row < M; row += gridDim.y*blockDim.y ) { for ( int col = blockIdx.x*blockDim.x + threadIdx.x; col < N; col += gridDim.x*blockDim.x ) { float acc = 0.f; for ( int k = 0; k < K; k++ ) acc += A[row*K + k] * B[k*N + col]; C[row*N + col] = acc; } }}sgemm_naive).This kernel reaches 886 GFLOP/s – about 7% of the card’s FP32 peak, and nowhere near cuBLAS. The reason is entirely in the memory traffic. Every one of the M×N threads reads a full row of A and a full column of B from global memory, so each element of A is fetched from DRAM once per output column and each element of B once per output row. The multiply-adds are cheap; the kernel spends almost all of its time waiting on memory. The arithmetic intensity that Section 17.1 promised is being discarded, because no data is reused.
Figure 17-1. The naive kernel computes each output C[i][j] from a full row of A and a full column of B, with no reuse between threads.