Developing a GPU Accelerated BLAS

I wrote these CUDA kernels to understand how BLAS accelerates linear algebra. Implementing these operations on a GPU involves distributing work across threads, accessing memory, and combining partial results. Writing these routines in CUDA gave me a way to connect the mathematics to the work the hardware actually performs.

BLAS, the Basic Linear Algebra Subprograms, provides the familiar operations underneath scientific computing and many machine-learning workloads. My goal was to understand the implementation choices that make those operations fast: how data is laid out, how work is divided, and where a value can be reused before another memory access is needed.

I worked through three routines: DOT, a vector reduction; GEMV, a matrix-vector product; and GEMM, a matrix-matrix product. The walkthrough follows my optimization process, including experiments that barely helped or made things slower.

The measurements in this article are saved Nsight Compute captures from an NVIDIA GeForce RTX 4060 using FP32. DOT uses 2²⁴ elements, GEMV uses a 4096 × 4096 matrix, and GEMM uses 2048 × 2048 matrices. These are historical kernel measurements, not repeated-run statistical estimates or end-to-end API timings. Supporting screenshots are expandable, and every screenshot can be enlarged. The animations play while they are on screen; use their controls to pause or scrub.

Profiler components

NVIDIA Nsight Compute is the kernel profiler used throughout this article. Its command-line tool, ncu, launches the program and collects GPU performance counters for selected kernel launches. ncu-ui opens the saved report as the tables and charts shown below. Counters describe what the hardware did, helping connect a code change to a possible bottleneck.

For example, from the BLAS repository, the project’s profiling workflow can be written as:

mkdir -p reports
ncu --kernel-name 'regex:gemm' --launch-skip 1 --launch-count 7 \
    --set full -o reports/gemm_example ./build/bench
ncu-ui reports/gemm_example.ncu-rep

The name filter selects matching kernels; the next flags skip one matching launch and collect up to seven. --set full requests a broad collection of metrics, and -o saves the report. This is an example invocation, not a reconstruction of every historical capture. See NVIDIA’s CLI guide for filtering and collection options.

The profiler reports several core components:

• Kernel and duration. The kernel name identifies the implementation being measured. Duration is its execution time, excluding allocation and host-transfer costs.
• Speed Of Light (SOL). Compute and memory throughput are expressed relative to hardware ceilings. The breakdown distinguishes activity in caches, shared memory, and DRAM; high utilization in one does not imply high utilization in the others.
• Memory chart. The chart shows traffic through L1/TEX, L2, and device memory. Requests count memory operations, transferred bytes measure traffic, and throughput measures the rate of transfer. Cache reuse can reduce DRAM traffic without changing the number of source-level loads.
• Launch, occupancy, and scheduler statistics. Blocks and threads describe how work is distributed. Occupancy measures resident warps, while eligible warps are ready to issue an instruction. Register use, shared-memory use, and stalls help explain the difference. NVIDIA’s profiling guide describes these metrics.

Comparison labels. “Current” is the selected result and “Baseline 1” is its comparison result. In the final GEMV comparison, Current is cuBLAS and Baseline 1 is the custom kernel. Percentage changes are relative to that selection. See NVIDIA’s baseline documentation.

Nsight Compute may replay a kernel to collect counters that cannot be gathered together. Profiling therefore takes longer than normal execution; timing the entire ncu command is not a kernel benchmark. Cache and clock controls also affect the measurement conditions. These reports help diagnose kernel behavior under comparable measurement conditions. Repeated unprofiled timings provide a separate performance check. See replay and measurement controls.

The mathematical model

Before optimizing a kernel, redundant work should be identified. Counting the required arithmetic and data movement provides a baseline for evaluating changes. Let \(F\) be the number of floating-point operations, \(Q\) the bytes crossing a chosen memory boundary, and \(t\) the kernel time. A multiply-add counts as two FLOPs. Arithmetic intensity and achieved performance are

With compute ceiling \(P_{\max}\) and bandwidth \(B\), executing the arithmetic takes at least \(F/P_{\max}\) seconds and moving the data takes at least \(Q/B\) seconds. Therefore,

This is the roofline model. Below the ridge point \(I_*\), the bandwidth ceiling is lower; above it, the compute ceiling is lower. Neither is a promise of performance. Contention, dependencies, and synchronization can keep a kernel well below both. NERSC’s roofline guide describes the measurement methodology.

The important qualification is which bytes we count. A source-level load may hit in L1 or L2 instead of reaching DRAM. Reducing load instructions can help even when DRAM traffic hardly changes. Conversely, a DRAM roofline cannot be built from a count of source-level loads without accounting for cache reuse.

Let \(s\) be bytes per element and let \(\delta_\beta=1\) when \(\beta\ne0\), otherwise zero. The ideal useful traffic counts each required input once and includes an output read only when the old output matters.

DOT stays near \(0.25\,\mathrm{FLOP/byte}\) in FP32. Large GEMV approaches \(0.5\,\mathrm{FLOP/byte}\). Square GEMM can increase its intensity with matrix size because each input participates in many output values. That reuse is the opportunity; the kernel still has to realize it.

The implementation also has to preserve the linear-algebra contract. Matrices are column-major, so \(A_{ij}=A[i+j\,ld]\). A leading dimension can exceed the logical row count. Transpose flags change the access pattern. When \(\beta=0\), an old NaN in the output must not leak into the result. These details stay relevant even as the kernels become more complicated.

I. DOT

The dot product is

The products can be computed independently, then reduced to a single sum. The reduction sequence follows the ideas in Mark Harris’s Optimizing Parallel Reduction in CUDA. The times below are measurements from this implementation, not numbers from that tutorial.

DOT 0. Atomic baseline

The first implementation lets each thread multiply a pair of elements and atomically add the result to one global scalar. That is easy to write, but every update competes for the same destination. At \(n=2^{24}\), we request roughly 16.8 million updates to one answer.

Every product reaches the same atomic destination:

__global__ void dot_kernel(const Vector x, const Vector y, real_t* result) {
    const real_t* __restrict__ xp = x.data;
    const real_t* __restrict__ yp = y.data;
    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i+= (long)gridDim.x * blockDim.x) {
	atomicAdd(result, xp[i * x.inc] * yp[i * y.inc]);
    }
}

Read the baseline memory and utilization evidence

DOT 1. Fixed block size

I replaced the occupancy API’s 768-thread block with a fixed 256-thread block. Achieved occupancy rose from about 64% to 77%, but the contended update was still there.

The launch uses a fixed block size of 256 threads. The snippet isolates this change from the reduction improvements that follow.

#define DOT_BLOCK_SIZE 256

    const int grid = (x.n + DOT_BLOCK_SIZE - 1) / DOT_BLOCK_SIZE;
    real_t* d_result = nullptr;
    CUDA_ERROR_CHECK(cudaMalloc(&d_result, sizeof(real_t)));
    CUDA_ERROR_CHECK(cudaMemset(d_result, 0, sizeof(real_t)));
    dot_kernel<<<grid, DOT_BLOCK_SIZE>>>(dx, dy, d_result);

This was a useful failed experiment. Occupancy tells me how many warps are resident, not how much independent work they can complete. More threads waiting for the same destination do not solve the serialization problem.

DOT 2. Block reduction

Each thread now accumulates a private sum, then the block combines its partial sums in shared memory. Only thread zero performs a global atomic update. With 256 threads per block, the initial configuration reduces the number of atomic updates by a factor of 256.

Each thread accumulates a private sum, then the block reduces those sums in shared memory and issues one atomic update. The full diff below also includes the contiguous-thread indexing explained in the next step.

    real_t sum = 0;
    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i += (long)gridDim.x * blockDim.x) {
	sum += xp[i * x.inc] * yp[i * y.inc];
    }
    sdata[tid] = sum;
    __syncthreads();

    if (tid == 0) {
	atomicAdd(result, sdata[0]);
    }

@@ -2,6 +2,8 @@

 #include <cmath>

+#define DOT_BLOCK_SIZE 256
+
 static long span_of(int n, int inc) {
     return n > 0 ? 1 + (long)(n - 1) * (inc < 0 ? -inc : inc) : 0;
 }
@@ -158,27 +160,42 @@ void swap(Vector& x, Vector& y) {
 __global__ void dot_kernel(const Vector x, const Vector y, real_t* result) {
     const real_t* __restrict__ xp = x.data;
     const real_t* __restrict__ yp = y.data;
-    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i+= (long)gridDim.x * blockDim.x) {
-	atomicAdd(result, xp[i * x.inc] * yp[i * y.inc]);
+    __shared__ real_t sdata[DOT_BLOCK_SIZE];
+    unsigned int tid = threadIdx.x;
+
+    real_t sum = 0;
+    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i += (long)gridDim.x * blockDim.x) {
+	sum += xp[i * x.inc] * yp[i * y.inc];
+    }
+    sdata[tid] = sum;
+    __syncthreads();
+
+    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
+	int index = 2 * s * tid;
+	if (index < blockDim.x) {
+	    sdata[index] += sdata[index + s];
+	}
+	__syncthreads();
+    }
+
+    if (tid == 0) {
+	atomicAdd(result, sdata[0]);
     }
 }

 real_t dot(const Vector& x, const Vector& y) {
     if (x.n <= 0 || y.n <= 0) return 0;

-    int min_grid = 0, block = 0;
-    GET_MAX_POTENTIAL_BLOCKS_SIZE(min_grid, block, dot_kernel);
-
     real_t* sx = nullptr;
     real_t* sy = nullptr;
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = (x.n + block - 1) / block;
+    const int grid = (x.n + DOT_BLOCK_SIZE - 1) / DOT_BLOCK_SIZE;
     real_t* d_result = nullptr;
     CUDA_ERROR_CHECK(cudaMalloc(&d_result, sizeof(real_t)));
     CUDA_ERROR_CHECK(cudaMemset(d_result, 0, sizeof(real_t)));
-    dot_kernel<<<grid, block>>>(dx, dy, d_result);
+    dot_kernel<<<grid, DOT_BLOCK_SIZE>>>(dx, dy, d_result);
     CUDA_ERROR_CHECK(cudaGetLastError());
     real_t result = 0;
     CUDA_ERROR_CHECK(cudaMemcpy(&result, d_result, sizeof(real_t), cudaMemcpyDeviceToHost));

The input bytes have barely changed. The expensive coordination has: instead of one global update per product, there is one update per block. The arithmetic intensity estimate did not identify this contention cost, but it still tells us where the kernel might end up after removing it.

DOT 3. Contiguous active threads

The first tree uses interleaved participation: threads whose indices are multiples of the current stride remain active. Reindexing the tree keeps the active threads contiguous instead of scattering them across each warp. Late stages still have fewer active lanes; this change reduces unnecessary divergence earlier in the tree.

The code below keeps participating threads contiguous. The earlier interleaved tree illustrates the experiment; the full kernel diff above includes this indexing improvement.

    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
	int index = 2 * s * tid;
	if (index < blockDim.x) {
	    sdata[index] += sdata[index + s];
	}
	__syncthreads();
    }

Inspect the newly exposed shared-memory bottleneck

DOT 4. Sequential shared-memory addressing

I reversed the tree traversal. At stride \(s\), thread \(i<s\) adds shared entries \(i\) and \(i+s\). Adjacent active threads now access adjacent shared words, reducing the bank conflicts found in the previous profile.

for (unsigned stride = blockDim.x / 2; stride > 0; stride >>= 1) {
    if (tid < stride)
        partial[tid] += partial[tid + stride];
    __syncthreads();
}

Every thread reaches the barrier, including threads whose arithmetic is inactive.

The reduction now starts with a large stride and halves it at each stage, keeping shared-memory accesses adjacent.

@@ -170,10 +170,9 @@ __global__ void dot_kernel(const Vector x, const Vector y, real_t* result) {
     sdata[tid] = sum;
     __syncthreads();

-    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
-	int index = 2 * s * tid;
-	if (index < blockDim.x) {
-	    sdata[index] += sdata[index + s];
+    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
+	if (tid < s) {
+	    sdata[tid] += sdata[tid + s];
 	}
 	__syncthreads();
     }

DOT 5. More elements per thread

The next change halves the grid. Each thread processes two input positions in its grid-stride loop before joining the reduction tree. A block now amortizes its tree, barriers, and final atomic update over twice as many products.

Only the grid calculation changes. The existing grid-stride loop makes a second pass at this input size without manually unrolling another load.

@@ -190,7 +190,7 @@ real_t dot(const Vector& x, const Vector& y) {
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = (x.n + DOT_BLOCK_SIZE - 1) / DOT_BLOCK_SIZE;
+    const int grid = (x.n + 2L * DOT_BLOCK_SIZE - 1) / (2L * DOT_BLOCK_SIZE);
     real_t* d_result = nullptr;
     CUDA_ERROR_CHECK(cudaMalloc(&d_result, sizeof(real_t)));
     CUDA_ERROR_CHECK(cudaMemset(d_result, 0, sizeof(real_t)));

See the duration and throughput measurements for two elements per thread

DOT 6. Warp reduction

Finally, each warp combines its partial sums with shuffle instructions. One lane per warp writes a partial result to shared memory; the first warp combines those results after a block barrier. The shared array shrinks from one entry per thread to one entry per warp.

The kernel stores one shared partial per warp and uses shuffle instructions for both reduction stages. The diff below shows the synchronization and final atomic update.

@@ -160,25 +160,32 @@ void swap(Vector& x, Vector& y) {
 __global__ void dot_kernel(const Vector x, const Vector y, real_t* result) {
     const real_t* __restrict__ xp = x.data;
     const real_t* __restrict__ yp = y.data;
-    __shared__ real_t sdata[DOT_BLOCK_SIZE];
-    unsigned int tid = threadIdx.x;
+    __shared__ real_t sdata[DOT_BLOCK_SIZE / 32];
+    unsigned int warp = threadIdx.x / 32;
+    unsigned int lane = threadIdx.x % 32;

     real_t sum = 0;
     for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i += (long)gridDim.x * blockDim.x) {
 	sum += xp[i * x.inc] * yp[i * y.inc];
     }
-    sdata[tid] = sum;
+    sum += __shfl_down_sync(0xffffffff, sum, 16);
+    sum += __shfl_down_sync(0xffffffff, sum, 8);
+    sum += __shfl_down_sync(0xffffffff, sum, 4);
+    sum += __shfl_down_sync(0xffffffff, sum, 2);
+    sum += __shfl_down_sync(0xffffffff, sum, 1);
+    if (lane == 0)
+	sdata[warp] = sum;
     __syncthreads();

-    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
-	if (tid < s) {
-	    sdata[tid] += sdata[tid + s];
-	}
-	__syncthreads();
-    }
-
-    if (tid == 0) {
-	atomicAdd(result, sdata[0]);
+    if (warp == 0) {
+	sum = lane < DOT_BLOCK_SIZE / 32 ? sdata[lane] : 0;
+	sum += __shfl_down_sync(0xffffffff, sum, 16);
+	sum += __shfl_down_sync(0xffffffff, sum, 8);
+	sum += __shfl_down_sync(0xffffffff, sum, 4);
+	sum += __shfl_down_sync(0xffffffff, sum, 2);
+	sum += __shfl_down_sync(0xffffffff, sum, 1);
+	if (lane == 0)
+	    atomicAdd(result, sum);
     }
 }

Using useful input bytes, the effective bandwidth is

The small final gain makes sense: the preceding version was already moving data efficiently. Across the entire sequence, DOT improved by about \(61.8\times\). Floating-point addition order changed along the way, so correctness checks use numerical tolerances rather than demand identical bits across different reduction trees.

Complete DOT implementation

Both the warp-shuffle kernel and the host entry point are shown below. The wrapper stages the inputs and copies the scalar result back; those transfers are outside the kernel timings above. The shared types, constants, and memory helpers are listed at the end.

__global__ void dot_kernel(const Vector x, const Vector y, real_t* result) {
    const real_t* __restrict__ xp = x.data;
    const real_t* __restrict__ yp = y.data;
    __shared__ real_t sdata[BLOCK_SIZE / 32];
    unsigned int warp = threadIdx.x / 32;
    unsigned int lane = threadIdx.x % 32;

    real_t sum = 0;
    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < x.n; i += (long)gridDim.x * blockDim.x) {
	sum += xp[i * x.inc] * yp[i * y.inc];
    }
    sum += __shfl_down_sync(0xffffffff, sum, 16);
    sum += __shfl_down_sync(0xffffffff, sum, 8);
    sum += __shfl_down_sync(0xffffffff, sum, 4);
    sum += __shfl_down_sync(0xffffffff, sum, 2);
    sum += __shfl_down_sync(0xffffffff, sum, 1);
    if (lane == 0)
	sdata[warp] = sum;
    __syncthreads();

    if (warp == 0) {
	sum = lane < BLOCK_SIZE / 32 ? sdata[lane] : 0;
	sum += __shfl_down_sync(0xffffffff, sum, 16);
	sum += __shfl_down_sync(0xffffffff, sum, 8);
	sum += __shfl_down_sync(0xffffffff, sum, 4);
	sum += __shfl_down_sync(0xffffffff, sum, 2);
	sum += __shfl_down_sync(0xffffffff, sum, 1);
	if (lane == 0)
	    atomicAdd(result, sum);
    }
}

real_t dot(const Vector& x, const Vector& y) {
    if (x.n <= 0 || y.n <= 0) return 0;

    real_t* sx = nullptr;
    real_t* sy = nullptr;
    const Vector dx = stage_vector(x, sx);
    const Vector dy = stage_vector(y, sy);

    const int grid = (x.n + 2L * BLOCK_SIZE - 1) / (2L * BLOCK_SIZE);
    real_t* d_result = nullptr;
    CUDA_ERROR_CHECK(cudaMalloc(&d_result, sizeof(real_t)));
    CUDA_ERROR_CHECK(cudaMemset(d_result, 0, sizeof(real_t)));
    dot_kernel<<<grid, BLOCK_SIZE>>>(dx, dy, d_result);
    CUDA_ERROR_CHECK(cudaGetLastError());
    real_t result = 0;
    CUDA_ERROR_CHECK(cudaMemcpy(&result, d_result, sizeof(real_t), cudaMemcpyDeviceToHost));
    CUDA_ERROR_CHECK(cudaFree(d_result));
    release(sx);
    release(sy);
    return result;
}

II. GEMV

GEMV computes

Each output is a dot product, but the matrix dominates input traffic. Reusing \(\boldsymbol{x}\) helps only so much: each matrix element still contributes to just one multiply-add. For large dimensions,

GEMV 0. Baseline launch

The first kernel launched only eight blocks of 256 threads for 4096 output rows. On a GPU with 24 streaming multiprocessors, eight blocks cannot occupy every multiprocessor at once. Threads loop over more rows, but that does not create more independent blocks.

A fixed block size of 256 produces just eight blocks for 4096 rows, replacing the occupancy-based launch.

@@ -309,11 +307,11 @@ void gemv(const char* trans, real_t alpha, const Matrix& A, const Vector& x, rea
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = (m + block - 1) / block;
+    const int grid = (m + 2L * BLOCK_SIZE - 1) / (2L * BLOCK_SIZE);
     if (notrans) {
-	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
+	gemv_kernel<<<grid, BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
     } else {
-	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
+	gemv_kernel<<<grid, BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
     }
     CUDA_ERROR_CHECK(cudaGetLastError());
     unstage_vector(y, sy);

GEMV reference. cuBLAS

The reference kernel took \(279.62\,\mu\mathrm{s}\) in the same profiling report, roughly \(2.76\times\) faster. It moved a similar amount of matrix data but approached the DRAM throughput ceiling. That made a bandwidth-oriented implementation a plausible target.

The benchmark calls cuBLAS on device buffers with the same matrix dimensions and scalar parameters.

    const real_t alpha = 1, beta = 0;
    for (int i = 0; i < ITERS; ++i)
        CUBLAS_GEMV(handle, CUBLAS_OP_N, GEMV_N, GEMV_N, &alpha, dA, GEMV_N, dv, 1, &beta, dw, 1);
    CUDA_ERROR_CHECK(cudaDeviceSynchronize());

Compare how much matrix data the two GEMV kernels move

GEMV 1. Shared-memory vector cache

I cooperatively loaded chunks of \(\boldsymbol{x}\) into shared memory so the threads could reuse them. This added barriers around each chunk, while keeping the same eight-block launch.

The code caches the input vector in shared memory. The full implementation uses 128 blocks of 512 threads; the screenshot measures the earlier eight-block configuration.

__device__ real_t row_dot(const Matrix& A, const Vector& x, long i, long m, long len, Access at)
{
  __shared__ real_t xs[GEMV_BLOCK_SIZE];
  real_t sum = 0;
  for (long tile = 0; tile < len; tile += GEMV_BLOCK_SIZE) {
      long j = tile + threadIdx.x;
      if (j < len)
	  xs[threadIdx.x] = x.data[j * x.inc];
      __syncthreads();
      const long lim = (len - tile < GEMV_BLOCK_SIZE) ? len - tile : GEMV_BLOCK_SIZE;
      if (i < m) {
	for (long k = 0; k < lim; k++) {
	    sum += at(A, i, tile+k) * xs[k];
	}
      }
      __syncthreads();
  }
  return sum;
}

This historical excerpt is shown unchanged. Its caller guards row_dot with i < m even though the helper contains block barriers. The displayed square benchmark avoids a partial row block; this is not a general tail-safe implementation. The later row/column split moves the barriers outside that guard.

Threads already read the same vector entries in a cache-friendly pattern. The extra staging did not save enough expensive traffic to pay for its synchronization. Meanwhile, the larger issue identified by the first profile remained unchanged.

GEMV 2. Launch geometry

I kept the shared-vector implementation and changed the number of blocks and threads. This isolates an important distinction: larger blocks are not the same as a larger supply of independent blocks.

The listing uses 128 blocks of 512 threads. The sweep below also tests other launch configurations, including the faster 512 × 256 trial.

@@ -307,11 +323,11 @@ void gemv(const char* trans, real_t alpha, const Matrix& A, const Vector& x, rea
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = (m + 2L * BLOCK_SIZE - 1) / (2L * BLOCK_SIZE);
+    const int grid = 128;
     if (notrans) {
-	gemv_kernel<<<grid, BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
+	gemv_kernel<<<grid, GEMV_BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
     } else {
-	gemv_kernel<<<grid, BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
+	gemv_kernel<<<grid, GEMV_BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
     }
     CUDA_ERROR_CHECK(cudaGetLastError());
     unstage_vector(y, sy);

Compare the larger-block launch and its cuBLAS reference

The best sweep result is about \(1.86\times\) faster than the shared-vector version’s \(794.30\,\mu\mathrm{s}\). It is also faster than the \(128\times512\) configuration retained at that intermediate point in the project. Both measurements show how much launch geometry affects performance.

GEMV 3. Row and column splitting

The final change uses a two-dimensional block. The x coordinate selects one of 32 adjacent rows, while the y coordinate splits the columns into 16 chunks. Each thread computes a partial sum for its row and column chunk; those partials are combined before writing the output.

This is not one warp reading across a row. At a fixed column, adjacent lanes read adjacent rows, matching column-major storage. Splitting the reduction across column chunks also shortens the serial dependency chain within each thread.

Each block covers 32 rows and 16 column splits. Every split writes a partial sum; after synchronization, split zero combines the partials for each row.

    __shared__ real_t partial[GEMV_SPLITS][GEMV_ROWS];
    real_t* __restrict__ yp = y.data;
    const long chunk = (n + GEMV_SPLITS - 1) / GEMV_SPLITS;
    const long c0 = threadIdx.y * chunk;
    const long c1 = (c0 + chunk < n) ? c0 + chunk : n;
    for (long base = blockIdx.x * (long)GEMV_ROWS; base < m; base += (long)gridDim.x * GEMV_ROWS) {
	const long i = base + threadIdx.x;
	partial[threadIdx.y][threadIdx.x] = (i < m) ? row_dot(A, x, i, c0, c1, at) : real_t(0);
	__syncthreads();
	if (threadIdx.y == 0 && i < m) {
	    real_t sum = 0;
	    for (int t = 0; t < GEMV_SPLITS; t++)
		sum += partial[t][threadIdx.x];
	    yp[i * y.inc] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * yp[i * y.inc]);
	}
	__syncthreads();
    }
}

@@ -3,7 +3,9 @@
 #include <cmath>

 #define BLOCK_SIZE 256
-#define GEMV_BLOCK_SIZE 512
+#define GEMV_BLOCK_SIZE BLOCK_SIZE
+#define GEMV_ROWS 32
+#define GEMV_SPLITS 16

 static long span_of(int n, int inc) {
     return n > 0 ? 1 + (long)(n - 1) * (inc < 0 ? -inc : inc) : 0;
@@ -275,35 +277,32 @@ real_t asum(const Vector& x) {


 template <typename Access>
-__device__ real_t row_dot(const Matrix& A, const Vector& x, long i, long m, long len, Access at)
+__device__ real_t row_dot(const Matrix& A, const Vector& x, long i, long c0, long c1, Access at)
 {
-  __shared__ real_t xs[GEMV_BLOCK_SIZE];
   real_t sum = 0;
-  for (long tile = 0; tile < len; tile += GEMV_BLOCK_SIZE) {
-      long j = tile + threadIdx.x;
-      if (j < len)
-	  xs[threadIdx.x] = x.data[j * x.inc];
-      __syncthreads();
-      const long lim = (len - tile < GEMV_BLOCK_SIZE) ? len - tile : GEMV_BLOCK_SIZE;
-      if (i < m) {
-	for (long k = 0; k < lim; k++) {
-	    sum += at(A, i, tile+k) * xs[k];
-	}
-      }
-      __syncthreads();
-  }
+  for (long j = c0; j < c1; j++)
+      sum += at(A, i, j) * x.data[j * x.inc];
   return sum;
 }

 template <typename Access>
 __global__ void gemv_kernel(real_t alpha, const Matrix A, const Vector x, real_t beta, Vector y, long m, long n, Access at) {
+    __shared__ real_t partial[GEMV_SPLITS][GEMV_ROWS];
     real_t* __restrict__ yp = y.data;
-    for (long base = blockIdx.x * (long)blockDim.x; base < m; base += (long)gridDim.x * blockDim.x) {
+    const long chunk = (n + GEMV_SPLITS - 1) / GEMV_SPLITS;
+    const long c0 = threadIdx.y * chunk;
+    const long c1 = (c0 + chunk < n) ? c0 + chunk : n;
+    for (long base = blockIdx.x * (long)GEMV_ROWS; base < m; base += (long)gridDim.x * GEMV_ROWS) {
 	const long i = base + threadIdx.x;
-	if (i < m) {
-	    real_t sum = row_dot(A, x, i, m, n, at);
+	partial[threadIdx.y][threadIdx.x] = (i < m) ? row_dot(A, x, i, c0, c1, at) : real_t(0);
+	__syncthreads();
+	if (threadIdx.y == 0 && i < m) {
+	    real_t sum = 0;
+	    for (int t = 0; t < GEMV_SPLITS; t++)
+		sum += partial[t][threadIdx.x];
 	    yp[i * y.inc] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * yp[i * y.inc]);
 	}
+	__syncthreads();
     }
 }

@@ -323,11 +322,12 @@ void gemv(const char* trans, real_t alpha, const Matrix& A, const Vector& x, rea
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = 128;
+    const dim3 block(GEMV_ROWS, GEMV_SPLITS);
+    const int grid = (m + GEMV_ROWS - 1) / GEMV_ROWS;
     if (notrans) {
-	gemv_kernel<<<grid, GEMV_BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
+	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
     } else {
-	gemv_kernel<<<grid, GEMV_BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
+	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
     }
     CUDA_ERROR_CHECK(cudaGetLastError());
     unstage_vector(y, sy);
@@ -339,7 +339,7 @@ template <typename Access>
 __global__ void symv_kernel(real_t alpha, const Matrix A, const Vector x, real_t beta, const Vector y, Access at) {
     real_t* __restrict__ yp = y.data;
     for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < A.rows; i+= (long)gridDim.x * blockDim.x) {
-	real_t sum = row_dot(A, x, i, A.rows, A.cols, at);
+	real_t sum = row_dot(A, x, i, 0, A.cols, at);
 	yp[i * y.inc] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * yp[i * y.inc]);
     }
 }
@@ -358,7 +358,7 @@ void symv(const char* uplo, real_t alpha, const Matrix& A, const Vector& x, real
     const Vector dx = stage_vector(x, sx);
     const Vector dy = stage_vector(y, sy);

-    const int grid = 128;
+    const int grid = 512;
     if (upper) {
 	symv_kernel<<<grid, GEMV_BLOCK_SIZE>>>(alpha, dA, dx, beta, dy, UpperAt());
     } else {

Verify that the final GEMV comparison moves similar data

This capture shows comparable kernel performance, not a statistically established win over cuBLAS. Relative to the original custom kernel, the recorded improvement is about \(2.78\times\). The successful change makes the matrix stream easier to serve and exposes more independent reduction work; it does not change GEMV’s limited arithmetic intensity.

Complete GEMV implementation

This includes the complete partial-row helper, the kernel that combines column chunks, and the host entry point with transpose dispatch. The shared types, constants, and memory helpers are listed at the end.

template <typename Access>
__device__ real_t row_dot(const Matrix& A, const Vector& x, long i, long c0, long c1, Access at)
{
  real_t sum = 0;
  for (long j = c0; j < c1; j++)
      sum += at(A, i, j) * x.data[j * x.inc];
  return sum;
}

template <typename Access>
__global__ void gemv_kernel(real_t alpha, const Matrix A, const Vector x, real_t beta, Vector y, long m, long n, Access at) {
    __shared__ real_t partial[GEMV_SPLITS][GEMV_ROWS];
    real_t* __restrict__ yp = y.data;
    const long chunk = (n + GEMV_SPLITS - 1) / GEMV_SPLITS;
    const long c0 = threadIdx.y * chunk;
    const long c1 = (c0 + chunk < n) ? c0 + chunk : n;
    for (long base = blockIdx.x * (long)GEMV_ROWS; base < m; base += (long)gridDim.x * GEMV_ROWS) {
	const long i = base + threadIdx.x;
	partial[threadIdx.y][threadIdx.x] = (i < m) ? row_dot(A, x, i, c0, c1, at) : real_t(0);
	__syncthreads();
	if (threadIdx.y == 0 && i < m) {
	    real_t sum = 0;
	    for (int t = 0; t < GEMV_SPLITS; t++)
		sum += partial[t][threadIdx.x];
	    yp[i * y.inc] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * yp[i * y.inc]);
	}
	__syncthreads();
    }
}


void gemv(const char* trans, real_t alpha, const Matrix& A, const Vector& x, real_t beta, Vector& y) {
    if (A.rows <= 0 || A.cols <= 0 || (alpha == real_t(0) && beta == real_t(1))) return;

    const bool notrans = (trans[0] == 'N' || trans[0] == 'n');
    const long m = notrans ? A.rows : A.cols;
    const long n = notrans ? A.cols : A.rows;


    real_t* sa = nullptr;
    real_t* sx = nullptr;
    real_t* sy = nullptr;
    const Matrix dA = stage_matrix(A, sa);
    const Vector dx = stage_vector(x, sx);
    const Vector dy = stage_vector(y, sy);

    const dim3 block(GEMV_ROWS, GEMV_SPLITS);
    const int grid = (m + GEMV_ROWS - 1) / GEMV_ROWS;
    if (notrans) {
	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, NoTransAt());
    } else {
	gemv_kernel<<<grid, block>>>(alpha, dA, dx, beta, dy, m, n, TransAt());
    }
    CUDA_ERROR_CHECK(cudaGetLastError());
    unstage_vector(y, sy);
    release(sx);
    release(sa);
}

III. GEMM

GEMM computes

Unlike GEMV, each value of \(A\) can contribute to many columns of \(C\), and each value of \(B\) to many rows. For square matrices of order \(N\), the ideal useful traffic model gives

At \(N=2048\) and \(s=4\), that is approximately \(341\,\mathrm{FLOP/byte}\) when \(\beta=0\). Achieving this input reuse is the implementation problem.

GEMM 0. Baseline

The first kernel assigns one \(C_{ij}\) to each thread and loops over the reduction dimension. For each multiply-add, it issues a load from \(A\) and a load from \(B\). Ignoring output traffic, the source-level input-load intensity is only

This is a load-request model, not a measured DRAM intensity: caches can reuse some of those values. However, the instructions requesting them still have to execute.

The baseline loads both inputs directly from global memory:

template <typename AccessA, typename AccessB>
__global__ void gemm_kernel(real_t alpha, const Matrix A, const Matrix B, real_t beta, Matrix C, long m, long n, long k, AccessA at_a, AccessB at_b) {
    real_t* __restrict__ cp = C.data;
    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < m; i += (long)gridDim.x * blockDim.x) {
	for (long j = blockIdx.y * (long)blockDim.y + threadIdx.y; j < n; j += (long)gridDim.y * blockDim.y) {
	    real_t sum = 0;
	    for (long l = 0; l < k; l++) {
		sum += at_a(A, i, l) * at_b(B, l, j);
	    }
	    cp[i + j * C.ld] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
	}
    }
}

GEMM 1. Shared-memory tiling

Each block now stages tiles of \(A\) and \(B\) in shared memory. Threads cooperate on the loads, synchronize, and reuse the staged values for their output elements. The first tiled version uses a \(16\times16\) output tile with one output per thread.

For a \(B_M\times B_N\) output tile and reduction depth \(B_K\),

For a square tile of width \(T\), this becomes \(T/s\). With \(T=16\) and FP32, the stage’s input-load intensity rises from roughly \(0.25\) to \(4\,\mathrm{FLOP/byte}\). Output traffic and edge waste are omitted in this stage-level estimate.

Two 16 × 16 shared tiles hold the inputs. Each stage loads guarded values, synchronizes, reuses the tiles, and synchronizes again before overwriting them.

@@ -373,14 +374,24 @@ void symv(const char* uplo, real_t alpha, const Matrix& A, const Vector& x, real

 template <typename AccessA, typename AccessB>
 __global__ void gemm_kernel(real_t alpha, const Matrix A, const Matrix B, real_t beta, Matrix C, long m, long n, long k, AccessA at_a, AccessB at_b) {
+    __shared__ real_t as[GEMM_TILE][GEMM_TILE];
+    __shared__ real_t bs[GEMM_TILE][GEMM_TILE];
     real_t* __restrict__ cp = C.data;
-    for (long i = blockIdx.x * (long)blockDim.x + threadIdx.x; i < m; i += (long)gridDim.x * blockDim.x) {
-	for (long j = blockIdx.y * (long)blockDim.y + threadIdx.y; j < n; j += (long)gridDim.y * blockDim.y) {
+    for (long ib = blockIdx.x * (long)GEMM_TILE; ib < m; ib += (long)gridDim.x * GEMM_TILE) {
+	for (long jb = blockIdx.y * (long)GEMM_TILE; jb < n; jb += (long)gridDim.y * GEMM_TILE) {
+	    const long i = ib + threadIdx.x;
+	    const long j = jb + threadIdx.y;
 	    real_t sum = 0;
-	    for (long l = 0; l < k; l++) {
-		sum += at_a(A, i, l) * at_b(B, l, j);
+	    for (long t = 0; t < k; t += GEMM_TILE) {
+		as[threadIdx.y][threadIdx.x] = (i < m && t + threadIdx.y < k) ? at_a(A, i, t + threadIdx.y) : real_t(0);
+		bs[threadIdx.y][threadIdx.x] = (t + threadIdx.x < k && j < n) ? at_b(B, t + threadIdx.x, j) : real_t(0);
+		__syncthreads();
+		for (int l = 0; l < GEMM_TILE; l++)
+		    sum += as[l][threadIdx.x] * bs[threadIdx.y][l];
+		__syncthreads();
 	    }
-	    cp[i + j * C.ld] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
+	    if (i < m && j < n)
+		cp[i + j * C.ld] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
 	}
     }
 }

Why did 16-fold reuse produce only a modest time improvement?

The runtime did not improve by 16 times because the eliminated requests were not all DRAM misses. Caches already provided some reuse. The kernel also now executes many shared-memory loads and barriers. Shared tiling creates the foundation for reuse, but one output per thread still consumes two shared operands for each multiply-add.

GEMM 2. 4 × 4 register tile

Each thread now owns 16 output values. At one reduction index, it loads four values from the A tile and four from the B tile into registers, then forms their outer product.

For an \(R_M\times R_N\) register tile, the scalar shared operands per FMA are

A \(1\times1\) tile needs two operands per FMA; a \(4\times4\) tile needs only one half. The block still has 256 threads, so its output footprint grows to \(64\times64\).

A 4 × 4 accumulator array expands each block’s output to 64 × 64. Threads share the input loads across this larger tile, then reuse the values in the inner loop below.

		for (int l = 0; l < GEMM_TILE; l++) {
		    real_t ra[GEMM_THREAD_TILE];
		    real_t rb[GEMM_THREAD_TILE];
		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
			ra[u] = as[l][threadIdx.x + u * GEMM_TILE];
		    for (int v = 0; v < GEMM_THREAD_TILE; v++)
			rb[v] = bs[threadIdx.y + v * GEMM_TILE][l];
		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
			for (int v = 0; v < GEMM_THREAD_TILE; v++)
			    acc[u][v] += ra[u] * rb[v];
		}

@@ -7,6 +7,8 @@
 #define GEMV_ROWS 32
 #define GEMV_SPLITS 16
 #define GEMM_TILE 16
+#define GEMM_THREAD_TILE 4
+#define GEMM_BLOCK_TILE (GEMM_TILE * GEMM_THREAD_TILE)

 static long span_of(int n, int inc) {
     return n > 0 ? 1 + (long)(n - 1) * (inc < 0 ? -inc : inc) : 0;
@@ -374,24 +376,45 @@ void symv(const char* uplo, real_t alpha, const Matrix& A, const Vector& x, real

 template <typename AccessA, typename AccessB>
 __global__ void gemm_kernel(real_t alpha, const Matrix A, const Matrix B, real_t beta, Matrix C, long m, long n, long k, AccessA at_a, AccessB at_b) {
-    __shared__ real_t as[GEMM_TILE][GEMM_TILE];
-    __shared__ real_t bs[GEMM_TILE][GEMM_TILE];
+    __shared__ real_t as[GEMM_TILE][GEMM_BLOCK_TILE];
+    __shared__ real_t bs[GEMM_BLOCK_TILE][GEMM_TILE];
     real_t* __restrict__ cp = C.data;
-    for (long ib = blockIdx.x * (long)GEMM_TILE; ib < m; ib += (long)gridDim.x * GEMM_TILE) {
-	for (long jb = blockIdx.y * (long)GEMM_TILE; jb < n; jb += (long)gridDim.y * GEMM_TILE) {
-	    const long i = ib + threadIdx.x;
-	    const long j = jb + threadIdx.y;
-	    real_t sum = 0;
+    const int tid = threadIdx.y * GEMM_TILE + threadIdx.x;
+    for (long ib = blockIdx.x * (long)GEMM_BLOCK_TILE; ib < m; ib += (long)gridDim.x * GEMM_BLOCK_TILE) {
+	for (long jb = blockIdx.y * (long)GEMM_BLOCK_TILE; jb < n; jb += (long)gridDim.y * GEMM_BLOCK_TILE) {
+	    real_t acc[GEMM_THREAD_TILE][GEMM_THREAD_TILE] = {};
 	    for (long t = 0; t < k; t += GEMM_TILE) {
-		as[threadIdx.y][threadIdx.x] = (i < m && t + threadIdx.y < k) ? at_a(A, i, t + threadIdx.y) : real_t(0);
-		bs[threadIdx.y][threadIdx.x] = (t + threadIdx.x < k && j < n) ? at_b(B, t + threadIdx.x, j) : real_t(0);
+		for (int q = 0; q < GEMM_THREAD_TILE; q++) {
+		    const int e = tid + q * GEMM_TILE * GEMM_TILE;
+		    const int ar = e % GEMM_BLOCK_TILE;
+		    const int ak = e / GEMM_BLOCK_TILE;
+		    as[ak][ar] = (ib + ar < m && t + ak < k) ? at_a(A, ib + ar, t + ak) : real_t(0);
+		    const int bk = e % GEMM_TILE;
+		    const int bc = e / GEMM_TILE;
+		    bs[bc][bk] = (t + bk < k && jb + bc < n) ? at_b(B, t + bk, jb + bc) : real_t(0);
+		}
 		__syncthreads();
-		for (int l = 0; l < GEMM_TILE; l++)
-		    sum += as[l][threadIdx.x] * bs[threadIdx.y][l];
+		for (int l = 0; l < GEMM_TILE; l++) {
+		    real_t ra[GEMM_THREAD_TILE];
+		    real_t rb[GEMM_THREAD_TILE];
+		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
+			ra[u] = as[l][threadIdx.x + u * GEMM_TILE];
+		    for (int v = 0; v < GEMM_THREAD_TILE; v++)
+			rb[v] = bs[threadIdx.y + v * GEMM_TILE][l];
+		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
+			for (int v = 0; v < GEMM_THREAD_TILE; v++)
+			    acc[u][v] += ra[u] * rb[v];
+		}
 		__syncthreads();
 	    }
-	    if (i < m && j < n)
-		cp[i + j * C.ld] = alpha * sum + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
+	    for (int u = 0; u < GEMM_THREAD_TILE; u++) {
+		const long i = ib + threadIdx.x + u * GEMM_TILE;
+		for (int v = 0; v < GEMM_THREAD_TILE; v++) {
+		    const long j = jb + threadIdx.y + v * GEMM_TILE;
+		    if (i < m && j < n)
+			cp[i + j * C.ld] = alpha * acc[u][v] + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
+		}
+	    }
 	}
     }
 }
@@ -414,8 +437,8 @@ void gemm(const char* transa, const char* transb, real_t alpha, const Matrix& A,
     if (m <= 0 || n <= 0 || (alpha == real_t(0) && beta == real_t(1))) return;

     const dim3 block(GEMM_TILE, GEMM_TILE);
-    long gx = (m + block.x - 1) / block.x;
-    long gy = (n + block.y - 1) / block.y;
+    long gx = (m + GEMM_BLOCK_TILE - 1) / GEMM_BLOCK_TILE;
+    long gy = (n + GEMM_BLOCK_TILE - 1) / GEMM_BLOCK_TILE;
     //temporary cap on grid size
     if (gx > 65535) gx = 65535;
     if (gy > 65535) gy = 65535;

Check the traffic reduction produced by register blocking

This was the biggest GEMM improvement in the sequence. Each warp performs more arithmetic with the values it already has, and independent accumulators give the compiler more work to schedule. Occupancy decreased, but useful throughput increased.

GEMM 3. 8 × 8 register tile

The next change doubles both register-tile dimensions. Each thread owns 64 outputs, and the block covers \(128\times128\). Sixteen shared operands now feed 64 FMAs, giving a scalar operand/FMA ratio of \(0.25\).

Changing one constant expands the thread tile to 8 × 8 and the block tile to 128 × 128. The indexing, shared arrays, and launch geometry all follow that constant.

@@ -7,7 +7,7 @@
 #define GEMV_ROWS 32
 #define GEMV_SPLITS 16
 #define GEMM_TILE 16
-#define GEMM_THREAD_TILE 4
+#define GEMM_THREAD_TILE 8
 #define GEMM_BLOCK_TILE (GEMM_TILE * GEMM_THREAD_TILE)

 static long span_of(int n, int inc) {

Compare the instruction and traffic effects of the larger tile

The compiled kernel uses about 170 registers per thread and supports only one of these blocks per multiprocessor in the recorded configuration. Theoretical occupancy falls to about 16.7%. The scheduler statistics show no eligible warp in roughly 45.65% of cycles.

The larger tile saves work, but all resident warps now belong to one block. When that block waits at a barrier or on dependencies, there is no second resident block to draw from. This is a plausible explanation for the limited gain, supported by the resource and scheduler measurements; the scalar reuse equation alone does not predict the result.

GEMM 4. Double buffering

I allocated two shared-memory buffers for each input tile. The kernel preloads the first tile, then loads the next tile into registers before computing on the current shared tile. After the current work, it writes those prefetched values to the alternate buffer, swaps buffers, and synchronizes. This reduces the loop from two block barriers per reduction stage to one.

The double-buffered shared-memory footprint is

This implementation uses ordinary register prefetching, not asynchronous copy instructions. The compiler and hardware determine how much load latency overlaps with independent arithmetic.

Two shared buffers alternate between tiles. The kernel primes buffer zero, prefetches the next tile into pa/pb registers, and stages it into the alternate buffer using ordinary loads.

		if (tn < k) {
		    for (int q = 0; q < GEMM_THREAD_TILE; q++) {
			const int e = tid + q * GEMM_TILE * GEMM_TILE;
			as[buf ^ 1][e / GEMM_BLOCK_TILE][e % GEMM_BLOCK_TILE] = pa[q];
			bs[buf ^ 1][e / GEMM_TILE][e % GEMM_TILE] = pb[q];
		    }
		    buf ^= 1;
		}
		__syncthreads();

@@ -376,35 +376,58 @@ void symv(const char* uplo, real_t alpha, const Matrix& A, const Vector& x, real

 template <typename AccessA, typename AccessB>
 __global__ void gemm_kernel(real_t alpha, const Matrix A, const Matrix B, real_t beta, Matrix C, long m, long n, long k, AccessA at_a, AccessB at_b) {
-    __shared__ real_t as[GEMM_TILE][GEMM_BLOCK_TILE];
-    __shared__ real_t bs[GEMM_BLOCK_TILE][GEMM_TILE];
+    __shared__ real_t as[2][GEMM_TILE][GEMM_BLOCK_TILE];
+    __shared__ real_t bs[2][GEMM_BLOCK_TILE][GEMM_TILE];
     real_t* __restrict__ cp = C.data;
     const int tid = threadIdx.y * GEMM_TILE + threadIdx.x;
     for (long ib = blockIdx.x * (long)GEMM_BLOCK_TILE; ib < m; ib += (long)gridDim.x * GEMM_BLOCK_TILE) {
 	for (long jb = blockIdx.y * (long)GEMM_BLOCK_TILE; jb < n; jb += (long)gridDim.y * GEMM_BLOCK_TILE) {
 	    real_t acc[GEMM_THREAD_TILE][GEMM_THREAD_TILE] = {};
+	    for (int q = 0; q < GEMM_THREAD_TILE; q++) {
+		const int e = tid + q * GEMM_TILE * GEMM_TILE;
+		const int ar = e % GEMM_BLOCK_TILE;
+		const int ak = e / GEMM_BLOCK_TILE;
+		as[0][ak][ar] = (ib + ar < m && ak < k) ? at_a(A, ib + ar, ak) : real_t(0);
+		const int bk = e % GEMM_TILE;
+		const int bc = e / GEMM_TILE;
+		bs[0][bc][bk] = (bk < k && jb + bc < n) ? at_b(B, bk, jb + bc) : real_t(0);
+	    }
+	    __syncthreads();
+	    int buf = 0;
 	    for (long t = 0; t < k; t += GEMM_TILE) {
-		for (int q = 0; q < GEMM_THREAD_TILE; q++) {
-		    const int e = tid + q * GEMM_TILE * GEMM_TILE;
-		    const int ar = e % GEMM_BLOCK_TILE;
-		    const int ak = e / GEMM_BLOCK_TILE;
-		    as[ak][ar] = (ib + ar < m && t + ak < k) ? at_a(A, ib + ar, t + ak) : real_t(0);
-		    const int bk = e % GEMM_TILE;
-		    const int bc = e / GEMM_TILE;
-		    bs[bc][bk] = (t + bk < k && jb + bc < n) ? at_b(B, t + bk, jb + bc) : real_t(0);
+		const long tn = t + GEMM_TILE;
+		real_t pa[GEMM_THREAD_TILE];
+		real_t pb[GEMM_THREAD_TILE];
+		if (tn < k) {
+		    for (int q = 0; q < GEMM_THREAD_TILE; q++) {
+			const int e = tid + q * GEMM_TILE * GEMM_TILE;
+			const int ar = e % GEMM_BLOCK_TILE;
+			const int ak = e / GEMM_BLOCK_TILE;
+			pa[q] = (ib + ar < m && tn + ak < k) ? at_a(A, ib + ar, tn + ak) : real_t(0);
+			const int bk = e % GEMM_TILE;
+			const int bc = e / GEMM_TILE;
+			pb[q] = (tn + bk < k && jb + bc < n) ? at_b(B, tn + bk, jb + bc) : real_t(0);
+		    }
 		}
-		__syncthreads();
 		for (int l = 0; l < GEMM_TILE; l++) {
 		    real_t ra[GEMM_THREAD_TILE];
 		    real_t rb[GEMM_THREAD_TILE];
 		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
-			ra[u] = as[l][threadIdx.x + u * GEMM_TILE];
+			ra[u] = as[buf][l][threadIdx.x + u * GEMM_TILE];
 		    for (int v = 0; v < GEMM_THREAD_TILE; v++)
-			rb[v] = bs[threadIdx.y + v * GEMM_TILE][l];
+			rb[v] = bs[buf][threadIdx.y + v * GEMM_TILE][l];
 		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
 			for (int v = 0; v < GEMM_THREAD_TILE; v++)
 			    acc[u][v] += ra[u] * rb[v];
 		}
+		if (tn < k) {
+		    for (int q = 0; q < GEMM_THREAD_TILE; q++) {
+			const int e = tid + q * GEMM_TILE * GEMM_TILE;
+			as[buf ^ 1][e / GEMM_BLOCK_TILE][e % GEMM_BLOCK_TILE] = pa[q];
+			bs[buf ^ 1][e / GEMM_TILE][e % GEMM_TILE] = pb[q];
+		    }
+		    buf ^= 1;
+		}
 		__syncthreads();
 	    }
 	    for (int u = 0; u < GEMM_THREAD_TILE; u++) {

Check whether double buffering reduced the amount of work

The scheduler measurements make the difference clearer: issue-slot utilization rises from about 52.7% to 60.4%, and cycles with no eligible warp fall from 45.65% to 37.67%. Registers rise from 170 to 176 per thread, with no recorded spills. The remaining idle cycles explain why prefetching improves performance without closing the entire cuBLAS gap.

Complete GEMM implementation

This includes the complete double-buffered kernel, transpose dispatch, and host entry point. This listing is the current working version, not a byte-for-byte copy of c11b4fd: the later GEMM_K_TILE and GEMM_LOADS refactor separates reduction depth from launch geometry and permits the smaller FP64 shared-memory footprint. The historical diffs above retain the original constants; the screenshots are not measurements of that later refactor. With the FP32 definitions below, each block has 256 threads and computes a 128 × 128 output tile. The shared types, constants, and memory helpers are listed at the end.

template <typename AccessA, typename AccessB>
__global__ void gemm_kernel(real_t alpha, const Matrix A, const Matrix B, real_t beta, Matrix C, long m, long n, long k, AccessA at_a, AccessB at_b) {
    __shared__ real_t as[2][GEMM_K_TILE][GEMM_BLOCK_TILE];
    __shared__ real_t bs[2][GEMM_BLOCK_TILE][GEMM_K_TILE];
    real_t* __restrict__ cp = C.data;
    const int tid = threadIdx.y * GEMM_TILE + threadIdx.x;
    for (long ib = blockIdx.x * (long)GEMM_BLOCK_TILE; ib < m; ib += (long)gridDim.x * GEMM_BLOCK_TILE) {
	for (long jb = blockIdx.y * (long)GEMM_BLOCK_TILE; jb < n; jb += (long)gridDim.y * GEMM_BLOCK_TILE) {
	    real_t acc[GEMM_THREAD_TILE][GEMM_THREAD_TILE] = {};
	    for (int q = 0; q < GEMM_LOADS; q++) {
		const int e = tid + q * GEMM_TILE * GEMM_TILE;
		const int ar = e % GEMM_BLOCK_TILE;
		const int ak = e / GEMM_BLOCK_TILE;
		as[0][ak][ar] = (ib + ar < m && ak < k) ? at_a(A, ib + ar, ak) : real_t(0);
		const int bk = e % GEMM_K_TILE;
		const int bc = e / GEMM_K_TILE;
		bs[0][bc][bk] = (bk < k && jb + bc < n) ? at_b(B, bk, jb + bc) : real_t(0);
	    }
	    __syncthreads();
	    int buf = 0;
	    for (long t = 0; t < k; t += GEMM_K_TILE) {
		const long tn = t + GEMM_K_TILE;
		real_t pa[GEMM_LOADS];
		real_t pb[GEMM_LOADS];
		if (tn < k) {
		    for (int q = 0; q < GEMM_LOADS; q++) {
			const int e = tid + q * GEMM_TILE * GEMM_TILE;
			const int ar = e % GEMM_BLOCK_TILE;
			const int ak = e / GEMM_BLOCK_TILE;
			pa[q] = (ib + ar < m && tn + ak < k) ? at_a(A, ib + ar, tn + ak) : real_t(0);
			const int bk = e % GEMM_K_TILE;
			const int bc = e / GEMM_K_TILE;
			pb[q] = (tn + bk < k && jb + bc < n) ? at_b(B, tn + bk, jb + bc) : real_t(0);
		    }
		}
		for (int l = 0; l < GEMM_K_TILE; l++) {
		    real_t ra[GEMM_THREAD_TILE];
		    real_t rb[GEMM_THREAD_TILE];
		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
			ra[u] = as[buf][l][threadIdx.x + u * GEMM_TILE];
		    for (int v = 0; v < GEMM_THREAD_TILE; v++)
			rb[v] = bs[buf][threadIdx.y + v * GEMM_TILE][l];
		    for (int u = 0; u < GEMM_THREAD_TILE; u++)
			for (int v = 0; v < GEMM_THREAD_TILE; v++)
			    acc[u][v] += ra[u] * rb[v];
		}
		if (tn < k) {
		    for (int q = 0; q < GEMM_LOADS; q++) {
			const int e = tid + q * GEMM_TILE * GEMM_TILE;
			as[buf ^ 1][e / GEMM_BLOCK_TILE][e % GEMM_BLOCK_TILE] = pa[q];
			bs[buf ^ 1][e / GEMM_K_TILE][e % GEMM_K_TILE] = pb[q];
		    }
		    buf ^= 1;
		}
		__syncthreads();
	    }
	    for (int u = 0; u < GEMM_THREAD_TILE; u++) {
		const long i = ib + threadIdx.x + u * GEMM_TILE;
		for (int v = 0; v < GEMM_THREAD_TILE; v++) {
		    const long j = jb + threadIdx.y + v * GEMM_TILE;
		    if (i < m && j < n)
			cp[i + j * C.ld] = alpha * acc[u][v] + (beta == real_t(0) ? real_t(0) : beta * cp[i + j * C.ld]);
		}
	    }
	}
    }
}

template <typename AccessA>
static void gemm_launch(const char* transb, real_t alpha, const Matrix& A, const Matrix& B, real_t beta, Matrix& C, long m, long n, long k, dim3 grid, dim3 block, AccessA at_a) {
    if (transb[0] == 'N' || transb[0] == 'n') {
	gemm_kernel<<<grid, block>>>(alpha, A, B, beta, C, m, n, k, at_a, NoTransAt());
    } else {
	gemm_kernel<<<grid, block>>>(alpha, A, B, beta, C, m, n, k, at_a, TransAt());
    }
}

void gemm(const char* transa, const char* transb, real_t alpha, const Matrix& A, const Matrix& B, real_t beta, Matrix& C) {
    const bool na = (transa[0] == 'N' || transa[0] == 'n');
    const long m = na ? A.rows : A.cols;
    const long k = na ? A.cols : A.rows;
    const long n = (transb[0] == 'N' || transb[0] == 'n') ? B.cols : B.rows;

    if (m <= 0 || n <= 0 || (alpha == real_t(0) && beta == real_t(1))) return;

    const dim3 block(GEMM_TILE, GEMM_TILE);
    long gx = (m + GEMM_BLOCK_TILE - 1) / GEMM_BLOCK_TILE;
    long gy = (n + GEMM_BLOCK_TILE - 1) / GEMM_BLOCK_TILE;
    //temporary cap on grid size
    if (gx > 65535) gx = 65535;
    if (gy > 65535) gy = 65535;
    const dim3 grid((unsigned)gx, (unsigned)gy);

    real_t* sa = nullptr;
    real_t* sb = nullptr;
    real_t* sc = nullptr;
    const Matrix dA = stage_matrix(A, sa);
    const Matrix dB = stage_matrix(B, sb);
    Matrix dC = stage_matrix(C, sc);

    if (na) {
	gemm_launch(transb, alpha, dA, dB, beta, dC, m, n, k, grid, block, NoTransAt());
    } else {
	gemm_launch(transb, alpha, dA, dB, beta, dC, m, n, k, grid, block, TransAt());
    }
    CUDA_ERROR_CHECK(cudaGetLastError());
    unstage_matrix(C, sc);
    release(sb);
    release(sa);
}

Scaling with input size

The captures above measure one size per routine. To see whether the conclusions hold elsewhere, I rebuilt the saved kernel implementations and timed them across a range of sizes on the same RTX4060 alongside cuBLAS. Each implementation uses its original grid and block configuration. Timing uses CUDA events around launches on device-resident data, with the L2 cache flushed before every repetition. At each size all implementations run round-robin, so clock changes affect them alike, and each point is the median of three passes of at least 20 repetitions. Every result was checked against cuBLAS.

These numbers are not interchangeable with the Nsight Compute captures. Nsight Compute locks the GPU to its base clocks, while these runs let it boost, so the compute-bound GEMM kernels run faster here: 4.46 ms instead of 4.91 ms for the final version at 2048. Under sustained load the GPU also reached its 40 W power limit and sometimes lowered the memory clock, which pushes the medians of the bandwidth-bound DOT and GEMV kernels 12 to 15% above the captures, although their fastest repetitions agree with the captures to within about 1%. Compare curves within a chart rather than reading single points against the earlier tables. The raw measurements are available as a CSV file.

For DOT, the final kernel matches cuBLAS at every length, from 7.2 versus 10.2 µs at 2¹⁶ elements to 2,456 versus 2,464 µs at 2²⁶. Below about a million elements both are limited by launch and reduction overhead rather than by bandwidth, which is why every curve starts low. The chart also shows that step DOT 5 was not specific to 2²⁴ elements: at every length above a million, the two versions that give each thread one element stay 20 to 40% below the ones that give it two. The atomic baseline does not improve with size at all: its time grows in proportion to n, so its bandwidth stays flat. Its FP32 result also drifts as n grows, because millions of single additions into one accumulator lose low-order bits; at 2²⁶ elements its relative error reached 22%, while the reduction trees stay within 4 × 10⁻⁶ of cuBLAS.

For GEMV, the row and column split and cuBLAS are indistinguishable across the whole range. The baseline’s gap depends strongly on size. Its eight-block launch was written for the 4096 problem; on smaller matrices it has even less independent work per SM, and on larger ones each thread’s longer row loop partly compensates.

GEMM shows the clearest size dependence. From order 1024 upward the final kernel stays about 1.6 to 1.7 times behind cuBLAS, a stable gap that matches the single capture. The ordering of the versions changes at the small end. At order 256 the double-buffered kernel launches only four 128 × 128 blocks and takes 55.3 µs, while the 4 × 4 version, with sixteen 64 × 64 blocks, takes 27.6 µs and cuBLAS takes 15.4 µs. The tile size that maximizes reuse on a large matrix starves the GPU of blocks on a small one, which is why libraries such as cuBLAS choose among several kernels by problem size.

What I learned about acceleration

The useful question changed at each stage. For DOT, it was how to combine results without serializing every thread. For GEMV, it was how to expose enough independent work while respecting column-major memory layout. For GEMM, it was how many times a value could contribute to arithmetic after each load. The same GPU needed three different strategies because the operations place different demands on it.

DOT made the cost of coordination especially clear. Millions of independent multiplications still ran slowly when every thread updated the same scalar. Reducing within a block removed most of that contention; giving each thread more work then reduced how often the kernel had to combine partial sums. Changing the reduction tree also changed floating-point rounding, so checking the answer remained part of evaluating each improvement.

GEMV showed that reducing memory accesses is only part of the problem. Caching the input vector helped, but an eight-block launch still left too little independent work for the GPU. Splitting columns exposed more parallelism while keeping neighboring rows contiguous in memory. The lesson was to consider data layout and work distribution together: a kernel can access the right data efficiently and still leave much of the machine idle.

GEMM made reuse tangible. Shared tiles let threads cooperate on input loads, and register tiles let each thread use those values to accumulate several outputs. But larger tiles also reduce the number of blocks. The size sweep showed the consequence: a configuration that performed well on a large matrix could lose on a small one because it no longer supplied enough parallel work. Tile size is a tradeoff between reuse, resource use, and the amount of work available.

The roofline model helped distinguish operations with little reuse from operations with considerable reuse potential. The profiles explained why a particular implementation fell short. A traffic equation could predict the effect of tiling on requested bytes, but it could not by itself account for atomic contention, shared-memory instruction pressure, or a block waiting at a barrier. A useful optimization process starts with a prediction about the bottleneck, tests a code change, and checks both runtime and the profile against that prediction.

Acceleration comes from coordinating the mathematical structure with the machine: combining locally, mapping contiguous data to neighboring lanes, accumulating many outputs in registers, and arranging the remaining loads around independent work. Once one bottleneck is removed, another can become the limiting factor, so each improvement changes what is worth investigating next.

These improvements are kernel-level results. The current host-pointer API also allocates memory, transfers inputs, and copies outputs back. Those costs must be measured separately when deciding whether an application benefits. A useful next step is to keep data on the GPU across several operations, so transfers and allocations can be spread over more computation. A faster kernel matters most when the surrounding application can make use of it.

Shared code and a complete source file

The functions above use the following shared definitions and staging helpers. These are included once to keep each kernel listing focused. The code is a snapshot of the working implementation, including its host-memory API, rather than a claim of complete BLAS conformance.

Download the complete CUDA source. It combines this prelude with all three implementations and can be compiled as a translation unit:

nvcc -O3 -arch=sm_89 -c cuda-blas-kernels.cu

#include <cuda_runtime.h>
#include <iostream>
#include <cstdlib>

#ifdef DOUBLE_PRECISION
    typedef double real_t;
#else
    typedef float real_t;
#endif

//Vector: a strided view of n elements at data[0], data[inc], data[2*inc], ...
//A matrix column is {data + j*ld, rows, 1}; a row is {data + i, cols, ld}.
struct Vector {
    real_t* data;
    int n;
    int inc;  // stride between elements, in units of real_t (1 = contiguous)
};

//Matrix
struct Matrix {
    real_t* data;
    int rows, cols, ld;
};


struct NoTransAt {
    __device__ real_t operator()(const Matrix& A, long i, long j) const {
        return A.data[i + j * A.ld];
    }
};

struct TransAt {
    __device__ real_t operator()(const Matrix& A, long i, long j) const {
        return A.data[j + i * A.ld];
    }
};

#define CUDA_ERROR_CHECK(call)                                                             \
    do {                                                                                  \
        cudaError_t status = call;                                                        \
        if (status != cudaSuccess) {                                                      \
            std::cerr << "CUDA error: " << cudaGetErrorString(status) << std::endl;       \
            std::exit(1);                                                                 \
        }                                                                                 \
    } while (0)

#define BLOCK_SIZE 256
#define GEMV_BLOCK_SIZE BLOCK_SIZE
#define GEMV_ROWS 32
#define GEMV_SPLITS 16
#ifndef GEMM_TILE
#define GEMM_TILE 16
#endif
#ifndef GEMM_K_TILE
// Keep static shared storage below the per-block 48 KiB limit in FP64.
#ifdef DOUBLE_PRECISION
#define GEMM_K_TILE 8
#else
#define GEMM_K_TILE 16
#endif
#endif
#define GEMM_LOADS (GEMM_K_TILE * GEMM_BLOCK_TILE / (GEMM_TILE * GEMM_TILE))
#define GEMM_THREAD_TILE 8
#define GEMM_BLOCK_TILE (GEMM_TILE * GEMM_THREAD_TILE)
static_assert(GEMM_K_TILE * GEMM_BLOCK_TILE % (GEMM_TILE * GEMM_TILE) == 0,
              "GEMM tile loads must divide evenly across the block");

static long span_of(int n, int inc) {
    return n > 0 ? 1 + (long)(n - 1) * (inc < 0 ? -inc : inc) : 0;
}

static long offset_of(int n, int inc) {
    return (inc < 0 && n > 0) ? (long)(n - 1) * -inc : 0;
}

static real_t* device_copy_in(const real_t* host, long count) {
    real_t* dev = nullptr;
    CUDA_ERROR_CHECK(cudaMalloc(&dev, (count > 0 ? count : 1) * sizeof(real_t)));
    if (count > 0)
	CUDA_ERROR_CHECK(cudaMemcpy(dev, host, count * sizeof(real_t), cudaMemcpyHostToDevice));
    return dev;
}

static Vector stage_vector(const Vector& v, real_t*& slab) {
    const long off = offset_of(v.n, v.inc);
    slab = device_copy_in(v.data - off, span_of(v.n, v.inc));
    return Vector{slab + off, v.n, v.inc};
}

static void unstage_vector(const Vector& v, real_t* slab) {
    const long count = span_of(v.n, v.inc);
    if (count > 0)
	CUDA_ERROR_CHECK(cudaMemcpy(v.data - offset_of(v.n, v.inc), slab, count * sizeof(real_t), cudaMemcpyDeviceToHost));
    CUDA_ERROR_CHECK(cudaFree(slab));
}

static Matrix stage_matrix(const Matrix& A, real_t*& slab) {
    slab = device_copy_in(A.data, (long)A.ld * A.cols);
    return Matrix{slab, A.rows, A.cols, A.ld};
}

static void unstage_matrix(const Matrix& A, real_t* slab) {
    const long count = (long)A.ld * A.cols;
    if (count > 0)
	CUDA_ERROR_CHECK(cudaMemcpy(A.data, slab, count * sizeof(real_t), cudaMemcpyDeviceToHost));
    CUDA_ERROR_CHECK(cudaFree(slab));
}

static void release(real_t* slab) {
    CUDA_ERROR_CHECK(cudaFree(slab));
}