under test, not sure no errors

upstream_ref/sgemm_edtallison/01_naive.cuh

@@ -0,0 +1,37 @@
+# pragma once
+
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+/*
+
+Matrix sizes:
+MxK * KxN = MxN
+
+*/
+
+__global__ void sgemm_naive(
+    int M, int N, int K, // sizes
+    float alpha, const float *A, const float *B, float beta, float *C // pointers used to point to matrices
+) {
+    // compute position in C that this thread is responsible for
+    // "which block" * "width of block" to get to start of block + "which thread"
+    const uint x = blockIdx.x * blockDim.x + threadIdx.x; // "which row?" (inverted from graphical intuition, confusingly)
+    const uint y = blockIdx.y * blockDim.y + threadIdx.y; // "which column?"
+
+    // if M or N are not multiples of 32, there will be "extra"/"remainder" threads on the last block in x/y.
+    // we don't want those leftover threads to do anything (tile quantisation)
+    if (x < M && y < N) {
+        float tmp = 0.0;
+        for (int i = 0; i < K; ++i) { // K is the size of the row in A, col in B i.e. the dot product
+            // A: x * K gives the start of relevant row, i enumerates across the row (col by col)
+            // B: y gives the relevant column, i * N enumerates down the column, (row by row)
+            tmp += A[x * K + i] * B[i * N + y];
+        }
+        // C = alpha*(A@B) + beta*C
+        // x * N takes to start of relevant row, y moves across to the relevant column
+        C[x * N + y] = alpha * tmp + beta * C[x * N + y];
+    }
+}

upstream_ref/sgemm_edtallison/02_kernel_global_mem_coalesce.cuh

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+template <const uint BLOCKSIZE>
+// __global__ is used to specify that the function is run on GPU, called by host (CPU)
+__global__ void sgemm_global_mem_coalesce(int M, int N, int K, float alpha,
+                                          const float *A, const float *B,
+                                          float beta, float *C) {
+    const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); // note that blockDim is now 1-dimensional
+    const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);
+
+    // if statement is necessary to make things work under tile quantization
+    if (cRow < M && cCol < N) {
+        float tmp = 0.0;
+        for (int i = 0; i < K; ++i) {
+            tmp += A[cRow * K + i] * B[i * N + cCol];
+        }
+        C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol];
+    }
+}

upstream_ref/sgemm_edtallison/03_kernel_shared_mem_blocking.cuh

@@ -0,0 +1,54 @@
+#pragma once
+
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BLOCKSIZE>
+__global__ void sgemm_shared_mem_block(int M, int N, int K, float alpha,
+                                       const float *A, const float *B,
+                                       float beta, float *C) {
+    // output C block we want to compute with this threadBlock
+    const uint cRow = blockIdx.x;
+    const uint cCol = blockIdx.y;
+
+    // allocate buffer for current block in fast SMEM (shared between all threads in block)
+    __shared__ float As[BLOCKSIZE * BLOCKSIZE];
+    __shared__ float Bs[BLOCKSIZE * BLOCKSIZE];
+
+    // the inner row and col that we are accessing in this specific thread
+    const uint threadRow = threadIdx.x / BLOCKSIZE; // note similarity to previous kernel
+    const uint threadCol = threadIdx.x % BLOCKSIZE;
+
+    // advance pointers to the starting positions (they are input as pointers to first elements in the matrices)
+    A += cRow * BLOCKSIZE * K;                    // row=cRow, col=0 (the start of the relevant row)
+    B += cCol * BLOCKSIZE;                        // row=0, col=cCol (top of relevant col)
+    C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE; // row=cRow, col=cCol
+
+    float tmp = 0.0;
+    for (int bkIdx=0; bkIdx < K; bkIdx+=BLOCKSIZE) { // shifting the whole block along the row of A and col of B
+        // have each thread load one of the elements in A and B
+        // make the threadCol (=threadIdx.x) the consecutive index
+        // to allow GMEM access coalescing
+        As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
+        Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
+
+        // ensure cache is fully populated
+        __syncthreads();
+        A += BLOCKSIZE; // for next iteration
+        B += BLOCKSIZE * N;
+
+        // execute the dotproduct on the currently cached block
+        for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx) {
+            tmp += As[threadRow * BLOCKSIZE + dotIdx] *
+                   Bs[dotIdx * BLOCKSIZE + threadCol];
+        }
+        // sync so faster threads don't fetch the next block into cache
+        _syncthreads();
+    }
+    C[threadRow * N + threadCol] = alpha * tmp + beta * C[threadRow * N + threadCol];
+}

upstream_ref/sgemm_edtallison/04_kernel_1D_blocktiling.cuh

@@ -0,0 +1,80 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM>
+__global__ void sgemm1DBlocktiling(int M, int N, int K, float alpha,
+                                   const float *A, const float *B, float beta,
+                                   float *C) {
+  // If we flip x and y here we get ~30% less performance for large matrices.
+  // The current, 30% faster configuration ensures that blocks with sequential
+  // blockIDs access columns of B sequentially, while sharing the same row of A.
+  // The slower configuration would share columns of A, but access into B would
+  // be non-sequential. So the faster configuration has better spatial locality
+  // and hence a greater L2 hit rate.
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // each warp will calculate 32*TM elements, with 32 being the columnar dim.
+  const int threadCol = threadIdx.x % BN;
+  const int threadRow = threadIdx.x / BN;
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // todo: adjust this to each thread to load multiple entries and
+  // better exploit the cache sizes
+  assert(BM * BK == blockDim.x);
+  assert(BN * BK == blockDim.x);
+  const uint innerColA = threadIdx.x % BK; // warp-level GMEM coalescing
+  const uint innerRowA = threadIdx.x / BK;
+  const uint innerColB = threadIdx.x % BN; // warp-level GMEM coalescing
+  const uint innerRowB = threadIdx.x / BN;
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM] = {0.0};
+
+  // outer loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
+    Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;
+    B += BK * N;
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // we make the dotproduct loop the outside loop, which facilitates
+      // reuse of the Bs entry, which we can cache in a tmp var.
+      float tmpB = Bs[dotIdx * BN + threadCol];
+      for (uint resIdx = 0; resIdx < TM; ++resIdx) {
+        threadResults[resIdx] +=
+            As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdx = 0; resIdx < TM; ++resIdx) {
+    C[(threadRow * TM + resIdx) * N + threadCol] =
+        alpha * threadResults[resIdx] +
+        beta * C[(threadRow * TM + resIdx) * N + threadCol];
+  }
+}

upstream_ref/sgemm_edtallison/05_kernel_2D_blocktiling.cuh

@@ -0,0 +1,102 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void __launch_bounds__((BM * BN) / (TM * TN), 1)
+    sgemm2DBlocktiling(int M, int N, int K, float alpha, const float *A,
+                       const float *B, float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  const uint totalResultsBlocktile = BM * BN;
+  // A thread is responsible for calculating TM*TN elements in the blocktile
+  const uint numThreadsBlocktile = totalResultsBlocktile / (TM * TN);
+
+  // ResultsPerBlock / ResultsPerThread == ThreadsPerBlock
+  assert(numThreadsBlocktile == blockDim.x);
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  const uint innerRowA = threadIdx.x / BK;
+  const uint innerColA = threadIdx.x % BK;
+  // calculates the number of rows of As that are being loaded in a single step
+  // by a single block
+  const uint strideA = numThreadsBlocktile / BK;
+  const uint innerRowB = threadIdx.x / BN;
+  const uint innerColB = threadIdx.x % BN;
+  // for both As and Bs we want each load to span the full column-width, for
+  // better GMEM coalescing (as opposed to spanning full row-width and iterating
+  // across columns)
+  const uint strideB = numThreadsBlocktile / BN;
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  // register caches for As and Bs
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
+      As[(innerRowA + loadOffset) * BK + innerColA] =
+          A[(innerRowA + loadOffset) * K + innerColA];
+    }
+    for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
+      Bs[(innerRowB + loadOffset) * BN + innerColB] =
+          B[(innerRowB + loadOffset) * N + innerColB];
+    }
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+    for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+      C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN] =
+          alpha * threadResults[resIdxM * TN + resIdxN] +
+          beta * C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN];
+    }
+  }
+}

upstream_ref/sgemm_edtallison/06_kernel_vectorize.cuh

@@ -0,0 +1,98 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmVectorize(int M, int N, int K, float alpha, float *A,
+                               float *B, float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
+        reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/07_kernel_resolve_bank_conflicts.cuh

@@ -0,0 +1,103 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmResolveBankConflicts(int M, int N, int K, float alpha,
+                                          float *A, float *B, float beta,
+                                          float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    // "linearize" Bs while storing it
+    tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 0) * 16 + innerColB / 2] = tmp.x;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 1) * 16 + innerColB / 2] = tmp.y;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 2) * 16 + innerColB / 2] = tmp.z;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 3) * 16 + innerColB / 2] = tmp.w;
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[(dotIdx * 8 + i) * 16 + threadCol];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/08_kernel_bank_extra_col.cuh

@@ -0,0 +1,103 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmResolveBankExtraCol(int M, int N, int K, float alpha,
+                                         float *A, float *B, float beta,
+                                         float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  const int extraCols = 5;
+  __shared__ float Bs[BK * (BN + extraCols)];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 0] = tmp.x;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 1] = tmp.y;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 2] = tmp.z;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 3] = tmp.w;
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * (BN + extraCols) + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/09_kernel_autotuned.cuh

@@ -0,0 +1,127 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+const int K9_NUM_THREADS = 256;
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void __launch_bounds__(K9_NUM_THREADS)
+    sgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
+                   float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // size of warptile
+  constexpr int WM = TM * 16;
+  constexpr int WN = TN * 16;
+  // iterations of warptile
+  constexpr int WMITER = CEIL_DIV(BM, WM);
+  constexpr int WNITER = CEIL_DIV(BN, WN);
+
+  // Placement of the thread in the warptile
+  const int threadCol = threadIdx.x % (WN / TN);
+  const int threadRow = threadIdx.x / (WN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (K9_NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = K9_NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * WNITER * TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+      float4 tmp = reinterpret_cast<float4 *>(
+          &A[(innerRowA + offset) * K + innerColA * 4])[0];
+      // transpose A while storing it
+      As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+      As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+      As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+      As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+    }
+
+    for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+      reinterpret_cast<float4 *>(
+          &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+          reinterpret_cast<float4 *>(
+              &B[(innerRowB + offset) * N + innerColB * 4])[0];
+    }
+    __syncthreads();
+
+    for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
+      for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
+        // calculate per-thread results
+        for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+          // block into registers
+          for (uint i = 0; i < TM; ++i) {
+            regM[i] = As[dotIdx * BM + (wmIdx * WM) + threadRow * TM + i];
+          }
+          for (uint i = 0; i < TN; ++i) {
+            regN[i] = Bs[dotIdx * BN + (wnIdx * WN) + threadCol * TN + i];
+          }
+          for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+            for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+              threadResults[(wmIdx * TM + resIdxM) * (WNITER * TN) +
+                            wnIdx * TN + resIdxN] +=
+                  regM[resIdxM] * regN[resIdxN];
+            }
+          }
+        }
+      }
+    }
+    __syncthreads();
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+  }
+
+  // write out the results
+  for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
+    for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
+      float *C_interim = C + (wmIdx * WM * N) + (wnIdx * WN);
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRow * TM + resIdxM) * N + threadCol * TN +
+                         resIdxN])[0];
+          // perform GEMM update in reg
+          const int i =
+              (wmIdx * TM + resIdxM) * (WNITER * TN) + wnIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(&C_interim[(threadRow * TM + resIdxM) * N +
+                                                threadCol * TN + resIdxN])[0] =
+              tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/10_kernel_warptiling.cuh

@@ -0,0 +1,187 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+const int WARPSIZE = 32; // warpSize is not constexpr
+
+namespace wt {
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB>
+__device__ void loadFromGmem(int N, int K, const float *A, const float *B,
+                             float *As, float *Bs, int innerRowA, int innerColA,
+                             int innerRowB, int innerColB) {
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    const float4 tmp = reinterpret_cast<const float4 *>(
+        &A[(innerRowA + offset) * K + innerColA * 4])[0];
+    // float4 tmp;
+    // asm("ld.global.nc.v4.f32 {%0, %1, %2, %3}, [%4];"
+    //     : "=f"(tmp.x), "=f"(tmp.y), "=f"(tmp.z), "=f"(tmp.w)
+    //     : "l"(&A[(innerRowA + offset) * K + innerColA * 4]));
+    As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    reinterpret_cast<float4 *>(
+        &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+        reinterpret_cast<const float4 *>(
+            &B[(innerRowB + offset) * N + innerColB * 4])[0];
+    // asm("ld.global.v4.f32 {%0, %1, %2, %3}, [%4];"
+    //     : "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 0]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 1]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 2]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 3])
+    //     : "l"(&B[(innerRowB + offset) * N + innerColB * 4]));
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace wt
+
+/*
+ * @tparam BM The threadblock size for M dimension SMEM caching.
+ * @tparam BN The threadblock size for N dimension SMEM caching.
+ * @tparam BK The threadblock size for K dimension SMEM caching.
+ * @tparam WM M dim of continuous tile computed by each warp
+ * @tparam WN N dim of continuous tile computed by each warp
+ * @tparam WMITER The number of subwarp tiling steps in M dimension.
+ * @tparam WNITER The number of subwarp tiling steps in N dimension.
+ * @tparam TM The per-thread tile size for M dimension.
+ * @tparam TN The per-thread tile size for N dimension.
+ */
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    sgemmWarptiling(int M, int N, int K, float alpha, float *A, float *B,
+                    float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    wt::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB);
+    __syncthreads();
+    wt::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                        TN>(regM, regN, threadResults, As, Bs, warpRow, warpCol,
+                            threadRowInWarp, threadColInWarp);
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/11_kernel_double_buffering.cuh

@@ -0,0 +1,220 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+namespace db {
+
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB>
+__device__ void loadFromGmem(const int N, const int K, float *A, float *B,
+                             float *As, float *Bs, const int innerRowA,
+                             const int innerColA, const int innerRowB,
+                             const int innerColB) {
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    float4 tmp = reinterpret_cast<float4 *>(
+        &A[(innerRowA + offset) * K + innerColA * 4])[0];
+    // transpose A while storing it
+    As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    reinterpret_cast<float4 *>(
+        &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+        reinterpret_cast<float4 *>(
+            &B[(innerRowB + offset) * N + innerColB * 4])[0];
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace db
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    sgemmDoubleBuffering(const int M, const int N, const int K,
+                         const float alpha, float *A, float *B, float beta,
+                         float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[2 * BM * BK];
+  __shared__ float Bs[2 * BK * BN];
+
+  // setup double buffering split
+  bool doubleBufferIdx = threadIdx.x >= (NUM_THREADS / 2);
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // for the loading, we're pretending like there's half as many threads
+  // as there actually are
+  const uint innerRowA = (threadIdx.x % (NUM_THREADS / 2)) / (BK / 4);
+  const uint innerColA = (threadIdx.x % (NUM_THREADS / 2)) % (BK / 4);
+  constexpr uint rowStrideA = ((NUM_THREADS / 2) * 4) / BK;
+  const uint innerRowB = (threadIdx.x % (NUM_THREADS / 2)) / (BN / 4);
+  const uint innerColB = (threadIdx.x % (NUM_THREADS / 2)) % (BN / 4);
+  constexpr uint rowStrideB = (NUM_THREADS / 2) / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  if (doubleBufferIdx == 0) {
+    // load first (B0)
+    db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB);
+  }
+  __syncthreads();
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += 2 * BK) {
+    if (doubleBufferIdx == 0) {
+      // process current (B0)
+      db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                          TN>(regM, regN, threadResults, As, Bs, warpRow,
+                              warpCol, threadRowInWarp, threadColInWarp);
+      __syncthreads();
+
+      // process current+1 (B1)
+      if (bkIdx + BK < K) {
+        db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN,
+                            TM, TN>(regM, regN, threadResults, As + (BM * BK),
+                                    Bs + (BK * BN), warpRow, warpCol,
+                                    threadRowInWarp, threadColInWarp);
+      }
+      __syncthreads();
+
+      // load current + 2 (B0)
+      if (bkIdx + 2 * BK < K) {
+        db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+            N, K, A + 2 * BK, B + 2 * BK * N, As, Bs, innerRowA, innerColA,
+            innerRowB, innerColB);
+      }
+    } else {
+      // load current + 1 (B1)
+      if (bkIdx + BK < K) {
+        db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+            N, K, A + BK, B + BK * N, As + (BM * BK), Bs + (BK * BN), innerRowA,
+            innerColA, innerRowB, innerColB);
+      }
+      __syncthreads();
+
+      // process current (B0)
+      db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                          TN>(regM, regN, threadResults, As, Bs, warpRow,
+                              warpCol, threadRowInWarp, threadColInWarp);
+      __syncthreads();
+
+      // process current+1 (B1)
+      if (bkIdx + BK < K) {
+        db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN,
+                            TM, TN>(regM, regN, threadResults, As + (BM * BK),
+                                    Bs + (BK * BN), warpRow, warpCol,
+                                    threadRowInWarp, threadColInWarp);
+      }
+    }
+
+    A += 2 * BK;     // move BK columns to right
+    B += 2 * BK * N; // move BK rows down
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/12_kernel_double_buffering.cuh

@@ -0,0 +1,229 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cooperative_groups.h>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda/barrier>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+namespace {
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB, typename T>
+__device__ void loadFromGmem(int N, int K, float *A, float *B, float *As,
+                             float *Bs, int innerRowA, int innerColA,
+                             int innerRowB, int innerColB, T &barrier) {
+
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    cuda::memcpy_async(&As[(innerColA * 4 + 0) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 1) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 1],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 2) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 2],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 3) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 3],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    cuda::memcpy_async(&Bs[(innerRowB + offset) * BN + innerColB * 4],
+                       &B[(innerRowB + offset) * N + innerColB * 4],
+                       cuda::aligned_size_t<sizeof(float4)>(sizeof(float4)),
+                       barrier);
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace
+
+/*
+ * @tparam BM The threadblock size for M dimension SMEM caching.
+ * @tparam BN The threadblock size for N dimension SMEM caching.
+ * @tparam BK The threadblock size for K dimension SMEM caching.
+ * @tparam WM M dim of continuous tile computed by each warp
+ * @tparam WN N dim of continuous tile computed by each warp
+ * @tparam WMITER The number of subwarp tiling steps in M dimension.
+ * @tparam WNITER The number of subwarp tiling steps in N dimension.
+ * @tparam TM The per-thread tile size for M dimension.
+ * @tparam TN The per-thread tile size for N dimension.
+ */
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    runSgemmDoubleBuffering2(int M, int N, int K, float alpha, float *A,
+                             float *B, float beta, float *C) {
+  auto block = cooperative_groups::this_thread_block();
+  __shared__ cuda::barrier<cuda::thread_scope::thread_scope_block> frontBarrier;
+  __shared__ cuda::barrier<cuda::thread_scope::thread_scope_block> backBarrier;
+  auto frontBarrierPtr = &frontBarrier;
+  auto backBarrierPtr = &backBarrier;
+  if (block.thread_rank() == 0) {
+    init(&frontBarrier, block.size());
+    init(&backBarrier, block.size());
+  }
+  __syncthreads();
+
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[2 * BM * BK];
+  __shared__ float Bs[2 * BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  int As_offset = 0;
+  int Bs_offset = 0;
+
+  // double-buffering: load first blocktile into SMEM
+  loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+      N, K, A, B, As + As_offset * BM * BK, Bs + Bs_offset * BK * BN, innerRowA,
+      innerColA, innerRowB, innerColB, (*frontBarrierPtr));
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K - BK; bkIdx += BK) {
+    // double-buffering: load next blocktile into SMEM
+    loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A + BK, B + BK * N, As + (1 - As_offset) * BM * BK,
+        Bs + (1 - Bs_offset) * BK * BN, innerRowA, innerColA, innerRowB,
+        innerColB, (*backBarrierPtr));
+
+    // compute the current blocktile
+    (*frontBarrierPtr).arrive_and_wait();
+    processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM, TN>(
+        regM, regN, threadResults, As + As_offset * BM * BK,
+        Bs + Bs_offset * BK * BN, warpRow, warpCol, threadRowInWarp,
+        threadColInWarp);
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    As_offset = 1 - As_offset;
+    Bs_offset = 1 - Bs_offset;
+    // swap the front and back barriers
+    auto tmp = frontBarrierPtr;
+    frontBarrierPtr = backBarrierPtr;
+    backBarrierPtr = tmp;
+
+    __syncthreads();
+  }
+
+  // compute the last blocktile
+  (*frontBarrierPtr).arrive_and_wait();
+  processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM, TN>(
+      regM, regN, threadResults, As + As_offset * BM * BK,
+      Bs + Bs_offset * BK * BN, warpRow, warpCol, threadRowInWarp,
+      threadColInWarp);
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 3.19)
+project(NVIDIA_SGEMM_PRACTICE LANGUAGES CXX CUDA)
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+find_package(CUDA REQUIRED)
+
+# ensure cuda is available
+include(CheckLanguage)
+check_language(CUDA)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CUDA_COMPUTE_CAPABILITY 75)
+
+# in debug mode, add debug symbols to device code
+# this disables most optimizations and kills performance
+add_compile_options("$<$<AND:$<CONFIG:Debug>,$<COMPILE_LANGUAGE:CUDA>>:-G;-src-in-ptx>")
+# add_compile_options("--ptxas-options=-v")
+
+# Configure header file search paths
+include_directories(${CUDA_INCLUDE_DIRS})
+include_directories(${PROJECT_SOURCE_DIR}/src)
+# Configure the source file path to be compiled
+aux_source_directory(${PROJECT_SOURCE_DIR}/src SRC)
+
+# generate executable
+add_executable(sgemm sgemm.cu ${SRC})
+set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
+target_link_libraries(sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
+
+add_executable(cuBLAS_sgemm cuBLAS_sgemm.cu )
+set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
+target_link_libraries(cuBLAS_sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
+
+add_executable(simplest_kernel simplest_kernel.cu)
+set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
+target_link_libraries(simplest_kernel ${CUDA_LIBRARIES})

upstream_ref/sgemm_edtallison/Makefile

@@ -0,0 +1,34 @@
+.PHONY: all build debug clean profile bench cuobjdump
+
+CMAKE := cmake
+
+BUILD_DIR := build
+BENCHMARK_DIR := benchmark_results
+
+all: build
+
+build:
+	@mkdir -p $(BUILD_DIR)
+	@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Release ..
+	@$(MAKE) -C $(BUILD_DIR)
+
+debug:
+	@mkdir -p $(BUILD_DIR)
+	@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Debug ..
+	@$(MAKE) -C $(BUILD_DIR)
+
+clean:
+	@rm -rf $(BUILD_DIR)
+
+FUNCTION := $$(cuobjdump -symbols build/sgemm | grep -i Warptiling | awk '{print $$NF}')
+
+cuobjdump: build
+	@cuobjdump -arch sm_86 -sass -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.sass
+	@cuobjdump -arch sm_86 -ptx -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.ptx
+
+# Usage: make profile KERNEL=<integer> PREFIX=<optional string>
+profile: build
+	@ncu --set full --export $(BENCHMARK_DIR)/$(PREFIX)kernel_$(KERNEL) --force-overwrite $(BUILD_DIR)/sgemm $(KERNEL)
+
+bench: build
+	@bash gen_benchmark_results.sh

upstream_ref/sgemm_edtallison/README.md

@@ -0,0 +1,78 @@
+Reimplementation of Simon Boehm's [CUDA SGEMM](https://github.com/siboehm/SGEMM_CUDA) kernels.
+
+Following the [article](https://siboehm.com/articles/22/CUDA-MMM), for my learning :).
+
+## Run on Google Colab
+
+[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/edtallison/sgemm-cuda/blob/master/run_on_colab.ipynb)
+
+Click the link above to open and run the project in a GPU-enabled Google Colab environment. No additional setup required.
+
+# Notes
+(also scattered throughout kernel code)
+
+## 1. Naive
+
+- **Three-level hierarchy of computation**
+    - Grid, block, thread. Assume grid and block are 2D, thread is the atomic unit of computation.
+    - Blocks can have up to 1024 threads.
+    - Threads within the same block share memory (SMEM).
+
+- **Grid and block indexing**
+    - `gridDim` specifies dimensions of the grid i.e. rows and columns of blocks.
+    - `blockDim` specifies dimensions of the block i.e. rows and columns of threads.
+    - `blockIdx.x/y/z` specifies the block's position in the grid.
+    - `threadIdx.x/y/z `specifies the thread's position in the block.
+    - When used within a kernel, these vars are automatically assigned by the CUDA runtime.
+
+- **Matrix Multiplication**
+    - Matrix multiplication: element ij of C is the dot product of row i of A and column j of B.
+    - In this kernel, each thread computes one element of C. This can obviously be done in parallel so no synchronisation is required.
+
+- **Kernel Launch**
+    - When the kernel is launched, we make the grid as big as necessary to cover all of C, depending on the block size.
+    - The kernel execution is launched asynchronously i.e. the function call on the host (CPU) returns immediately.
+
+- **Memory Access Pattern**
+    - Threads within the same block e.g. `threadIds` (0, 0) and (0, 1) use the same column of B.
+    - They each load the whole column from global memory. Hmmm this seems inefficient...
+
+## 2. Global Memory Coalescing
+
+- **Warps**
+    - In execution, within a block, threads are grouped into "warps" of 32 threads.
+    - Each streaming multiprocessor (SM) has four warp schedulers - physical cores that execute instructions.
+    - Each warp is assigned to a warp scheduler, based on a consecutive `threadId` (x, y, z).
+    - Threads with neighbouring `threadId` become part of the same warp.
+
+- **Global Memory Coalescing**
+    - Sequential memory acceses by threads in the same warp can be grouped and executed as one.
+    - Important to keep in mind when optimising GMEM memory access.
+    - For coalescing, the memory addresses need to be consecutive, but the within-warp accesses don't need to be consecutive.
+    - GPU supports 32B, 64B, and 128B memory accesses.
+
+- **Memory Access Pattern** (this part took me some time to get my head around)
+    - In naive kernel, iterating threads with `threadIdx.x` (which aligns with consecutive `threadId`) actually leads to consecutive threads operating on consecutive rows of A, and the same row of B
+    - If, instead, the threads operated on the same row of A but consecutive columns of B, this accessing of the B values could be coalesced.
+    - This is achieved simply by changing the x and y position indices of the C element computed by each thread.
+    - Note that in either case, we can use within-warp broadcasting as the same row of A or col of B is being accessed by the threads.
+
+## 3. Shared Memory Cache-Blocking
+
+- **SMEM in GPU Memory Architecture**
+    - GPU has global memory GMEM.
+    - Each Streaming Multiprocessor (SM) has a much smaller memory called shared memory (SMEM).
+    - This SMEM is partitioned among the blocks.
+    - Each block of threads runs on a single SM. Multiple blocks can be assigned to the same SM.
+    - A thread can communicate with the other threads in its block via the SMEM chunk.
+    - SMEM, being located on-chip, has much lower latency and higher bandwidth than GMEM.
+
+- **Kernel Memory Access**
+    - Load a chunk of A and a chunk of B from GMEM into SMEM.
+    - Perform as much work as possible on the chunks.
+    - Perform partial sums on C, moving the chunks along the columns of A (same row) and rows of B (same col) until result fully computed.
+    - I.e. in this kernel, each block of threads computes one `BLOCKSIZE*BLOCKSIZE` tile of C.
+
+- **Improvement**
+    - For this kernel, resources mostly spent in waiting for SMEM accesses to return.
+    - Need to make the kernel issue less SMEM instructions to improve efficiency.

upstream_ref/sgemm_edtallison/cuBLAS_sgemm.cu

@@ -0,0 +1,108 @@
+#include <cstdio>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+/*
+ * A stand-alone script to invoke & benchmark standard cuBLAS SGEMM performance
+ */
+
+int main(int argc, char *argv[]) {
+  int m = 2;
+  int k = 3;
+  int n = 4;
+  int print = 1;
+  cudaError_t cudaStat;  // cudaMalloc status
+  cublasStatus_t stat;   // cuBLAS functions status
+  cublasHandle_t handle; // cuBLAS context
+
+  int i, j;
+
+  float *a, *b, *c;
+
+  // malloc for a,b,c...
+  a = (float *)malloc(m * k * sizeof(float));
+  b = (float *)malloc(k * n * sizeof(float));
+  c = (float *)malloc(m * n * sizeof(float));
+
+  int ind = 11;
+  for (j = 0; j < m * k; j++) {
+    a[j] = (float)ind++;
+  }
+
+  ind = 11;
+  for (j = 0; j < k * n; j++) {
+    b[j] = (float)ind++;
+  }
+
+  ind = 11;
+  for (j = 0; j < m * n; j++) {
+    c[j] = (float)ind++;
+  }
+
+  // DEVICE
+  float *d_a, *d_b, *d_c;
+
+  // cudaMalloc for d_a, d_b, d_c...
+  cudaMalloc((void **)&d_a, m * k * sizeof(float));
+  cudaMalloc((void **)&d_b, k * n * sizeof(float));
+  cudaMalloc((void **)&d_c, m * n * sizeof(float));
+
+  stat = cublasCreate(&handle); // initialize CUBLAS context
+
+  cudaMemcpy(d_a, a, m * k * sizeof(float), cudaMemcpyHostToDevice);
+  cudaMemcpy(d_b, b, k * n * sizeof(float), cudaMemcpyHostToDevice);
+  cudaMemcpy(d_c, c, m * n * sizeof(float), cudaMemcpyHostToDevice);
+
+  float alpha = 1.0f;
+  float beta = 0.5f;
+
+  if (print == 1) {
+    printf("alpha = %4.0f, beta = %4.0f\n", alpha, beta);
+    printf("A = (mxk: %d x %d)\n", m, k);
+    for (i = 0; i < m; i++) {
+      for (j = 0; j < k; j++) {
+        printf("%4.1f ", a[i * m + j]);
+      }
+      printf("\n");
+    }
+    printf("B = (kxn: %d x %d)\n", k, n);
+    for (i = 0; i < k; i++) {
+      for (j = 0; j < n; j++) {
+        printf("%4.1f ", b[i * n + j]);
+      }
+      printf("\n");
+    }
+    printf("C = (mxn: %d x %d)\n", m, n);
+    for (i = 0; i < m; i++) {
+      for (j = 0; j < n; j++) {
+        printf("%4.1f ", c[i * n + j]);
+      }
+      printf("\n");
+    }
+  }
+
+  stat = cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, d_b, n,
+                     d_a, k, &beta, d_c, n);
+
+  cudaMemcpy(c, d_c, m * n * sizeof(float), cudaMemcpyDeviceToHost);
+
+  if (print == 1) {
+    printf("\nC after SGEMM = \n");
+    for (i = 0; i < m; i++) {
+      for (j = 0; j < n; j++) {
+        printf("%4.1f ", c[i * n + j]);
+      }
+      printf("\n");
+    }
+  }
+
+  cudaFree(d_a);
+  cudaFree(d_b);
+  cudaFree(d_c);
+  cublasDestroy(handle); // destroy CUBLAS context
+  free(a);
+  free(b);
+  free(c);
+
+  return EXIT_SUCCESS;
+}

upstream_ref/sgemm_edtallison/kernels.cuh

@@ -0,0 +1,14 @@
+#pragma once
+
+#include "kernels/01_naive.cuh"
+#include "kernels/02_kernel_global_mem_coalesce.cuh"
+#include "kernels/03_kernel_shared_mem_blocking.cuh"
+#include "kernels/04_kernel_1D_blocktiling.cuh"
+#include "kernels/05_kernel_2D_blocktiling.cuh"
+#include "kernels/06_kernel_vectorize.cuh"
+#include "kernels/07_kernel_resolve_bank_conflicts.cuh"
+#include "kernels/08_kernel_bank_extra_col.cuh"
+#include "kernels/09_kernel_autotuned.cuh"
+#include "kernels/10_kernel_warptiling.cuh"
+#include "kernels/11_kernel_double_buffering.cuh"
+#include "kernels/12_kernel_double_buffering.cuh"

upstream_ref/sgemm_edtallison/run_on_colab.ipynb

@@ -0,0 +1,71 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "id": "b9326784",
+      "metadata": {},
+      "source": [
+        "This notebook installs dependencies, builds the project, and runs a selected kernel on Colab's GPU.\n",
+        "\n",
+        "**Note:** Ensure Colab runtime type is set to GPU (Runtime → Change runtime type → GPU)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "35ef02c6",
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "# install system dependencies\n",
+        "!apt-get update && apt-get install -y cmake ninja-build"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "37a03120",
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "# clone\n",
+        "!git clone https://github.com/edtallison/sgemm-cuda.git\n",
+        "%cd sgemm-cuda"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "cd1855de",
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "# build\n",
+        "!mkdir -p build && cd build && cmake -G Ninja .. && ninja"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "093a0b6d",
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "# run for kernel 1\n",
+        "!cd build && ./sgemm 1"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "name": "python"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}

upstream_ref/sgemm_edtallison/runner.cu

@@ -0,0 +1,549 @@
+#include "kernels.cuh"
+#include "runner.cuh"
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <iomanip>
+
+float get_sec() {
+  struct timeval time;
+  gettimeofday(&time, NULL);
+  return (1e6 * time.tv_sec + time.tv_usec);
+}
+
+float cpu_elapsed_time(float &beg, float &end) { return 1.0e-6 * (end - beg); }
+
+void cudaCheck(cudaError_t error, const char *file, int line) {
+  if (error != cudaSuccess) {
+    printf("[CUDA ERROR] at file %s:%d:\n%s\n", file, line,
+           cudaGetErrorString(error));
+    exit(EXIT_FAILURE);
+  }
+};
+
+void CudaDeviceInfo() {
+  int deviceId;
+
+  cudaGetDevice(&deviceId);
+
+  cudaDeviceProp props{};
+  cudaGetDeviceProperties(&props, deviceId);
+
+  printf("Device ID: %d\n\
+    Name: %s\n\
+    Compute Capability: %d.%d\n\
+    memoryBusWidth: %d\n\
+    maxThreadsPerBlock: %d\n\
+    maxThreadsPerMultiProcessor: %d\n\
+    maxRegsPerBlock: %d\n\
+    maxRegsPerMultiProcessor: %d\n\
+    totalGlobalMem: %zuMB\n\
+    sharedMemPerBlock: %zuKB\n\
+    sharedMemPerMultiprocessor: %zuKB\n\
+    totalConstMem: %zuKB\n\
+    multiProcessorCount: %d\n\
+    Warp Size: %d\n",
+         deviceId, props.name, props.major, props.minor, props.memoryBusWidth,
+         props.maxThreadsPerBlock, props.maxThreadsPerMultiProcessor,
+         props.regsPerBlock, props.regsPerMultiprocessor,
+         props.totalGlobalMem / 1024 / 1024, props.sharedMemPerBlock / 1024,
+         props.sharedMemPerMultiprocessor / 1024, props.totalConstMem / 1024,
+         props.multiProcessorCount, props.warpSize);
+};
+
+void randomize_matrix(float *mat, int N) {
+  // NOTICE: Use gettimeofday instead of srand((unsigned)time(NULL)); the time
+  // precision is too low and the same random number is generated.
+  struct timeval time {};
+  gettimeofday(&time, nullptr);
+  srand(time.tv_usec);
+  for (int i = 0; i < N; i++) {
+    float tmp = (float)(rand() % 5) + 0.01 * (rand() % 5);
+    tmp = (rand() % 2 == 0) ? tmp : tmp * (-1.);
+    mat[i] = tmp;
+  }
+}
+
+void range_init_matrix(float *mat, int N) {
+  for (int i = 0; i < N; i++) {
+    mat[i] = i;
+  }
+}
+
+void zero_init_matrix(float *mat, int N) {
+  for (int i = 0; i < N; i++) {
+    mat[i] = 0.0;
+  }
+}
+
+void copy_matrix(const float *src, float *dest, int N) {
+  int i;
+  for (i = 0; src + i && dest + i && i < N; i++)
+    *(dest + i) = *(src + i);
+  if (i != N)
+    printf("copy failed at %d while there are %d elements in total.\n", i, N);
+}
+
+void print_matrix(const float *A, int M, int N, std::ofstream &fs) {
+  int i;
+  fs << std::setprecision(2)
+     << std::fixed; // Set floating-point precision and fixed notation
+  fs << "[";
+  for (i = 0; i < M * N; i++) {
+    if ((i + 1) % N == 0)
+      fs << std::setw(5) << A[i]; // Set field width and write the value
+    else
+      fs << std::setw(5) << A[i] << ", ";
+    if ((i + 1) % N == 0) {
+      if (i + 1 < M * N)
+        fs << ";\n";
+    }
+  }
+  fs << "]\n";
+}
+
+bool verify_matrix(float *matRef, float *matOut, int N) {
+  double diff = 0.0;
+  int i;
+  for (i = 0; i < N; i++) {
+    diff = std::fabs(matRef[i] - matOut[i]);
+    if (diff > 0.01) {
+      printf("Divergence! Should %5.2f, Is %5.2f (Diff %5.2f) at %d\n",
+             matRef[i], matOut[i], diff, i);
+      return false;
+    }
+  }
+  return true;
+}
+
+int div_ceil(int numerator, int denominator) {
+  std::div_t res = std::div(numerator, denominator);
+  return res.rem ? (res.quot + 1) : res.quot;
+}
+
+void runCublasFP32(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // cuBLAS uses column-major order. So we change the order of our row-major A &
+  // B, since (B^T*A^T)^T = (A*B)
+  // This runs cuBLAS in full fp32 mode
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N, CUBLAS_COMPUTE_32F,
+               CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void runCublasBF16(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // This runs cuBLAS with mixed precision (performing the mul with operands
+  // downcast to bf16), which is ~4x faster
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N,
+               CUBLAS_COMPUTE_32F_FAST_16BF, CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void runCublasTF32(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // This runs cuBLAS with mixed precision (performing the mul with operands
+  // downcast to bf16), which is ~4x faster
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N,
+               CUBLAS_COMPUTE_32F_FAST_TF32, CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void run_sgemm_naive(int M, int N, int K, float alpha, float *A, float *B,
+                     float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32, 32);
+  sgemm_naive<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_sgemm_coalesce(int M, int N, int K, float alpha, float *A, float *B,
+                        float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32 * 32);
+  sgemm_global_mem_coalesce<32>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_sgemm_shared_mem_block(int M, int N, int K, float alpha, float *A,
+                                float *B, float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32 * 32);
+  // L1 cache becomes useless, since we access GMEM only via SMEM, so we carve
+  // out all of L1 to SMEM. This doesn't currently make a difference, since
+  // occupancy is limited by reg and thread count, but it's good to do anyway.
+  cudaFuncSetAttribute(sgemm_shared_mem_block<32>,
+                       cudaFuncAttributePreferredSharedMemoryCarveout,
+                       cudaSharedmemCarveoutMaxShared);
+  sgemm_shared_mem_block<32>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemm1DBlocktiling(int M, int N, int K, float alpha, float *A, float *B,
+                           float beta, float *C) {
+  const uint BM = 64;
+  const uint BN = 64;
+  const uint BK = 8;
+  const uint TM = 8;
+  dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+  dim3 blockDim((BM * BN) / TM);
+  sgemm1DBlocktiling<BM, BN, BK, TM>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemm2DBlocktiling(int M, int N, int K, float alpha, float *A, float *B,
+                           float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemm2DBlocktiling<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemm2DBlocktiling<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmVectorize(int M, int N, int K, float alpha, float *A, float *B,
+                       float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmVectorize<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmVectorize<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmResolveBankConflicts(int M, int N, int K, float alpha, float *A,
+                                  float *B, float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankConflicts<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankConflicts<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmResolveBankExtraCol(int M, int N, int K, float alpha, float *A,
+                                 float *B, float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankExtraCol<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankExtraCol<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
+                       float beta, float *C) {
+  // A100
+  // const uint K9_BK = 16;
+  // const uint K9_TM = 4;
+  // const uint K9_TN = 4;
+  // const uint K9_BM = 64;
+  // const uint K9_BN = 64;
+  // A6000
+  const uint K9_BK = 16;
+  const uint K9_TM = 8;
+  const uint K9_TN = 8;
+  const uint K9_BM = 128;
+  const uint K9_BN = 128;
+  dim3 blockDim(K9_NUM_THREADS);
+
+  static_assert(
+      (K9_NUM_THREADS * 4) % K9_BK == 0,
+      "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization issues "
+      "during GMEM->SMEM tiling (loading only parts of the final row of Bs "
+      "during each iteraion)");
+  static_assert(
+      (K9_NUM_THREADS * 4) % K9_BN == 0,
+      "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization issues "
+      "during GMEM->SMEM tiling (loading only parts of the final row of As "
+      "during each iteration)");
+  static_assert(
+      K9_BN % (16 * K9_TN) == 0,
+      "K9_BN must be a multiple of 16*K9_TN to avoid quantization effects");
+  static_assert(
+      K9_BM % (16 * K9_TM) == 0,
+      "K9_BM must be a multiple of 16*K9_TM to avoid quantization effects");
+  static_assert((K9_BM * K9_BK) % (4 * K9_NUM_THREADS) == 0,
+                "K9_BM*K9_BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K9_BN * K9_BK) % (4 * K9_NUM_THREADS) == 0,
+                "K9_BN*K9_BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K9_BN), CEIL_DIV(M, K9_BM));
+  sgemmAutotuned<K9_BM, K9_BN, K9_BK, K9_TM, K9_TN>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmWarptiling(int M, int N, int K, float alpha, float *A, float *B,
+                        float beta, float *C) {
+  // Settings for A100
+  // const uint K10_NUM_THREADS = 128;
+  // const uint K10_BN = 128;
+  // const uint K10_BM = 64;
+  // const uint K10_BK = 16;
+  // const uint K10_WN = 64;
+  // const uint K10_WM = 32;
+  // const uint K10_WNITER = 1;
+  // const uint K10_TN = 4;
+  // const uint K10_TM = 4;
+  // Settings for A6000
+  const uint K10_NUM_THREADS = 128;
+  const uint K10_BN = 128;
+  const uint K10_BM = 128;
+  const uint K10_BK = 16;
+  const uint K10_WN = 64;
+  const uint K10_WM = 64;
+  const uint K10_WNITER = 4;
+  const uint K10_TN = 4;
+  const uint K10_TM = 8;
+  dim3 blockDim(K10_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K10_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K10_BN % K10_WN == 0) and (K10_BM % K10_WM == 0));
+  static_assert((K10_BN / K10_WN) * (K10_BM / K10_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K10_WM * K10_WN) % (WARPSIZE * K10_TM * K10_TN * K10_WNITER) ==
+                0);
+  constexpr uint K10_WMITER =
+      (K10_WM * K10_WN) / (32 * K10_TM * K10_TN * K10_WNITER);
+  // warpsubtile in warptile
+  static_assert((K10_WM % K10_WMITER == 0) and (K10_WN % K10_WNITER == 0));
+
+  static_assert((K10_NUM_THREADS * 4) % K10_BK == 0,
+                "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K10_NUM_THREADS * 4) % K10_BN == 0,
+                "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K10_BN % (16 * K10_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K10_BM % (16 * K10_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K10_BM * K10_BK) % (4 * K10_NUM_THREADS) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K10_BN * K10_BK) % (4 * K10_NUM_THREADS) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K10_BN), CEIL_DIV(M, K10_BM));
+  sgemmWarptiling<K10_BM, K10_BN, K10_BK, K10_WM, K10_WN, K10_WNITER, K10_TM,
+                  K10_TN, K10_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmDoubleBuffering(int M, int N, int K, float alpha, float *A,
+                             float *B, float beta, float *C) {
+  // Settings for A100
+  // const uint K11_NUM_THREADS = 256;
+  // const uint K11_BN = 128;
+  // const uint K11_BM = 64;
+  // const uint K11_BK = 16;
+  // const uint K11_WN = 32;
+  // const uint K11_WM = 32;
+  // const uint K11_WNITER = 2;
+  // const uint K11_TN = 4;
+  // const uint K11_TM = 4;
+  // Settings for A6000
+  const uint K11_NUM_THREADS = 256;
+  const uint K11_BN = 256;
+  const uint K11_BM = 128;
+  const uint K11_BK = 16;
+  const uint K11_WN = 32;
+  const uint K11_WM = 128;
+  const uint K11_WNITER = 1;
+  const uint K11_TN = 8;
+  const uint K11_TM = 8;
+  dim3 blockDim(K11_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K11_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K11_BN % K11_WN == 0) and (K11_BM % K11_WM == 0));
+  static_assert((K11_BN / K11_WN) * (K11_BM / K11_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K11_WM * K11_WN) % (WARPSIZE * K11_TM * K11_TN * K11_WNITER) ==
+                0);
+  constexpr uint K11_WMITER =
+      (K11_WM * K11_WN) / (32 * K11_TM * K11_TN * K11_WNITER);
+  // warpsubtile in warptile
+  static_assert((K11_WM % K11_WMITER == 0) and (K11_WN % K11_WNITER == 0));
+
+  static_assert((K11_NUM_THREADS / 2 * 4) % K11_BK == 0,
+                "NUM_THREADS*4 must be multiple of BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K11_NUM_THREADS / 2 * 4) % K11_BN == 0,
+                "NUM_THREADS*4 must be multiple of BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K11_BN % (16 * K11_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K11_BM % (16 * K11_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K11_BM * K11_BK) % (4 * K11_NUM_THREADS / 2) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K11_BN * K11_BK) % (4 * K11_NUM_THREADS / 2) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K11_BN), CEIL_DIV(M, K11_BM));
+  sgemmDoubleBuffering<K11_BM, K11_BN, K11_BK, K11_WM, K11_WN, K11_WNITER,
+                       K11_TM, K11_TN, K11_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmDoubleBuffering2(int M, int N, int K, float alpha, float *A,
+                              float *B, float beta, float *C) {
+  // Settings for A6000
+  const uint K12_NUM_THREADS = 128;
+  const uint K12_BN = 128;
+  const uint K12_BM = 128;
+  const uint K12_BK = 16;
+  const uint K12_WN = 64;
+  const uint K12_WM = 64;
+  const uint K12_WNITER = 4;
+  const uint K12_TN = 4;
+  const uint K12_TM = 8;
+  dim3 blockDim(K12_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K12_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K12_BN % K12_WN == 0) and (K12_BM % K12_WM == 0));
+  static_assert((K12_BN / K12_WN) * (K12_BM / K12_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K12_WM * K12_WN) % (WARPSIZE * K12_TM * K12_TN * K12_WNITER) ==
+                0);
+  constexpr uint K12_WMITER =
+      (K12_WM * K12_WN) / (32 * K12_TM * K12_TN * K12_WNITER);
+  // warpsubtile in warptile
+  static_assert((K12_WM % K12_WMITER == 0) and (K12_WN % K12_WNITER == 0));
+
+  static_assert((K12_NUM_THREADS * 4) % K12_BK == 0,
+                "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K12_NUM_THREADS * 4) % K12_BN == 0,
+                "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K12_BN % (16 * K12_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K12_BM % (16 * K12_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K12_BM * K12_BK) % (4 * K12_NUM_THREADS) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K12_BN * K12_BK) % (4 * K12_NUM_THREADS) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K12_BN), CEIL_DIV(M, K12_BM));
+  runSgemmDoubleBuffering2<K12_BM, K12_BN, K12_BK, K12_WM, K12_WN, K12_WNITER,
+                           K12_TM, K12_TN, K12_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_kernel(int kernel_num, int M, int N, int K, float alpha, float *A,
+                float *B, float beta, float *C, cublasHandle_t handle) {
+  switch (kernel_num) {
+  case 0:
+    runCublasFP32(handle, M, N, K, alpha, A, B, beta, C);
+    break;
+  case 1:
+    run_sgemm_naive(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 2:
+    run_sgemm_coalesce(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 3:
+    run_sgemm_shared_mem_block(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 4:
+    runSgemm1DBlocktiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 5:
+    runSgemm2DBlocktiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 6:
+    runSgemmVectorize(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 7:
+    runSgemmResolveBankConflicts(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 8:
+    runSgemmResolveBankExtraCol(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 9:
+    runSgemmAutotuned(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 10:
+    runSgemmWarptiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 11:
+    runSgemmDoubleBuffering(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 12:
+    runSgemmDoubleBuffering2(M, N, K, alpha, A, B, beta, C);
+    break;
+  default:
+    throw std::invalid_argument("Unknown kernel number");
+  }
+}

upstream_ref/sgemm_edtallison/runner.cuh

@@ -0,0 +1,26 @@
+#pragma once
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+#include <fstream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <time.h>
+#include <unistd.h>
+
+void cudaCheck(cudaError_t error, const char *file,
+               int line); // CUDA error check
+void CudaDeviceInfo();    // print CUDA information
+
+void range_init_matrix(float *mat, int N);
+void randomize_matrix(float *mat, int N);
+void zero_init_matrix(float *mat, int N);
+void copy_matrix(const float *src, float *dest, int N);
+void print_matrix(const float *A, int M, int N, std::ofstream &fs);
+bool verify_matrix(float *mat1, float *mat2, int N);
+
+float get_current_sec();                        // Get the current moment
+float cpu_elapsed_time(float &beg, float &end); // Calculate time difference
+
+void run_kernel(int kernel_num, int m, int n, int k, float alpha, float *A,
+                float *B, float beta, float *C, cublasHandle_t handle);

upstream_ref/sgemm_edtallison/scripts/bank_calc.py

@@ -0,0 +1,33 @@
+banks_naive = lambda r, c: (r * 32 + c) % 32
+banks_one_extra = lambda r, c: (r * 33 + c) % 32
+
+ITEMS_PER_WARP = 8
+
+
+def printBankConflicts(bank_fun):
+    for c in range(1):
+        banks = []
+        for i in range(32):
+            row = (i * ITEMS_PER_WARP) // 16 
+            col = (i * ITEMS_PER_WARP + c) % 16
+            banks.append((i, row, col, bank_fun(row, col)))
+        print("Step", c, "\n", "\n".join(["(" + ",".join(str(x) for x in i) + ")" for i in banks]))
+        d = {k: 0 for k in range(32)}
+        for i in banks:
+            d[i[-1]] += 1
+
+        count = 0
+        for key, val in d.items():
+            if val > 0:
+                count += 1
+
+        print(
+            f"Bank conflicts (Step {c}): {sorted(d.items(), key=lambda item: item[1], reverse=True)[0][1]}, banks accessed: {count}/32\n"
+        )
+
+
+print("---NAIVE---")
+printBankConflicts(banks_naive, 32)
+
+print("\n---EXTRA COL---")
+printBankConflicts(banks_one_extra, 33)

upstream_ref/sgemm_edtallison/scripts/kernel_10_autotuner.sh

@@ -0,0 +1,115 @@
+#!/usr/bin/env bash
+
+set -u 
+
+# Define the range of values for each parameter
+BK_VALUES=(8 16 32 64)
+BM_VALUES=(64 128 256)
+BN_VALUES=(64 128 256)
+WM_VALUES=(32 64 128 256)
+WN_VALUES=(32 64 128 256)
+WNITER_VALUES=(1 2 4 8)
+TM_VALUES=(4 8 16 32)
+TN_VALUES=(4 8 16 32)
+NUM_THREADS_VALUES=(128 256)
+
+cd "$(dirname "$0")"
+cd "../build"
+
+RUNNER="../src/runner.cu"
+OUTPUT="../benchmark_results/kernel_10_autotune_results.txt"
+
+# Clear the output file
+echo "" > $OUTPUT
+
+# Set GPU to use
+export DEVICE="0"
+WARPSIZE=32
+
+
+TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
+CONFIG_NUM=0
+
+# Loop through all combinations of parameters
+for BK in "${BK_VALUES[@]}"; do
+for BM in "${BM_VALUES[@]}"; do
+for BN in "${BN_VALUES[@]}"; do
+for WM in "${WM_VALUES[@]}"; do
+for WN in "${WN_VALUES[@]}"; do
+for WN_ITER in "${WNITER_VALUES[@]}"; do
+for TM in "${TM_VALUES[@]}"; do
+for TN in "${TN_VALUES[@]}"; do
+for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
+echo ""
+CONFIG_NUM=$(( CONFIG_NUM + 1 ))
+# skip configurations that don't fullfil preconditions
+NUM_WARPS=$(( NUM_THREADS / 32 ))
+if ! (( BN % WN == 0 && BM % WM == 0 )); then
+  echo "Error: BN % WN must be 0 and BM % WM must be 0."
+  continue
+fi
+if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
+  echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
+  continue
+fi
+if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
+  echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
+  continue
+fi
+WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
+if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
+  echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
+  continue
+fi
+if ! (( (NUM_THREADS * 4) % BK == 0 )); then
+  echo "Error: (NUM_THREADS * 4) % BK must be 0."
+  continue
+fi
+if ! (( (NUM_THREADS * 4) % BN == 0 )); then
+  echo "Error: (NUM_THREADS * 4) % BN must be 0."
+  continue
+fi
+if ! (( BN % (16 * TN) == 0 )); then
+  echo "Error: BN must be a multiple of 16 * TN."
+  continue
+fi
+if ! (( BM % (16 * TM) == 0 )); then
+  echo "Error: BM must be a multiple of 16 * TM."
+  continue
+fi
+if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
+  echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
+  continue
+fi
+if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
+  echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
+  continue
+fi
+
+# Update the parameters in the source code
+sed -i "s/const uint K10_NUM_THREADS = .*/const uint K10_NUM_THREADS = $NUM_THREADS;/" $RUNNER
+sed -i "s/const uint K10_BN = .*/const uint K10_BN = $BN;/" $RUNNER
+sed -i "s/const uint K10_BM = .*/const uint K10_BM = $BM;/" $RUNNER
+sed -i "s/const uint K10_BK = .*/const uint K10_BK = $BK;/" $RUNNER
+sed -i "s/const uint K10_WM = .*/const uint K10_WM = $WM;/" $RUNNER
+sed -i "s/const uint K10_WN = .*/const uint K10_WN = $WN;/" $RUNNER
+sed -i "s/const uint K10_WNITER = .*/const uint K10_WNITER = $WN_ITER;/" $RUNNER
+sed -i "s/const uint K10_TM = .*/const uint K10_TM = $TM;/" $RUNNER
+sed -i "s/const uint K10_TN = .*/const uint K10_TN = $TN;/" $RUNNER
+
+# Rebuild the program
+make
+
+echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
+# Run the benchmark and get the result
+# Kill the program after 4 seconds if it doesn't finish
+timeout -v 8 ./sgemm 10 | tee -a $OUTPUT
+done
+done
+done
+done
+done
+done
+done
+done
+done

upstream_ref/sgemm_edtallison/scripts/kernel_11_autotuner.sh

@@ -0,0 +1,115 @@
+#!/usr/bin/env bash
+
+set -u 
+
+# Define the range of values for each parameter
+BK_VALUES=(8 16 32 64)
+BM_VALUES=(64 128 256)
+BN_VALUES=(64 128 256)
+WM_VALUES=(32 64 128 256)
+WN_VALUES=(32 64 128 256)
+WNITER_VALUES=(1 2 4 8)
+TM_VALUES=(4 8 16 32)
+TN_VALUES=(4 8 16 32)
+NUM_THREADS_VALUES=(128 256)
+
+cd "$(dirname "$0")"
+cd "../build"
+
+RUNNER="../src/runner.cu"
+OUTPUT="../benchmark_results/kernel_11_autotune_results.txt"
+
+# Clear the output file
+echo "" > $OUTPUT
+
+# Set GPU to use
+export DEVICE="0"
+WARPSIZE=32
+
+
+TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
+CONFIG_NUM=0
+
+# Loop through all combinations of parameters
+for BK in "${BK_VALUES[@]}"; do
+for BM in "${BM_VALUES[@]}"; do
+for BN in "${BN_VALUES[@]}"; do
+for WM in "${WM_VALUES[@]}"; do
+for WN in "${WN_VALUES[@]}"; do
+for WN_ITER in "${WNITER_VALUES[@]}"; do
+for TM in "${TM_VALUES[@]}"; do
+for TN in "${TN_VALUES[@]}"; do
+for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
+echo ""
+CONFIG_NUM=$(( CONFIG_NUM + 1 ))
+# skip configurations that don't fullfil preconditions
+NUM_WARPS=$(( NUM_THREADS / 32 ))
+if ! (( BN % WN == 0 && BM % WM == 0 )); then
+  echo "Error: BN % WN must be 0 and BM % WM must be 0."
+  continue
+fi
+if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
+  echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
+  continue
+fi
+if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
+  echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
+  continue
+fi
+WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
+if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
+  echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
+  continue
+fi
+if ! (( (NUM_THREADS * 4) % BK == 0 )); then
+  echo "Error: (NUM_THREADS * 4) % BK must be 0."
+  continue
+fi
+if ! (( (NUM_THREADS * 4) % BN == 0 )); then
+  echo "Error: (NUM_THREADS * 4) % BN must be 0."
+  continue
+fi
+if ! (( BN % (16 * TN) == 0 )); then
+  echo "Error: BN must be a multiple of 16 * TN."
+  continue
+fi
+if ! (( BM % (16 * TM) == 0 )); then
+  echo "Error: BM must be a multiple of 16 * TM."
+  continue
+fi
+if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
+  echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
+  continue
+fi
+if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
+  echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
+  continue
+fi
+
+# Update the parameters in the source code
+sed -i "s/const uint K11_NUM_THREADS = .*/const uint K11_NUM_THREADS = $NUM_THREADS;/" $RUNNER
+sed -i "s/const uint K11_BN = .*/const uint K11_BN = $BN;/" $RUNNER
+sed -i "s/const uint K11_BM = .*/const uint K11_BM = $BM;/" $RUNNER
+sed -i "s/const uint K11_BK = .*/const uint K11_BK = $BK;/" $RUNNER
+sed -i "s/const uint K11_WM = .*/const uint K11_WM = $WM;/" $RUNNER
+sed -i "s/const uint K11_WN = .*/const uint K11_WN = $WN;/" $RUNNER
+sed -i "s/const uint K11_WNITER = .*/const uint K11_WNITER = $WN_ITER;/" $RUNNER
+sed -i "s/const uint K11_TM = .*/const uint K11_TM = $TM;/" $RUNNER
+sed -i "s/const uint K11_TN = .*/const uint K11_TN = $TN;/" $RUNNER
+
+# Rebuild the program
+make
+
+echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
+# Run the benchmark and get the result
+# Kill the program after 8 seconds if it doesn't finish
+timeout -v 8 ./sgemm 11 | tee -a $OUTPUT
+done
+done
+done
+done
+done
+done
+done
+done
+done

upstream_ref/sgemm_edtallison/scripts/kernel_9_autotuner.sh

@@ -0,0 +1,86 @@
+#!/usr/bin/env bash
+
+set -u
+
+# Define the range of values for each parameter
+BK_VALUES=(8 16 32 64)
+TM_VALUES=(4 8 16 32)
+TN_VALUES=(4 8 16 32)
+BM_VALUES=(64 128 256)
+BN_VALUES=(64 128 256)
+NUM_THREADS_VALUES=(256)
+
+cd "$(dirname "$0")"
+cd "../build"
+
+RUNNER="../src/runner.cu"
+KERNEL="../src/kernels/9_kernel_autotuned.cuh"
+OUTPUT="../benchmark_results/kernel_9_autotune_results.txt"
+
+# Clear the output file
+echo "" > $OUTPUT
+
+# Set GPU to use
+export DEVICE="2"
+
+TOTAL_CONFIGS="$(( ${#NUM_THREADS_VALUES[@]} * ${#BK_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} ))"
+CONFIG_NUM=0
+
+# Loop through all combinations of parameters
+for bk in ${BK_VALUES[@]}; do
+  for tm in ${TM_VALUES[@]}; do
+    for tn in ${TN_VALUES[@]}; do
+      for bm in ${BM_VALUES[@]}; do
+        for bn in ${BN_VALUES[@]}; do
+          for nt in ${NUM_THREADS_VALUES[@]}; do
+            echo ""
+            CONFIG_NUM=$(( $CONFIG_NUM + 1 ))
+
+            # skip configurations that don't fullfil preconditions
+            config="BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NT=$nt"
+            if [[ $(( ($nt * 4) % bk )) -ne 0 ]]; then
+              echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BK = $(( ($nt * 4) % bk )) != 0))"
+              continue
+            fi
+            if [[ $(( ($nt * 4) % bn )) -ne 0 ]]; then
+              echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BN = $(( ($nt * 4) % bn )) != 0))"
+              continue
+            fi
+            if [[ $(( $bn % (16 * $tn ) )) -ne 0 ]]; then
+              echo "QUANTIZATION: Skipping $config because BN % (16 * TN) = $(( $bn % (16 * $tn ) )) != 0))"
+              continue
+            fi
+            if [[ $(( $bm % (16 * $tm ) )) -ne 0 ]]; then
+              echo "QUANTIZATION: Skipping $config because BM % (16 * TM) = $(( $bm % (16 * $tm ) )) != 0))"
+              continue
+            fi
+            if [[ $(( ($bm * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
+              echo "VECTORIZE: Skipping $config because (BM * BK) % (4 * NUM_THREADS) = $(( ($bm * $bk) % ( 4 * 256 ) )) != 0))"
+              continue
+            fi
+            if [[ $(( ($bn * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
+              echo "VECTORIZE: Skipping $config because (BN * BK) % (4 * NUM_THREADS) = $(( ($bn * $bk) % ( 4 * 256 ) )) != 0))"
+              continue
+            fi
+
+            # Update the parameters in the source code
+            sed -i "s/const uint K9_BK = .*/const uint K9_BK = $bk;/" $RUNNER
+            sed -i "s/const uint K9_TM = .*/const uint K9_TM = $tm;/" $RUNNER
+            sed -i "s/const uint K9_TN = .*/const uint K9_TN = $tn;/" $RUNNER
+            sed -i "s/const uint K9_BM = .*/const uint K9_BM = $bm;/" $RUNNER
+            sed -i "s/const uint K9_BN = .*/const uint K9_BN = $bn;/" $RUNNER
+            sed -i "s/const int K9_NUM_THREADS = .*/const int K9_NUM_THREADS = $nt;/" $KERNEL
+            
+            # Rebuild the program
+            make 
+
+            echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NUM_THREADS=$nt" |& tee -a $OUTPUT
+            # Run the benchmark and get the result
+            # Kill the program after 4 seconds if it doesn't finish
+            timeout -v 4 ./sgemm 9 | tee -a $OUTPUT
+          done
+        done
+      done
+    done
+  done
+done

upstream_ref/sgemm_edtallison/sgemm.cu

@@ -0,0 +1,168 @@
+#include <cstdio>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <runner.cuh>
+#include <vector>
+
+#define cudaCheck(err) (cudaCheck(err, __FILE__, __LINE__))
+
+const std::string errLogFile = "matrixValidationFailure.txt";
+
+int main(int argc, char **argv) {
+  if (argc != 2) {
+    std::cerr << "Please select a kernel (range 0 - 12, 0 for NVIDIA cuBLAS)"
+              << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // get kernel number
+  int kernel_num = std::stoi(argv[1]);
+  if (kernel_num < 0 || kernel_num > 12) {
+    std::cerr << "Please enter a valid kernel number (0-12)" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // get environment variable for device
+  int deviceIdx = 0;
+  if (getenv("DEVICE") != NULL) {
+    deviceIdx = atoi(getenv("DEVICE"));
+  }
+  cudaCheck(cudaSetDevice(deviceIdx));
+
+  printf("Running kernel %d on device %d.\n", kernel_num, deviceIdx);
+
+  // print some device info
+  // CudaDeviceInfo();
+
+  // Declare the handle, create the handle, cublasCreate will return a value of
+  // type cublasStatus_t to determine whether the handle was created
+  // successfully (the value is 0)
+  cublasHandle_t handle;
+  if (cublasCreate(&handle)) {
+    std::cerr << "Create cublas handle error." << std::endl;
+    exit(EXIT_FAILURE);
+  };
+
+  // Using cudaEvent for gpu stream timing, cudaEvent is equivalent to
+  // publishing event tasks in the target stream
+  float elapsed_time;
+  cudaEvent_t beg, end;
+  cudaEventCreate(&beg);
+  cudaEventCreate(&end);
+
+  // cuBLAS FLOPs ceiling is reached at 8192
+  std::vector<int> SIZE = {128, 256, 512, 1024, 2048, 4096};
+
+  long m, n, k, max_size;
+  max_size = SIZE[SIZE.size() - 1];
+  std::cout << "Max size: " << max_size << std::endl;
+
+  float alpha = 0.5, beta = 3.0; // GEMM input parameters, C=α*AB+β*C
+
+  float *A = nullptr, *B = nullptr, *C = nullptr,
+        *C_ref = nullptr; // host matrices
+  float *dA = nullptr, *dB = nullptr, *dC = nullptr,
+        *dC_ref = nullptr; // device matrices
+
+  A = (float *)malloc(sizeof(float) * max_size * max_size);
+  B = (float *)malloc(sizeof(float) * max_size * max_size);
+  C = (float *)malloc(sizeof(float) * max_size * max_size);
+  C_ref = (float *)malloc(sizeof(float) * max_size * max_size);
+
+  randomize_matrix(A, max_size * max_size);
+  randomize_matrix(B, max_size * max_size);
+  randomize_matrix(C, max_size * max_size);
+
+  cudaCheck(cudaMalloc((void **)&dA, sizeof(float) * max_size * max_size));
+  cudaCheck(cudaMalloc((void **)&dB, sizeof(float) * max_size * max_size));
+  cudaCheck(cudaMalloc((void **)&dC, sizeof(float) * max_size * max_size));
+  cudaCheck(cudaMalloc((void **)&dC_ref, sizeof(float) * max_size * max_size));
+
+  cudaCheck(cudaMemcpy(dA, A, sizeof(float) * max_size * max_size,
+                       cudaMemcpyHostToDevice));
+  cudaCheck(cudaMemcpy(dB, B, sizeof(float) * max_size * max_size,
+                       cudaMemcpyHostToDevice));
+  cudaCheck(cudaMemcpy(dC, C, sizeof(float) * max_size * max_size,
+                       cudaMemcpyHostToDevice));
+  cudaCheck(cudaMemcpy(dC_ref, C, sizeof(float) * max_size * max_size,
+                       cudaMemcpyHostToDevice));
+
+  int repeat_times = 50;
+  for (int size : SIZE) {
+    m = n = k = size;
+
+    std::cout << "dimensions(m=n=k) " << m << ", alpha: " << alpha
+              << ", beta: " << beta << std::endl;
+    // Verify the correctness of the calculation, and execute it once before the
+    // kernel function timing to avoid cold start errors
+    if (kernel_num != 0) {
+      run_kernel(0, m, n, k, alpha, dA, dB, beta, dC_ref,
+                 handle); // cuBLAS
+      run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC,
+                 handle); // Executes the kernel, modifies the result matrix
+      cudaCheck(cudaDeviceSynchronize());
+      cudaCheck(cudaGetLastError()); // Check for async errors during kernel run
+      cudaMemcpy(C, dC, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
+      cudaMemcpy(C_ref, dC_ref, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
+
+      if (!verify_matrix(C_ref, C, m * n)) {
+        std::cout
+            << "Failed to pass the correctness verification against NVIDIA "
+               "cuBLAS."
+            << std::endl;
+        if (m <= 128) {
+          std::cout << " Logging faulty output into " << errLogFile << "\n";
+          std::ofstream fs;
+          fs.open(errLogFile);
+          fs << "A:\n";
+          print_matrix(A, m, n, fs);
+          fs << "B:\n";
+          print_matrix(B, m, n, fs);
+          fs << "C:\n";
+          print_matrix(C, m, n, fs);
+          fs << "Should:\n";
+          print_matrix(C_ref, m, n, fs);
+        }
+        exit(EXIT_FAILURE);
+      }
+    }
+
+    cudaEventRecord(beg);
+    for (int j = 0; j < repeat_times; j++) {
+      // We don't reset dC between runs to save time
+      run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC, handle);
+    }
+    cudaEventRecord(end);
+    cudaEventSynchronize(beg);
+    cudaEventSynchronize(end);
+    cudaEventElapsedTime(&elapsed_time, beg, end);
+    elapsed_time /= 1000.; // Convert to seconds
+
+    long flops = 2 * m * n * k;
+    printf(
+        "Average elapsed time: (%7.6f) s, performance: (%7.1f) GFLOPS. size: "
+        "(%ld).\n",
+        elapsed_time / repeat_times,
+        (repeat_times * flops * 1e-9) / elapsed_time, m);
+    fflush(stdout);
+    // make dC and dC_ref equal again (we modified dC while calling our kernel
+    // for benchmarking)
+    cudaCheck(cudaMemcpy(dC, dC_ref, sizeof(float) * m * n,
+                         cudaMemcpyDeviceToDevice));
+  }
+
+  // Free up CPU and GPU space
+  free(A);
+  free(B);
+  free(C);
+  free(C_ref);
+  cudaFree(dA);
+  cudaFree(dB);
+  cudaFree(dC);
+  cudaFree(dC_ref);
+  cublasDestroy(handle);
+
+  return 0;
+};

upstream_ref/sgemm_edtallison/simplest_kernel.cu

@@ -0,0 +1,46 @@
+#include <cuda_runtime.h>
+#include <iostream>
+#include <vector>
+
+__global__ void kernel(uint *A, uint *B, int row) {
+  auto x = threadIdx.x / 4;
+  auto y = threadIdx.x % 4;
+  A[x * row + y] = x;
+  B[x * row + y] = y;
+}
+
+int main(int argc, char **argv) {
+  uint *Xs, *Ys;
+  uint *Xs_d, *Ys_d;
+
+  uint SIZE = 4;
+
+  Xs = (uint *)malloc(SIZE * SIZE * sizeof(uint));
+  Ys = (uint *)malloc(SIZE * SIZE * sizeof(uint));
+
+  cudaMalloc((void **)&Xs_d, SIZE * SIZE * sizeof(uint));
+  cudaMalloc((void **)&Ys_d, SIZE * SIZE * sizeof(uint));
+
+  dim3 grid_size(1, 1, 1);
+  dim3 block_size(4 * 4);
+
+  kernel<<<grid_size, block_size>>>(Xs_d, Ys_d, 4);
+
+  cudaMemcpy(Xs, Xs_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
+  cudaMemcpy(Ys, Ys_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
+
+  cudaDeviceSynchronize();
+
+  for (int row = 0; row < SIZE; ++row) {
+    for (int col = 0; col < SIZE; ++col) {
+      std::cout << "[" << Xs[row * SIZE + col] << "|" << Ys[row * SIZE + col]
+                << "] ";
+    }
+    std::cout << "\n";
+  }
+
+  cudaFree(Xs_d);
+  cudaFree(Ys_d);
+  free(Xs);
+  free(Ys);
+}

upstream_ref/sgemm_edtallison/src/kernels.cuh

@@ -0,0 +1,14 @@
+#pragma once
+
+#include "kernels/01_naive.cuh"
+#include "kernels/02_kernel_global_mem_coalesce.cuh"
+#include "kernels/03_kernel_shared_mem_blocking.cuh"
+#include "kernels/04_kernel_1D_blocktiling.cuh"
+#include "kernels/05_kernel_2D_blocktiling.cuh"
+#include "kernels/06_kernel_vectorize.cuh"
+#include "kernels/07_kernel_resolve_bank_conflicts.cuh"
+#include "kernels/08_kernel_bank_extra_col.cuh"
+#include "kernels/09_kernel_autotuned.cuh"
+#include "kernels/10_kernel_warptiling.cuh"
+#include "kernels/11_kernel_double_buffering.cuh"
+#include "kernels/12_kernel_double_buffering.cuh"

upstream_ref/sgemm_edtallison/src/kernels/01_naive.cuh

@@ -0,0 +1,37 @@
+# pragma once
+
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+/*
+
+Matrix sizes:
+MxK * KxN = MxN
+
+*/
+
+__global__ void sgemm_naive(
+    int M, int N, int K, // sizes
+    float alpha, const float *A, const float *B, float beta, float *C // pointers used to point to matrices
+) {
+    // compute position in C that this thread is responsible for
+    // "which block" * "width of block" to get to start of block + "which thread"
+    const uint x = blockIdx.x * blockDim.x + threadIdx.x; // "which row?" (inverted from graphical intuition, confusingly)
+    const uint y = blockIdx.y * blockDim.y + threadIdx.y; // "which column?"
+
+    // if M or N are not multiples of 32, there will be "extra"/"remainder" threads on the last block in x/y.
+    // we don't want those leftover threads to do anything (tile quantisation)
+    if (x < M && y < N) {
+        float tmp = 0.0;
+        for (int i = 0; i < K; ++i) { // K is the size of the row in A, col in B i.e. the dot product
+            // A: x * K gives the start of relevant row, i enumerates across the row (col by col)
+            // B: y gives the relevant column, i * N enumerates down the column, (row by row)
+            tmp += A[x * K + i] * B[i * N + y];
+        }
+        // C = alpha*(A@B) + beta*C
+        // x * N takes to start of relevant row, y moves across to the relevant column
+        C[x * N + y] = alpha * tmp + beta * C[x * N + y];
+    }
+}

upstream_ref/sgemm_edtallison/src/kernels/02_kernel_global_mem_coalesce.cuh

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+template <const uint BLOCKSIZE>
+// __global__ is used to specify that the function is run on GPU, called by host (CPU)
+__global__ void sgemm_global_mem_coalesce(int M, int N, int K, float alpha,
+                                          const float *A, const float *B,
+                                          float beta, float *C) {
+    const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); // note that blockDim is now 1-dimensional
+    const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);
+
+    // if statement is necessary to make things work under tile quantization
+    if (cRow < M && cCol < N) {
+        float tmp = 0.0;
+        for (int i = 0; i < K; ++i) {
+            tmp += A[cRow * K + i] * B[i * N + cCol];
+        }
+        C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol];
+    }
+}

upstream_ref/sgemm_edtallison/src/kernels/03_kernel_shared_mem_blocking.cuh

@@ -0,0 +1,54 @@
+#pragma once
+
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BLOCKSIZE>
+__global__ void sgemm_shared_mem_block(int M, int N, int K, float alpha,
+                                       const float *A, const float *B,
+                                       float beta, float *C) {
+    // output C block we want to compute with this threadBlock
+    const uint cRow = blockIdx.x;
+    const uint cCol = blockIdx.y;
+
+    // allocate buffer for current block in fast SMEM (shared between all threads in block)
+    __shared__ float As[BLOCKSIZE * BLOCKSIZE];
+    __shared__ float Bs[BLOCKSIZE * BLOCKSIZE];
+
+    // the inner row and col that we are accessing in this specific thread
+    const uint threadRow = threadIdx.x / BLOCKSIZE; // note similarity to previous kernel
+    const uint threadCol = threadIdx.x % BLOCKSIZE;
+
+    // advance pointers to the starting positions (they are input as pointers to first elements in the matrices)
+    A += cRow * BLOCKSIZE * K;                    // row=cRow, col=0 (the start of the relevant row)
+    B += cCol * BLOCKSIZE;                        // row=0, col=cCol (top of relevant col)
+    C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE; // row=cRow, col=cCol
+
+    float tmp = 0.0;
+    for (int bkIdx=0; bkIdx < K; bkIdx+=BLOCKSIZE) { // shifting the whole block along the row of A and col of B
+        // have each thread load one of the elements in A and B
+        // make the threadCol (=threadIdx.x) the consecutive index
+        // to allow GMEM access coalescing
+        As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
+        Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
+
+        // ensure cache is fully populated
+        __syncthreads();
+        A += BLOCKSIZE; // for next iteration
+        B += BLOCKSIZE * N;
+
+        // execute the dotproduct on the currently cached block
+        for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx) {
+            tmp += As[threadRow * BLOCKSIZE + dotIdx] *
+                   Bs[dotIdx * BLOCKSIZE + threadCol];
+        }
+        // sync so faster threads don't fetch the next block into cache
+        _syncthreads();
+    }
+    C[threadRow * N + threadCol] = alpha * tmp + beta * C[threadRow * N + threadCol];
+}

upstream_ref/sgemm_edtallison/src/kernels/04_kernel_1D_blocktiling.cuh

@@ -0,0 +1,80 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM>
+__global__ void sgemm1DBlocktiling(int M, int N, int K, float alpha,
+                                   const float *A, const float *B, float beta,
+                                   float *C) {
+  // If we flip x and y here we get ~30% less performance for large matrices.
+  // The current, 30% faster configuration ensures that blocks with sequential
+  // blockIDs access columns of B sequentially, while sharing the same row of A.
+  // The slower configuration would share columns of A, but access into B would
+  // be non-sequential. So the faster configuration has better spatial locality
+  // and hence a greater L2 hit rate.
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // each warp will calculate 32*TM elements, with 32 being the columnar dim.
+  const int threadCol = threadIdx.x % BN;
+  const int threadRow = threadIdx.x / BN;
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // todo: adjust this to each thread to load multiple entries and
+  // better exploit the cache sizes
+  assert(BM * BK == blockDim.x);
+  assert(BN * BK == blockDim.x);
+  const uint innerColA = threadIdx.x % BK; // warp-level GMEM coalescing
+  const uint innerRowA = threadIdx.x / BK;
+  const uint innerColB = threadIdx.x % BN; // warp-level GMEM coalescing
+  const uint innerRowB = threadIdx.x / BN;
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM] = {0.0};
+
+  // outer loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
+    Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;
+    B += BK * N;
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // we make the dotproduct loop the outside loop, which facilitates
+      // reuse of the Bs entry, which we can cache in a tmp var.
+      float tmpB = Bs[dotIdx * BN + threadCol];
+      for (uint resIdx = 0; resIdx < TM; ++resIdx) {
+        threadResults[resIdx] +=
+            As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdx = 0; resIdx < TM; ++resIdx) {
+    C[(threadRow * TM + resIdx) * N + threadCol] =
+        alpha * threadResults[resIdx] +
+        beta * C[(threadRow * TM + resIdx) * N + threadCol];
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/05_kernel_2D_blocktiling.cuh

@@ -0,0 +1,102 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void __launch_bounds__((BM * BN) / (TM * TN), 1)
+    sgemm2DBlocktiling(int M, int N, int K, float alpha, const float *A,
+                       const float *B, float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  const uint totalResultsBlocktile = BM * BN;
+  // A thread is responsible for calculating TM*TN elements in the blocktile
+  const uint numThreadsBlocktile = totalResultsBlocktile / (TM * TN);
+
+  // ResultsPerBlock / ResultsPerThread == ThreadsPerBlock
+  assert(numThreadsBlocktile == blockDim.x);
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  const uint innerRowA = threadIdx.x / BK;
+  const uint innerColA = threadIdx.x % BK;
+  // calculates the number of rows of As that are being loaded in a single step
+  // by a single block
+  const uint strideA = numThreadsBlocktile / BK;
+  const uint innerRowB = threadIdx.x / BN;
+  const uint innerColB = threadIdx.x % BN;
+  // for both As and Bs we want each load to span the full column-width, for
+  // better GMEM coalescing (as opposed to spanning full row-width and iterating
+  // across columns)
+  const uint strideB = numThreadsBlocktile / BN;
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  // register caches for As and Bs
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
+      As[(innerRowA + loadOffset) * BK + innerColA] =
+          A[(innerRowA + loadOffset) * K + innerColA];
+    }
+    for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
+      Bs[(innerRowB + loadOffset) * BN + innerColB] =
+          B[(innerRowB + loadOffset) * N + innerColB];
+    }
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+    for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+      C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN] =
+          alpha * threadResults[resIdxM * TN + resIdxN] +
+          beta * C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN];
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/06_kernel_vectorize.cuh

@@ -0,0 +1,98 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmVectorize(int M, int N, int K, float alpha, float *A,
+                               float *B, float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
+        reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/07_kernel_resolve_bank_conflicts.cuh

@@ -0,0 +1,103 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmResolveBankConflicts(int M, int N, int K, float alpha,
+                                          float *A, float *B, float beta,
+                                          float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    // "linearize" Bs while storing it
+    tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 0) * 16 + innerColB / 2] = tmp.x;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 1) * 16 + innerColB / 2] = tmp.y;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 2) * 16 + innerColB / 2] = tmp.z;
+    Bs[((innerColB % 2) * 4 + innerRowB * 8 + 3) * 16 + innerColB / 2] = tmp.w;
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[(dotIdx * 8 + i) * 16 + threadCol];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/08_kernel_bank_extra_col.cuh

@@ -0,0 +1,103 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void sgemmResolveBankExtraCol(int M, int N, int K, float alpha,
+                                         float *A, float *B, float beta,
+                                         float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // BN/TN are the number of threads to span a column
+  const int threadCol = threadIdx.x % (BN / TN);
+  const int threadRow = threadIdx.x / (BN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  const int extraCols = 5;
+  __shared__ float Bs[BK * (BN + extraCols)];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    // transpose A while loading it
+    float4 tmp =
+        reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
+    As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
+
+    tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 0] = tmp.x;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 1] = tmp.y;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 2] = tmp.z;
+    Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 3] = tmp.w;
+    __syncthreads();
+
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    // calculate per-thread results
+    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+      // block into registers
+      for (uint i = 0; i < TM; ++i) {
+        regM[i] = As[dotIdx * BM + threadRow * TM + i];
+      }
+      for (uint i = 0; i < TN; ++i) {
+        regN[i] = Bs[dotIdx * (BN + extraCols) + threadCol * TN + i];
+      }
+      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+          threadResults[resIdxM * TN + resIdxN] +=
+              regM[resIdxM] * regN[resIdxN];
+        }
+      }
+    }
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+    for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+      // load C vector into registers
+      float4 tmp = reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
+      // perform GEMM update in reg
+      tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
+      tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
+      tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
+      tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
+      // write back
+      reinterpret_cast<float4 *>(
+          &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
+          tmp;
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/09_kernel_autotuned.cuh

@@ -0,0 +1,127 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+const int K9_NUM_THREADS = 256;
+
+template <const int BM, const int BN, const int BK, const int TM, const int TN>
+__global__ void __launch_bounds__(K9_NUM_THREADS)
+    sgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
+                   float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // size of warptile
+  constexpr int WM = TM * 16;
+  constexpr int WN = TN * 16;
+  // iterations of warptile
+  constexpr int WMITER = CEIL_DIV(BM, WM);
+  constexpr int WNITER = CEIL_DIV(BN, WN);
+
+  // Placement of the thread in the warptile
+  const int threadCol = threadIdx.x % (WN / TN);
+  const int threadRow = threadIdx.x / (WN / TN);
+
+  // allocate space for the current blocktile in smem
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  C += cRow * BM * N + cCol * BN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (K9_NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = K9_NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * WNITER * TM * TN] = {0.0};
+  float regM[TM] = {0.0};
+  float regN[TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    // populate the SMEM caches
+    for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+      float4 tmp = reinterpret_cast<float4 *>(
+          &A[(innerRowA + offset) * K + innerColA * 4])[0];
+      // transpose A while storing it
+      As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+      As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+      As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+      As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+    }
+
+    for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+      reinterpret_cast<float4 *>(
+          &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+          reinterpret_cast<float4 *>(
+              &B[(innerRowB + offset) * N + innerColB * 4])[0];
+    }
+    __syncthreads();
+
+    for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
+      for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
+        // calculate per-thread results
+        for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+          // block into registers
+          for (uint i = 0; i < TM; ++i) {
+            regM[i] = As[dotIdx * BM + (wmIdx * WM) + threadRow * TM + i];
+          }
+          for (uint i = 0; i < TN; ++i) {
+            regN[i] = Bs[dotIdx * BN + (wnIdx * WN) + threadCol * TN + i];
+          }
+          for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+            for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+              threadResults[(wmIdx * TM + resIdxM) * (WNITER * TN) +
+                            wnIdx * TN + resIdxN] +=
+                  regM[resIdxM] * regN[resIdxN];
+            }
+          }
+        }
+      }
+    }
+    __syncthreads();
+    // advance blocktile
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+  }
+
+  // write out the results
+  for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
+    for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
+      float *C_interim = C + (wmIdx * WM * N) + (wnIdx * WN);
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRow * TM + resIdxM) * N + threadCol * TN +
+                         resIdxN])[0];
+          // perform GEMM update in reg
+          const int i =
+              (wmIdx * TM + resIdxM) * (WNITER * TN) + wnIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(&C_interim[(threadRow * TM + resIdxM) * N +
+                                                threadCol * TN + resIdxN])[0] =
+              tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/10_kernel_warptiling.cuh

@@ -0,0 +1,187 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+const int WARPSIZE = 32; // warpSize is not constexpr
+
+namespace wt {
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB>
+__device__ void loadFromGmem(int N, int K, const float *A, const float *B,
+                             float *As, float *Bs, int innerRowA, int innerColA,
+                             int innerRowB, int innerColB) {
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    const float4 tmp = reinterpret_cast<const float4 *>(
+        &A[(innerRowA + offset) * K + innerColA * 4])[0];
+    // float4 tmp;
+    // asm("ld.global.nc.v4.f32 {%0, %1, %2, %3}, [%4];"
+    //     : "=f"(tmp.x), "=f"(tmp.y), "=f"(tmp.z), "=f"(tmp.w)
+    //     : "l"(&A[(innerRowA + offset) * K + innerColA * 4]));
+    As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    reinterpret_cast<float4 *>(
+        &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+        reinterpret_cast<const float4 *>(
+            &B[(innerRowB + offset) * N + innerColB * 4])[0];
+    // asm("ld.global.v4.f32 {%0, %1, %2, %3}, [%4];"
+    //     : "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 0]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 1]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 2]),
+    //       "=f"(Bs[(innerRowB + offset) * BN + innerColB * 4 + 3])
+    //     : "l"(&B[(innerRowB + offset) * N + innerColB * 4]));
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace wt
+
+/*
+ * @tparam BM The threadblock size for M dimension SMEM caching.
+ * @tparam BN The threadblock size for N dimension SMEM caching.
+ * @tparam BK The threadblock size for K dimension SMEM caching.
+ * @tparam WM M dim of continuous tile computed by each warp
+ * @tparam WN N dim of continuous tile computed by each warp
+ * @tparam WMITER The number of subwarp tiling steps in M dimension.
+ * @tparam WNITER The number of subwarp tiling steps in N dimension.
+ * @tparam TM The per-thread tile size for M dimension.
+ * @tparam TN The per-thread tile size for N dimension.
+ */
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    sgemmWarptiling(int M, int N, int K, float alpha, float *A, float *B,
+                    float beta, float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[BM * BK];
+  __shared__ float Bs[BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
+    wt::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB);
+    __syncthreads();
+    wt::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                        TN>(regM, regN, threadResults, As, Bs, warpRow, warpCol,
+                            threadRowInWarp, threadColInWarp);
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/11_kernel_double_buffering.cuh

@@ -0,0 +1,220 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+namespace db {
+
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB>
+__device__ void loadFromGmem(const int N, const int K, float *A, float *B,
+                             float *As, float *Bs, const int innerRowA,
+                             const int innerColA, const int innerRowB,
+                             const int innerColB) {
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    float4 tmp = reinterpret_cast<float4 *>(
+        &A[(innerRowA + offset) * K + innerColA * 4])[0];
+    // transpose A while storing it
+    As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
+    As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
+    As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
+    As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    reinterpret_cast<float4 *>(
+        &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
+        reinterpret_cast<float4 *>(
+            &B[(innerRowB + offset) * N + innerColB * 4])[0];
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace db
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    sgemmDoubleBuffering(const int M, const int N, const int K,
+                         const float alpha, float *A, float *B, float beta,
+                         float *C) {
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[2 * BM * BK];
+  __shared__ float Bs[2 * BK * BN];
+
+  // setup double buffering split
+  bool doubleBufferIdx = threadIdx.x >= (NUM_THREADS / 2);
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // for the loading, we're pretending like there's half as many threads
+  // as there actually are
+  const uint innerRowA = (threadIdx.x % (NUM_THREADS / 2)) / (BK / 4);
+  const uint innerColA = (threadIdx.x % (NUM_THREADS / 2)) % (BK / 4);
+  constexpr uint rowStrideA = ((NUM_THREADS / 2) * 4) / BK;
+  const uint innerRowB = (threadIdx.x % (NUM_THREADS / 2)) / (BN / 4);
+  const uint innerColB = (threadIdx.x % (NUM_THREADS / 2)) % (BN / 4);
+  constexpr uint rowStrideB = (NUM_THREADS / 2) / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  if (doubleBufferIdx == 0) {
+    // load first (B0)
+    db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB);
+  }
+  __syncthreads();
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K; bkIdx += 2 * BK) {
+    if (doubleBufferIdx == 0) {
+      // process current (B0)
+      db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                          TN>(regM, regN, threadResults, As, Bs, warpRow,
+                              warpCol, threadRowInWarp, threadColInWarp);
+      __syncthreads();
+
+      // process current+1 (B1)
+      if (bkIdx + BK < K) {
+        db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN,
+                            TM, TN>(regM, regN, threadResults, As + (BM * BK),
+                                    Bs + (BK * BN), warpRow, warpCol,
+                                    threadRowInWarp, threadColInWarp);
+      }
+      __syncthreads();
+
+      // load current + 2 (B0)
+      if (bkIdx + 2 * BK < K) {
+        db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+            N, K, A + 2 * BK, B + 2 * BK * N, As, Bs, innerRowA, innerColA,
+            innerRowB, innerColB);
+      }
+    } else {
+      // load current + 1 (B1)
+      if (bkIdx + BK < K) {
+        db::loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+            N, K, A + BK, B + BK * N, As + (BM * BK), Bs + (BK * BN), innerRowA,
+            innerColA, innerRowB, innerColB);
+      }
+      __syncthreads();
+
+      // process current (B0)
+      db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM,
+                          TN>(regM, regN, threadResults, As, Bs, warpRow,
+                              warpCol, threadRowInWarp, threadColInWarp);
+      __syncthreads();
+
+      // process current+1 (B1)
+      if (bkIdx + BK < K) {
+        db::processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN,
+                            TM, TN>(regM, regN, threadResults, As + (BM * BK),
+                                    Bs + (BK * BN), warpRow, warpCol,
+                                    threadRowInWarp, threadColInWarp);
+      }
+    }
+
+    A += 2 * BK;     // move BK columns to right
+    B += 2 * BK * N; // move BK rows down
+    __syncthreads();
+  }
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/kernels/12_kernel_double_buffering.cuh

@@ -0,0 +1,229 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cooperative_groups.h>
+#include <cstdio>
+#include <cstdlib>
+#include <cublas_v2.h>
+#include <cuda/barrier>
+#include <cuda_runtime.h>
+
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+namespace {
+template <const int BM, const int BN, const int BK, const int rowStrideA,
+          const int rowStrideB, typename T>
+__device__ void loadFromGmem(int N, int K, float *A, float *B, float *As,
+                             float *Bs, int innerRowA, int innerColA,
+                             int innerRowB, int innerColB, T &barrier) {
+
+  for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
+    cuda::memcpy_async(&As[(innerColA * 4 + 0) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 1) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 1],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 2) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 2],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+    cuda::memcpy_async(&As[(innerColA * 4 + 3) * BM + innerRowA + offset],
+                       &A[(innerRowA + offset) * K + innerColA * 4 + 3],
+                       cuda::aligned_size_t<sizeof(float)>(sizeof(float)),
+                       barrier);
+  }
+
+  for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
+    cuda::memcpy_async(&Bs[(innerRowB + offset) * BN + innerColB * 4],
+                       &B[(innerRowB + offset) * N + innerColB * 4],
+                       cuda::aligned_size_t<sizeof(float4)>(sizeof(float4)),
+                       barrier);
+  }
+}
+
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WMITER, const int WNITER, const int WSUBM, const int WSUBN,
+          const int TM, const int TN>
+__device__ void
+processFromSmem(float *regM, float *regN, float *threadResults, const float *As,
+                const float *Bs, const uint warpRow, const uint warpCol,
+                const uint threadRowInWarp, const uint threadColInWarp) {
+  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
+    // populate registers for whole warptile
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint i = 0; i < TM; ++i) {
+        regM[wSubRowIdx * TM + i] =
+            As[(dotIdx * BM) + warpRow * WM + wSubRowIdx * WSUBM +
+               threadRowInWarp * TM + i];
+      }
+    }
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      for (uint i = 0; i < TN; ++i) {
+        regN[wSubColIdx * TN + i] =
+            Bs[(dotIdx * BN) + warpCol * WN + wSubColIdx * WSUBN +
+               threadColInWarp * TN + i];
+      }
+    }
+
+    // execute warptile matmul
+    for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+      for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+        // calculate per-thread results
+        for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
+          for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
+            threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                          (wSubColIdx * TN) + resIdxN] +=
+                regM[wSubRowIdx * TM + resIdxM] *
+                regN[wSubColIdx * TN + resIdxN];
+          }
+        }
+      }
+    }
+  }
+}
+
+} // namespace
+
+/*
+ * @tparam BM The threadblock size for M dimension SMEM caching.
+ * @tparam BN The threadblock size for N dimension SMEM caching.
+ * @tparam BK The threadblock size for K dimension SMEM caching.
+ * @tparam WM M dim of continuous tile computed by each warp
+ * @tparam WN N dim of continuous tile computed by each warp
+ * @tparam WMITER The number of subwarp tiling steps in M dimension.
+ * @tparam WNITER The number of subwarp tiling steps in N dimension.
+ * @tparam TM The per-thread tile size for M dimension.
+ * @tparam TN The per-thread tile size for N dimension.
+ */
+template <const int BM, const int BN, const int BK, const int WM, const int WN,
+          const int WNITER, const int TM, const int TN, const int NUM_THREADS>
+__global__ void __launch_bounds__(NUM_THREADS)
+    runSgemmDoubleBuffering2(int M, int N, int K, float alpha, float *A,
+                             float *B, float beta, float *C) {
+  auto block = cooperative_groups::this_thread_block();
+  __shared__ cuda::barrier<cuda::thread_scope::thread_scope_block> frontBarrier;
+  __shared__ cuda::barrier<cuda::thread_scope::thread_scope_block> backBarrier;
+  auto frontBarrierPtr = &frontBarrier;
+  auto backBarrierPtr = &backBarrier;
+  if (block.thread_rank() == 0) {
+    init(&frontBarrier, block.size());
+    init(&backBarrier, block.size());
+  }
+  __syncthreads();
+
+  const uint cRow = blockIdx.y;
+  const uint cCol = blockIdx.x;
+
+  // Placement of the warp in the threadblock tile
+  const uint warpIdx = threadIdx.x / WARPSIZE; // the warp this thread is in
+  const uint warpCol = warpIdx % (BN / WN);
+  const uint warpRow = warpIdx / (BN / WN);
+
+  // size of the warp subtile
+  constexpr uint WMITER = (WM * WN) / (WARPSIZE * TM * TN * WNITER);
+  constexpr uint WSUBM = WM / WMITER; // 64/2=32
+  constexpr uint WSUBN = WN / WNITER; // 32/2=16
+
+  // Placement of the thread in the warp subtile
+  const uint threadIdxInWarp = threadIdx.x % WARPSIZE;         // [0, 31]
+  const uint threadColInWarp = threadIdxInWarp % (WSUBN / TN); // i%(16/4)
+  const uint threadRowInWarp = threadIdxInWarp / (WSUBN / TN); // i/4
+
+  // allocate space for the current blocktile in SMEM
+  __shared__ float As[2 * BM * BK];
+  __shared__ float Bs[2 * BK * BN];
+
+  // Move blocktile to beginning of A's row and B's column
+  A += cRow * BM * K;
+  B += cCol * BN;
+  // Move C_ptr to warp's output tile
+  C += (cRow * BM + warpRow * WM) * N + cCol * BN + warpCol * WN;
+
+  // calculating the indices that this thread will load into SMEM
+  // we'll load 128bit / 32bit = 4 elements per thread at each step
+  const uint innerRowA = threadIdx.x / (BK / 4);
+  const uint innerColA = threadIdx.x % (BK / 4);
+  constexpr uint rowStrideA = (NUM_THREADS * 4) / BK;
+  const uint innerRowB = threadIdx.x / (BN / 4);
+  const uint innerColB = threadIdx.x % (BN / 4);
+  constexpr uint rowStrideB = NUM_THREADS / (BN / 4);
+
+  // allocate thread-local cache for results in registerfile
+  float threadResults[WMITER * TM * WNITER * TN] = {0.0};
+  // we cache into registers on the warptile level
+  float regM[WMITER * TM] = {0.0};
+  float regN[WNITER * TN] = {0.0};
+
+  int As_offset = 0;
+  int Bs_offset = 0;
+
+  // double-buffering: load first blocktile into SMEM
+  loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+      N, K, A, B, As + As_offset * BM * BK, Bs + Bs_offset * BK * BN, innerRowA,
+      innerColA, innerRowB, innerColB, (*frontBarrierPtr));
+
+  // outer-most loop over block tiles
+  for (uint bkIdx = 0; bkIdx < K - BK; bkIdx += BK) {
+    // double-buffering: load next blocktile into SMEM
+    loadFromGmem<BM, BN, BK, rowStrideA, rowStrideB>(
+        N, K, A + BK, B + BK * N, As + (1 - As_offset) * BM * BK,
+        Bs + (1 - Bs_offset) * BK * BN, innerRowA, innerColA, innerRowB,
+        innerColB, (*backBarrierPtr));
+
+    // compute the current blocktile
+    (*frontBarrierPtr).arrive_and_wait();
+    processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM, TN>(
+        regM, regN, threadResults, As + As_offset * BM * BK,
+        Bs + Bs_offset * BK * BN, warpRow, warpCol, threadRowInWarp,
+        threadColInWarp);
+    A += BK;     // move BK columns to right
+    B += BK * N; // move BK rows down
+
+    As_offset = 1 - As_offset;
+    Bs_offset = 1 - Bs_offset;
+    // swap the front and back barriers
+    auto tmp = frontBarrierPtr;
+    frontBarrierPtr = backBarrierPtr;
+    backBarrierPtr = tmp;
+
+    __syncthreads();
+  }
+
+  // compute the last blocktile
+  (*frontBarrierPtr).arrive_and_wait();
+  processFromSmem<BM, BN, BK, WM, WN, WMITER, WNITER, WSUBM, WSUBN, TM, TN>(
+      regM, regN, threadResults, As + As_offset * BM * BK,
+      Bs + Bs_offset * BK * BN, warpRow, warpCol, threadRowInWarp,
+      threadColInWarp);
+
+  // write out the results
+  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
+    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
+      // move C pointer to current warp subtile
+      float *C_interim = C + (wSubRowIdx * WSUBM) * N + wSubColIdx * WSUBN;
+      for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
+        for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
+          // load C vector into registers
+          float4 tmp = reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0];
+          // perform GEMM update in reg
+          const int i = (wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
+                        wSubColIdx * TN + resIdxN;
+          tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
+          tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
+          tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
+          tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
+          // write back
+          reinterpret_cast<float4 *>(
+              &C_interim[(threadRowInWarp * TM + resIdxM) * N +
+                         threadColInWarp * TN + resIdxN])[0] = tmp;
+        }
+      }
+    }
+  }
+}

upstream_ref/sgemm_edtallison/src/runner.cu

@@ -0,0 +1,549 @@
+#include "kernels.cuh"
+#include "runner.cuh"
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <iomanip>
+
+float get_sec() {
+  struct timeval time;
+  gettimeofday(&time, NULL);
+  return (1e6 * time.tv_sec + time.tv_usec);
+}
+
+float cpu_elapsed_time(float &beg, float &end) { return 1.0e-6 * (end - beg); }
+
+void cudaCheck(cudaError_t error, const char *file, int line) {
+  if (error != cudaSuccess) {
+    printf("[CUDA ERROR] at file %s:%d:\n%s\n", file, line,
+           cudaGetErrorString(error));
+    exit(EXIT_FAILURE);
+  }
+};
+
+void CudaDeviceInfo() {
+  int deviceId;
+
+  cudaGetDevice(&deviceId);
+
+  cudaDeviceProp props{};
+  cudaGetDeviceProperties(&props, deviceId);
+
+  printf("Device ID: %d\n\
+    Name: %s\n\
+    Compute Capability: %d.%d\n\
+    memoryBusWidth: %d\n\
+    maxThreadsPerBlock: %d\n\
+    maxThreadsPerMultiProcessor: %d\n\
+    maxRegsPerBlock: %d\n\
+    maxRegsPerMultiProcessor: %d\n\
+    totalGlobalMem: %zuMB\n\
+    sharedMemPerBlock: %zuKB\n\
+    sharedMemPerMultiprocessor: %zuKB\n\
+    totalConstMem: %zuKB\n\
+    multiProcessorCount: %d\n\
+    Warp Size: %d\n",
+         deviceId, props.name, props.major, props.minor, props.memoryBusWidth,
+         props.maxThreadsPerBlock, props.maxThreadsPerMultiProcessor,
+         props.regsPerBlock, props.regsPerMultiprocessor,
+         props.totalGlobalMem / 1024 / 1024, props.sharedMemPerBlock / 1024,
+         props.sharedMemPerMultiprocessor / 1024, props.totalConstMem / 1024,
+         props.multiProcessorCount, props.warpSize);
+};
+
+void randomize_matrix(float *mat, int N) {
+  // NOTICE: Use gettimeofday instead of srand((unsigned)time(NULL)); the time
+  // precision is too low and the same random number is generated.
+  struct timeval time {};
+  gettimeofday(&time, nullptr);
+  srand(time.tv_usec);
+  for (int i = 0; i < N; i++) {
+    float tmp = (float)(rand() % 5) + 0.01 * (rand() % 5);
+    tmp = (rand() % 2 == 0) ? tmp : tmp * (-1.);
+    mat[i] = tmp;
+  }
+}
+
+void range_init_matrix(float *mat, int N) {
+  for (int i = 0; i < N; i++) {
+    mat[i] = i;
+  }
+}
+
+void zero_init_matrix(float *mat, int N) {
+  for (int i = 0; i < N; i++) {
+    mat[i] = 0.0;
+  }
+}
+
+void copy_matrix(const float *src, float *dest, int N) {
+  int i;
+  for (i = 0; src + i && dest + i && i < N; i++)
+    *(dest + i) = *(src + i);
+  if (i != N)
+    printf("copy failed at %d while there are %d elements in total.\n", i, N);
+}
+
+void print_matrix(const float *A, int M, int N, std::ofstream &fs) {
+  int i;
+  fs << std::setprecision(2)
+     << std::fixed; // Set floating-point precision and fixed notation
+  fs << "[";
+  for (i = 0; i < M * N; i++) {
+    if ((i + 1) % N == 0)
+      fs << std::setw(5) << A[i]; // Set field width and write the value
+    else
+      fs << std::setw(5) << A[i] << ", ";
+    if ((i + 1) % N == 0) {
+      if (i + 1 < M * N)
+        fs << ";\n";
+    }
+  }
+  fs << "]\n";
+}
+
+bool verify_matrix(float *matRef, float *matOut, int N) {
+  double diff = 0.0;
+  int i;
+  for (i = 0; i < N; i++) {
+    diff = std::fabs(matRef[i] - matOut[i]);
+    if (diff > 0.01) {
+      printf("Divergence! Should %5.2f, Is %5.2f (Diff %5.2f) at %d\n",
+             matRef[i], matOut[i], diff, i);
+      return false;
+    }
+  }
+  return true;
+}
+
+int div_ceil(int numerator, int denominator) {
+  std::div_t res = std::div(numerator, denominator);
+  return res.rem ? (res.quot + 1) : res.quot;
+}
+
+void runCublasFP32(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // cuBLAS uses column-major order. So we change the order of our row-major A &
+  // B, since (B^T*A^T)^T = (A*B)
+  // This runs cuBLAS in full fp32 mode
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N, CUBLAS_COMPUTE_32F,
+               CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void runCublasBF16(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // This runs cuBLAS with mixed precision (performing the mul with operands
+  // downcast to bf16), which is ~4x faster
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N,
+               CUBLAS_COMPUTE_32F_FAST_16BF, CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void runCublasTF32(cublasHandle_t handle, int M, int N, int K, float alpha,
+                   float *A, float *B, float beta, float *C) {
+  // This runs cuBLAS with mixed precision (performing the mul with operands
+  // downcast to bf16), which is ~4x faster
+  cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, CUDA_R_32F,
+               N, A, CUDA_R_32F, K, &beta, C, CUDA_R_32F, N,
+               CUBLAS_COMPUTE_32F_FAST_TF32, CUBLAS_GEMM_DEFAULT_TENSOR_OP);
+}
+
+void run_sgemm_naive(int M, int N, int K, float alpha, float *A, float *B,
+                     float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32, 32);
+  sgemm_naive<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_sgemm_coalesce(int M, int N, int K, float alpha, float *A, float *B,
+                        float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32 * 32);
+  sgemm_global_mem_coalesce<32>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_sgemm_shared_mem_block(int M, int N, int K, float alpha, float *A,
+                                float *B, float beta, float *C) {
+  dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+  dim3 blockDim(32 * 32);
+  // L1 cache becomes useless, since we access GMEM only via SMEM, so we carve
+  // out all of L1 to SMEM. This doesn't currently make a difference, since
+  // occupancy is limited by reg and thread count, but it's good to do anyway.
+  cudaFuncSetAttribute(sgemm_shared_mem_block<32>,
+                       cudaFuncAttributePreferredSharedMemoryCarveout,
+                       cudaSharedmemCarveoutMaxShared);
+  sgemm_shared_mem_block<32>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemm1DBlocktiling(int M, int N, int K, float alpha, float *A, float *B,
+                           float beta, float *C) {
+  const uint BM = 64;
+  const uint BN = 64;
+  const uint BK = 8;
+  const uint TM = 8;
+  dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+  dim3 blockDim((BM * BN) / TM);
+  sgemm1DBlocktiling<BM, BN, BK, TM>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemm2DBlocktiling(int M, int N, int K, float alpha, float *A, float *B,
+                           float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemm2DBlocktiling<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemm2DBlocktiling<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmVectorize(int M, int N, int K, float alpha, float *A, float *B,
+                       float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmVectorize<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmVectorize<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmResolveBankConflicts(int M, int N, int K, float alpha, float *A,
+                                  float *B, float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankConflicts<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankConflicts<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmResolveBankExtraCol(int M, int N, int K, float alpha, float *A,
+                                 float *B, float beta, float *C) {
+  const uint BK = 8;
+  const uint TM = 8;
+  const uint TN = 8;
+  if (M >= 128 and N >= 128) {
+    const uint BM = 128;
+    const uint BN = 128;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankExtraCol<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  } else {
+    // this is a hacky solution to the underlying problem
+    // of not having proper bounds checking in the kernel
+    const uint BM = 64;
+    const uint BN = 64;
+    dim3 gridDim(CEIL_DIV(N, BN), CEIL_DIV(M, BM));
+    dim3 blockDim((BM * BN) / (TM * TN));
+    sgemmResolveBankExtraCol<BM, BN, BK, TM, TN>
+        <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+  }
+}
+
+void runSgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
+                       float beta, float *C) {
+  // A100
+  // const uint K9_BK = 16;
+  // const uint K9_TM = 4;
+  // const uint K9_TN = 4;
+  // const uint K9_BM = 64;
+  // const uint K9_BN = 64;
+  // A6000
+  const uint K9_BK = 16;
+  const uint K9_TM = 8;
+  const uint K9_TN = 8;
+  const uint K9_BM = 128;
+  const uint K9_BN = 128;
+  dim3 blockDim(K9_NUM_THREADS);
+
+  static_assert(
+      (K9_NUM_THREADS * 4) % K9_BK == 0,
+      "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization issues "
+      "during GMEM->SMEM tiling (loading only parts of the final row of Bs "
+      "during each iteraion)");
+  static_assert(
+      (K9_NUM_THREADS * 4) % K9_BN == 0,
+      "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization issues "
+      "during GMEM->SMEM tiling (loading only parts of the final row of As "
+      "during each iteration)");
+  static_assert(
+      K9_BN % (16 * K9_TN) == 0,
+      "K9_BN must be a multiple of 16*K9_TN to avoid quantization effects");
+  static_assert(
+      K9_BM % (16 * K9_TM) == 0,
+      "K9_BM must be a multiple of 16*K9_TM to avoid quantization effects");
+  static_assert((K9_BM * K9_BK) % (4 * K9_NUM_THREADS) == 0,
+                "K9_BM*K9_BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K9_BN * K9_BK) % (4 * K9_NUM_THREADS) == 0,
+                "K9_BN*K9_BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K9_BN), CEIL_DIV(M, K9_BM));
+  sgemmAutotuned<K9_BM, K9_BN, K9_BK, K9_TM, K9_TN>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmWarptiling(int M, int N, int K, float alpha, float *A, float *B,
+                        float beta, float *C) {
+  // Settings for A100
+  // const uint K10_NUM_THREADS = 128;
+  // const uint K10_BN = 128;
+  // const uint K10_BM = 64;
+  // const uint K10_BK = 16;
+  // const uint K10_WN = 64;
+  // const uint K10_WM = 32;
+  // const uint K10_WNITER = 1;
+  // const uint K10_TN = 4;
+  // const uint K10_TM = 4;
+  // Settings for A6000
+  const uint K10_NUM_THREADS = 128;
+  const uint K10_BN = 128;
+  const uint K10_BM = 128;
+  const uint K10_BK = 16;
+  const uint K10_WN = 64;
+  const uint K10_WM = 64;
+  const uint K10_WNITER = 4;
+  const uint K10_TN = 4;
+  const uint K10_TM = 8;
+  dim3 blockDim(K10_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K10_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K10_BN % K10_WN == 0) and (K10_BM % K10_WM == 0));
+  static_assert((K10_BN / K10_WN) * (K10_BM / K10_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K10_WM * K10_WN) % (WARPSIZE * K10_TM * K10_TN * K10_WNITER) ==
+                0);
+  constexpr uint K10_WMITER =
+      (K10_WM * K10_WN) / (32 * K10_TM * K10_TN * K10_WNITER);
+  // warpsubtile in warptile
+  static_assert((K10_WM % K10_WMITER == 0) and (K10_WN % K10_WNITER == 0));
+
+  static_assert((K10_NUM_THREADS * 4) % K10_BK == 0,
+                "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K10_NUM_THREADS * 4) % K10_BN == 0,
+                "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K10_BN % (16 * K10_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K10_BM % (16 * K10_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K10_BM * K10_BK) % (4 * K10_NUM_THREADS) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K10_BN * K10_BK) % (4 * K10_NUM_THREADS) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K10_BN), CEIL_DIV(M, K10_BM));
+  sgemmWarptiling<K10_BM, K10_BN, K10_BK, K10_WM, K10_WN, K10_WNITER, K10_TM,
+                  K10_TN, K10_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmDoubleBuffering(int M, int N, int K, float alpha, float *A,
+                             float *B, float beta, float *C) {
+  // Settings for A100
+  // const uint K11_NUM_THREADS = 256;
+  // const uint K11_BN = 128;
+  // const uint K11_BM = 64;
+  // const uint K11_BK = 16;
+  // const uint K11_WN = 32;
+  // const uint K11_WM = 32;
+  // const uint K11_WNITER = 2;
+  // const uint K11_TN = 4;
+  // const uint K11_TM = 4;
+  // Settings for A6000
+  const uint K11_NUM_THREADS = 256;
+  const uint K11_BN = 256;
+  const uint K11_BM = 128;
+  const uint K11_BK = 16;
+  const uint K11_WN = 32;
+  const uint K11_WM = 128;
+  const uint K11_WNITER = 1;
+  const uint K11_TN = 8;
+  const uint K11_TM = 8;
+  dim3 blockDim(K11_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K11_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K11_BN % K11_WN == 0) and (K11_BM % K11_WM == 0));
+  static_assert((K11_BN / K11_WN) * (K11_BM / K11_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K11_WM * K11_WN) % (WARPSIZE * K11_TM * K11_TN * K11_WNITER) ==
+                0);
+  constexpr uint K11_WMITER =
+      (K11_WM * K11_WN) / (32 * K11_TM * K11_TN * K11_WNITER);
+  // warpsubtile in warptile
+  static_assert((K11_WM % K11_WMITER == 0) and (K11_WN % K11_WNITER == 0));
+
+  static_assert((K11_NUM_THREADS / 2 * 4) % K11_BK == 0,
+                "NUM_THREADS*4 must be multiple of BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K11_NUM_THREADS / 2 * 4) % K11_BN == 0,
+                "NUM_THREADS*4 must be multiple of BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K11_BN % (16 * K11_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K11_BM % (16 * K11_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K11_BM * K11_BK) % (4 * K11_NUM_THREADS / 2) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K11_BN * K11_BK) % (4 * K11_NUM_THREADS / 2) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K11_BN), CEIL_DIV(M, K11_BM));
+  sgemmDoubleBuffering<K11_BM, K11_BN, K11_BK, K11_WM, K11_WN, K11_WNITER,
+                       K11_TM, K11_TN, K11_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void runSgemmDoubleBuffering2(int M, int N, int K, float alpha, float *A,
+                              float *B, float beta, float *C) {
+  // Settings for A6000
+  const uint K12_NUM_THREADS = 128;
+  const uint K12_BN = 128;
+  const uint K12_BM = 128;
+  const uint K12_BK = 16;
+  const uint K12_WN = 64;
+  const uint K12_WM = 64;
+  const uint K12_WNITER = 4;
+  const uint K12_TN = 4;
+  const uint K12_TM = 8;
+  dim3 blockDim(K12_NUM_THREADS);
+
+  constexpr uint NUM_WARPS = K12_NUM_THREADS / 32;
+
+  // warptile in threadblocktile
+  static_assert((K12_BN % K12_WN == 0) and (K12_BM % K12_WM == 0));
+  static_assert((K12_BN / K12_WN) * (K12_BM / K12_WM) == NUM_WARPS);
+
+  // threads in warpsubtile
+  static_assert((K12_WM * K12_WN) % (WARPSIZE * K12_TM * K12_TN * K12_WNITER) ==
+                0);
+  constexpr uint K12_WMITER =
+      (K12_WM * K12_WN) / (32 * K12_TM * K12_TN * K12_WNITER);
+  // warpsubtile in warptile
+  static_assert((K12_WM % K12_WMITER == 0) and (K12_WN % K12_WNITER == 0));
+
+  static_assert((K12_NUM_THREADS * 4) % K12_BK == 0,
+                "NUM_THREADS*4 must be multiple of K9_BK to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of Bs during each iteraion)");
+  static_assert((K12_NUM_THREADS * 4) % K12_BN == 0,
+                "NUM_THREADS*4 must be multiple of K9_BN to avoid quantization "
+                "issues during GMEM->SMEM tiling (loading only parts of the "
+                "final row of As during each iteration)");
+  static_assert(K12_BN % (16 * K12_TN) == 0,
+                "BN must be a multiple of 16*TN to avoid quantization effects");
+  static_assert(K12_BM % (16 * K12_TM) == 0,
+                "BM must be a multiple of 16*TM to avoid quantization effects");
+  static_assert((K12_BM * K12_BK) % (4 * K12_NUM_THREADS) == 0,
+                "BM*BK must be a multiple of 4*256 to vectorize loads");
+  static_assert((K12_BN * K12_BK) % (4 * K12_NUM_THREADS) == 0,
+                "BN*BK must be a multiple of 4*256 to vectorize loads");
+
+  dim3 gridDim(CEIL_DIV(N, K12_BN), CEIL_DIV(M, K12_BM));
+  runSgemmDoubleBuffering2<K12_BM, K12_BN, K12_BK, K12_WM, K12_WN, K12_WNITER,
+                           K12_TM, K12_TN, K12_NUM_THREADS>
+      <<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void run_kernel(int kernel_num, int M, int N, int K, float alpha, float *A,
+                float *B, float beta, float *C, cublasHandle_t handle) {
+  switch (kernel_num) {
+  case 0:
+    runCublasFP32(handle, M, N, K, alpha, A, B, beta, C);
+    break;
+  case 1:
+    run_sgemm_naive(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 2:
+    run_sgemm_coalesce(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 3:
+    run_sgemm_shared_mem_block(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 4:
+    runSgemm1DBlocktiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 5:
+    runSgemm2DBlocktiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 6:
+    runSgemmVectorize(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 7:
+    runSgemmResolveBankConflicts(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 8:
+    runSgemmResolveBankExtraCol(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 9:
+    runSgemmAutotuned(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 10:
+    runSgemmWarptiling(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 11:
+    runSgemmDoubleBuffering(M, N, K, alpha, A, B, beta, C);
+    break;
+  case 12:
+    runSgemmDoubleBuffering2(M, N, K, alpha, A, B, beta, C);
+    break;
+  default:
+    throw std::invalid_argument("Unknown kernel number");
+  }
+}

upstream_ref/sgemm_edtallison/src/runner.cuh

@@ -0,0 +1,26 @@
+#pragma once
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+#include <fstream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <time.h>
+#include <unistd.h>
+
+void cudaCheck(cudaError_t error, const char *file,
+               int line); // CUDA error check
+void CudaDeviceInfo();    // print CUDA information
+
+void range_init_matrix(float *mat, int N);
+void randomize_matrix(float *mat, int N);
+void zero_init_matrix(float *mat, int N);
+void copy_matrix(const float *src, float *dest, int N);
+void print_matrix(const float *A, int M, int N, std::ofstream &fs);
+bool verify_matrix(float *mat1, float *mat2, int N);
+
+float get_current_sec();                        // Get the current moment
+float cpu_elapsed_time(float &beg, float &end); // Calculate time difference
+
+void run_kernel(int kernel_num, int m, int n, int k, float alpha, float *A,
+                float *B, float beta, float *C, cublasHandle_t handle);