data: complete SGEMM upstream from 3 repos (siboehm+wangzyon+edtallison) + xllm fused_qknorm_rope + xattention kernels

SGEMM repos (upstream_ref/sgemm_cuda/, 41 files): siboehm/SGEMM_CUDA: kernel 1-12, runner, CMake, cuBLAS benchmark wangzyon/NVIDIA_SGEMM_PRACTICE: kernel 1-7 (Chinese comments), utils edtallison/sgemm-cuda: kernel 01-09 (learning notes), Makefile xllm kernels (ex_engine/xllm_kernels/cuda/): fused_qknorm_rope.cu + bind — saves 128 kernel launches/fwd xattention/ — 6 files from upstream xllm headers: corex_compat_utils.h, topk_last_dim.cuh ilu/CMakeLists.txt SO_BUILD_MANIFEST.md — complete .so inventory and call chain analysis
2026-08-15 07:00:04 +00:00
parent 7cfa87b5ac
commit 36676f2d1b
42 changed files with 6808 additions and 0 deletions
--- a/upstream_ref/sgemm_cuda/01_naive.cuh
+++ b/upstream_ref/sgemm_cuda/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];
+    }
+}
--- a/upstream_ref/sgemm_cuda/02_kernel_global_mem_coalesce.cuh
+++ b/upstream_ref/sgemm_cuda/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];
+    }
+}
--- a/upstream_ref/sgemm_cuda/03_kernel_shared_mem_blocking.cuh
+++ b/upstream_ref/sgemm_cuda/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];
+}
--- a/upstream_ref/sgemm_cuda/04_kernel_1D_blocktiling.cuh
+++ b/upstream_ref/sgemm_cuda/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];
+  }
+}
--- a/upstream_ref/sgemm_cuda/05_kernel_2D_blocktiling.cuh
+++ b/upstream_ref/sgemm_cuda/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];
+    }
+  }
+}
--- a/upstream_ref/sgemm_cuda/06_kernel_vectorize.cuh
+++ b/upstream_ref/sgemm_cuda/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;
+    }
+  }
+}
--- a/upstream_ref/sgemm_cuda/07_kernel_resolve_bank_conflicts.cuh
+++ b/upstream_ref/sgemm_cuda/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;
+    }
+  }
+}
--- a/upstream_ref/sgemm_cuda/08_kernel_bank_extra_col.cuh
+++ b/upstream_ref/sgemm_cuda/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;
+    }
+  }
+}
--- a/upstream_ref/sgemm_cuda/09_kernel_autotuned.cuh
+++ b/upstream_ref/sgemm_cuda/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;
+        }
+      }
+    }
+  }
+}
--- a/upstream_ref/sgemm_cuda/CMakeLists_edtallison.txt
+++ b/upstream_ref/sgemm_cuda/CMakeLists_edtallison.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})
--- a/upstream_ref/sgemm_cuda/CMakeLists_siboehm.txt
+++ b/upstream_ref/sgemm_cuda/CMakeLists_siboehm.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 86)
+
+# 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})
--- a/upstream_ref/sgemm_cuda/Makefile_edtallison
+++ b/upstream_ref/sgemm_cuda/Makefile_edtallison
@@ -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
--- a/upstream_ref/sgemm_cuda/cuBLAS_sgemm.cu
+++ b/upstream_ref/sgemm_cuda/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;
+}
--- a/upstream_ref/sgemm_cuda/kernel_1.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_1.cuh
@@ -0,0 +1,19 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+__global__ __launch_bounds__(1024) void
+mysgemm_v1(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+
+    int gx = blockIdx.x * blockDim.x + threadIdx.x; // 全局x
+    int gy = blockIdx.y * blockDim.y + threadIdx.y; // 全局y
+
+    float tmp = 0.;
+    for (int i = 0; i < K; i++) {
+        tmp += A[gy * K + i] * B[i * N + gx]; // 两次全局内存访问和一次FMA（累加乘）
+    }
+    C[gy * N + gx] = alpha * tmp + beta * C[gy * N + gx];
+}
--- a/upstream_ref/sgemm_cuda/kernel_2.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_2.cuh
@@ -0,0 +1,45 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+template<const int BLOCK_SIZE>
+__global__ void mysgemm_v2(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+
+    const int BM = BLOCK_SIZE;
+    const int BN = BLOCK_SIZE;
+    const int BK = BLOCK_SIZE;
+
+    int tx = threadIdx.x % BN;
+    int ty = threadIdx.x / BN;
+
+    // 申请共享内存空间
+    __shared__ float As[BM * BK];
+    __shared__ float Bs[BK * BN];
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+    float tmp = 0.;
+    for (int k = 0; k < K; k += BK) {
+        // 缓存A_tile和B_tile
+        As[ty * BK + tx] = A[ty * K + tx];
+        Bs[ty * BN + tx] = B[ty * N + tx];
+        // 同步所有线程缓存完成
+        __syncthreads();
+        A += BK;
+        B += BK * N;
+        for (int i = 0; i < BK; i++) {
+            tmp += As[ty * BK + i] * Bs[i * BN + tx];
+        }
+        // FMA计算需要读取缓存数据，在新一轮写入缓存前进行同步，确保所有线程计算完成
+        __syncthreads();
+    }
+    C[ty * N + tx] = alpha * tmp + beta * C[ty * N + tx];
+}
--- a/upstream_ref/sgemm_cuda/kernel_3.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_3.cuh
@@ -0,0 +1,71 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+template<const int BM,
+        const int BN,
+        const int BK,
+        const int TM>
+__global__ void mysgemm_v3(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+    int thread_num = BM * BN / TM; // 一个线程负责block中计算TM个元素
+
+    int tx = threadIdx.x % BN;
+    int ty = threadIdx.x / BN * TM;
+
+    __shared__ float As[BM * BK];
+    __shared__ float Bs[BK * BN];
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+    /*
+    当前线程负责搬运全局内存中第a_tile_row行，第a_tile_col列元素至共享内存第a_tile_row行，第a_tile_col列
+    a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存；
+
+    若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM，表示每个线程搬运一轮即可完成所需元素的搬运;
+    若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
+    */
+    int a_tile_row = threadIdx.x / BK;
+    int a_tile_col = threadIdx.x % BK;
+    int a_tile_stride = thread_num / BK;
+
+    int b_tile_row = threadIdx.x / BN;
+    int b_tile_col = threadIdx.x % BN;
+    int b_tile_stride = thread_num / BN;
+
+    float tmp[TM + 1] = {0.}; // 每个线程负责TM个元素，则需要申请TM个寄存器保存累加值，额外的一个寄存器用于缓存；
+    #pragma unroll
+    for (int k = 0; k < K; k += BK) {
+        #pragma unroll
+        for (int i = 0; i < BM; i += a_tile_stride) {
+            As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
+        }
+        #pragma unroll
+        for (int i = 0; i < BK; i += b_tile_stride) {
+            Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
+        }
+        __syncthreads();
+        A += BK;
+        B += BK * N;
+        #pragma unroll
+        for (int i = 0; i < BK; i++) {
+            tmp[TM] = Bs[tx + i * BN]; // 额外的一个寄存器，避免反复从共享内存中读取Bs[tx + i * BN]
+            #pragma unroll  // 循环展开，增加指令并行度
+            for (int j = 0; j < TM; j++) {
+                tmp[j] += As[(ty + j) * BK + i] * tmp[TM];
+            }
+        }
+        __syncthreads();
+    }
+    #pragma unroll
+    for (int j = 0; j < TM; j++) {
+        C[(ty + j) * N + tx] = alpha * tmp[j] + beta * C[(ty + j) * N + tx];
+    }
+}
--- a/upstream_ref/sgemm_cuda/kernel_4.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_4.cuh
@@ -0,0 +1,76 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+template<const int BM,
+        const int BN,
+        const int BK,
+        const int TM,
+        const int TN>
+__global__ void mysgemm_v4(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+
+    int block_row_thread = BN / TN;
+    int block_col_thread = BM / TM;
+    int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
+
+    int tx = (threadIdx.x % block_row_thread) * TN;
+    int ty = (threadIdx.x / block_row_thread) * TM;
+
+    __shared__ float As[BM * BK];
+    __shared__ float Bs[BK * BN];
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+    /*
+    当前线程负责搬运全局内存中第a_tile_row行，第a_tile_col列元素至共享内存第a_tile_row行，第a_tile_col列
+    a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存；
+
+    若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM，表示每个线程搬运一轮即可完成所需元素的搬运;
+    若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
+    */
+    int a_tile_row = threadIdx.x / BK;
+    int a_tile_col = threadIdx.x % BK;
+    int a_tile_stride = thread_num / BK;
+
+    int b_tile_row = threadIdx.x / BN;
+    int b_tile_col = threadIdx.x % BN;
+    int b_tile_stride = thread_num / BN;
+
+    float tmp[TM][TN] = {0.}; // 每个线程负责TM*TN个元素，则需要申请TM*TN个寄存器保存累加值，额外的一个寄存器用于缓存；
+    #pragma unroll
+    for (int k = 0; k < K; k += BK) {
+        #pragma unroll
+        for (int i = 0; i < BM; i += a_tile_stride) {
+            As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
+        }
+        #pragma unroll
+        for (int i = 0; i < BK; i += b_tile_stride) {
+            Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
+        }
+        __syncthreads();
+        A += BK;
+        B += BK * N;
+        #pragma unroll
+        for (int i = 0; i < BK; i++) {
+            #pragma unroll  // 循环展开，增加指令并行度
+            for (int j = 0; j < TM; j++) {
+                for (int l = 0; l < TN; l++)
+                    tmp[j][l] += As[(ty + j) * BK + i] * Bs[tx + l + i * BN];
+            }
+        }
+        __syncthreads();
+    }
+    #pragma unroll
+    for (int j = 0; j < TM; j++) {
+        for (int l = 0; l < TN; l++)
+            C[(ty + j) * N + tx + l] = alpha * tmp[j][l] + beta * C[(ty + j) * N + tx + l];
+    }
+}
--- a/upstream_ref/sgemm_cuda/kernel_5.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_5.cuh
@@ -0,0 +1,88 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+template<const int BM,
+        const int BN,
+        const int BK,
+        const int TM,
+        const int TN>
+__global__ void mysgemm_v5(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+
+    int block_row_thread = BN / TN;
+    int block_col_thread = BM / TM;
+    int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
+
+    int tx = (threadIdx.x % block_row_thread) * TN;
+    int ty = (threadIdx.x / block_row_thread) * TM;
+
+    __shared__ float As[BM * BK];
+    __shared__ float Bs[BK * BN];
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+    /*
+    当前线程负责搬运全局内存中第a_tile_row行，第a_tile_col列元素至共享内存第a_tile_row行，第a_tile_col列
+    a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存；
+
+    若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM，表示每个线程搬运一轮即可完成所需元素的搬运;
+    若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
+    */
+    int a_tile_row = threadIdx.x / BK;
+    int a_tile_col = threadIdx.x % BK;
+    int a_tile_stride = thread_num / BK;
+
+    int b_tile_row = threadIdx.x / BN;
+    int b_tile_col = threadIdx.x % BN;
+    int b_tile_stride = thread_num / BN;
+
+    float tmp[TM][TN] = {0.}; // 每个线程负责TM*TN个元素，则需要申请TM*TN个寄存器保存累加值，额外的一个寄存器用于缓存；
+    float a_frag[TM] = {0.};
+    float b_frag[TN] = {0.};
+
+    #pragma unroll
+    for (int k = 0; k < K; k += BK) {
+        #pragma unroll
+        for (int i = 0; i < BM; i += a_tile_stride) {
+            As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
+        }
+        #pragma unroll
+        for (int i = 0; i < BK; i += b_tile_stride) {
+            Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
+        }
+        __syncthreads();
+        A += BK;
+        B += BK * N;
+        #pragma unroll
+        for (int i = 0; i < BK; i++) {
+            #pragma unroll
+            for (int j = 0; j < TM; j++) {
+                a_frag[j] = As[(ty + j) * BK + i];
+            }
+            #pragma unroll
+            for (int l = 0; l < TN; l++) {
+                b_frag[l] = Bs[tx + l + i * BN];
+            }
+            #pragma unroll
+            for (int j = 0; j < TM; j++) {
+                #pragma unroll
+                for (int l = 0; l < TN; l++)
+                    tmp[j][l] += a_frag[j] * b_frag[l];
+            }
+        }
+        __syncthreads();
+    }
+    #pragma unroll
+    for (int j = 0; j < TM; j++) {
+        for (int l = 0; l < TN; l++)
+            C[(ty + j) * N + tx + l] = alpha * tmp[j][l] + beta * C[(ty + j) * N + tx + l];
+    }
+}
--- a/upstream_ref/sgemm_cuda/kernel_6.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_6.cuh
@@ -0,0 +1,110 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define OFFSET(row, col, ld) ((row)*(ld)+(col))
+#define FETCH_FLOAT4(pointer) (reinterpret_cast<float4*>(&(pointer))[0])
+
+template<const int BM,
+        const int BN,
+        const int BK,
+        const int TM,
+        const int TN>
+__global__ void mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+
+    const int block_row_thread = BN / TN;
+    const int block_col_thread = BM / TM;
+    const int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
+
+    // 当前线程对应thread tile的左上角元素在block中的位置
+    int tx = (threadIdx.x % block_row_thread) * TN;
+    int ty = (threadIdx.x / block_row_thread) * TM;
+
+    __shared__ float As[BK * BM];
+    __shared__ float Bs[BK * BN];
+
+
+    const int ldg_a_num = BK * BM / thread_num / 4; // 每个线程搬运4个浮点数，完成搬运至As需要所有线程搬运ldg_a_num轮
+    const int ldg_b_num = BK * BN / thread_num / 4; // 每个线程搬运4个浮点数，完成搬运至Bs需要所有线程搬运ldg_b_num轮
+
+    int a_tile_row = threadIdx.x / (BK / 4); // 每行4个字节作为一个内存块，当前线程负责第a_tile_row行的第a_tile_col个内存块的搬运
+    int a_tile_col = threadIdx.x % (BK / 4) * 4;
+    int a_tile_stride = BM / ldg_a_num; // 一共BM行，搬运ldg_a_num轮，每论搬运a_tile_stride行
+
+    int b_tile_row = threadIdx.x / (BN / 4); // 每行4个字节作为一个内存块，当前线程负责第b_tile_row行的第b_tile_col个内存块的搬运
+    int b_tile_col = threadIdx.x % (BN / 4) * 4;
+    int b_tile_stride = BK / ldg_b_num; // 一共BK行，搬运ldg_b_num轮，每论搬运b_tile_stride行
+
+    float accum[TM][TN] = {0.}; // 每个线程负责TM*TN个元素，则需要申请TM*TN个寄存器保存累加值，额外的一个寄存器用于缓存；
+
+    // 计算ldg_a_num的所有参数必须全部是const，否则不能用来申明数组大小
+    float ldg_a_reg[4 * ldg_a_num] = {0.}; // 每个线程搬运ldg_a_num轮，寄存器缓存ldg_a_num个float4元素，用于转置As矩阵
+
+    float a_frag[TM];  // 缓存As共享内存
+    float b_frag[TN];  // 缓存Bs共享内存
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+#pragma unroll
+    for (int k = 0; k < K; k += BK) {
+#pragma unroll
+        for (int i = 0; i < BM; i += a_tile_stride) {
+            int ldg_index = i / a_tile_stride * 4;  // 第ldg_index轮
+            FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
+                    FETCH_FLOAT4(A[OFFSET(a_tile_row + i, a_tile_col, K)]);
+            // As转置存，其中ldg_a_reg做中间缓存，目的是读取时可以按FLOAT4读取
+            As[OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
+            As[OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
+            As[OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
+            As[OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
+        }
+#pragma unroll
+        for (int i = 0; i < BK; i += b_tile_stride) {
+            FETCH_FLOAT4(Bs[OFFSET(b_tile_row + i, b_tile_col, BN)]) =
+                    FETCH_FLOAT4(B[OFFSET(b_tile_row + i, b_tile_col, N)]); // 不需要转置
+        }
+        __syncthreads();
+        A += BK;
+        B += BK * N;
+#pragma unroll
+        for (int i = 0; i < BK; i++) {
+#pragma unroll
+            for (int m = 0; m < TM; m += 4) {
+                FETCH_FLOAT4(a_frag[m]) = FETCH_FLOAT4(As[OFFSET(i, ty + m, BM)]); // 偏移到当前thread tile
+            }
+#pragma unroll
+            for (int n = 0; n < TN; n += 4) {
+                FETCH_FLOAT4(b_frag[n]) = FETCH_FLOAT4(Bs[OFFSET(i, tx + n, BN)]); // 偏移到当前thread tile
+            }
+#pragma unroll
+            for (int m = 0; m < TM; m++) {
+#pragma unroll
+                for (int n = 0; n < TN; n++) {
+                    accum[m][n] += a_frag[m] * b_frag[n];
+                }
+            }
+        }
+        __syncthreads();
+    }
+#pragma unroll
+    for (int m = 0; m < TM; m++) {
+#pragma unroll
+        for (int n = 0; n < TN; n += 4) {
+            float4 ctmp = FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]);
+            //float4 atmp = FETCH_FLOAT4(accum[m][n]);
+            ctmp.x = alpha * accum[m][n] + beta * ctmp.x;
+            ctmp.y = alpha * accum[m][n + 1] + beta * ctmp.y;
+            ctmp.z = alpha * accum[m][n + 2] + beta * ctmp.z;
+            ctmp.w = alpha * accum[m][n + 3] + beta * ctmp.w;
+            FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]) = ctmp;
+        }
+    }
+}
--- a/upstream_ref/sgemm_cuda/kernel_7.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_7.cuh
@@ -0,0 +1,180 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define OFFSET(row, col, ld) ((row)*(ld)+(col))
+#define FETCH_FLOAT4(pointer) (reinterpret_cast<float4*>(&(pointer))[0])
+
+template<const int BM,
+        const int BN,
+        const int BK,
+        const int TM,
+        const int TN>
+__global__ void mysgemm_v7(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    int bx = blockIdx.x;
+    int by = blockIdx.y;
+
+    const int block_row_thread = BN / TN;
+    const int block_col_thread = BM / TM;
+    const int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
+
+    // 当前线程对应thread tile的左上角元素在block中的位置
+    int tx = (threadIdx.x % block_row_thread) * TN;
+    int ty = (threadIdx.x / block_row_thread) * TM;
+
+    __shared__ float As[2][BK * BM]; // 增加一倍共享内存大小用于缓存
+    __shared__ float Bs[2][BK * BN];
+
+
+    const int ldg_a_num = BK * BM / thread_num / 4; // 每个线程搬运4个浮点数，完成搬运至As需要所有线程搬运ldg_a_num轮
+    const int ldg_b_num = BK * BN / thread_num / 4; // 每个线程搬运4个浮点数，完成搬运至Bs需要所有线程搬运ldg_b_num轮
+
+    int a_tile_row = threadIdx.x / (BK / 4); // 每行4个字节作为一个内存块，当前线程负责第a_tile_row行的第a_tile_col个内存块的搬运
+    int a_tile_col = threadIdx.x % (BK / 4) * 4;
+    int a_tile_stride = BM / ldg_a_num; // 一共BM行，搬运ldg_a_num轮，每论搬运a_tile_stride行
+
+    int b_tile_row = threadIdx.x / (BN / 4); // 每行4个字节作为一个内存块，当前线程负责第b_tile_row行的第b_tile_col个内存块的搬运
+    int b_tile_col = threadIdx.x % (BN / 4) * 4;
+    int b_tile_stride = BK / ldg_b_num; // 一共BK行，搬运ldg_b_num轮，每论搬运b_tile_stride行
+
+    float accum[TM][TN] = {0.}; // 每个线程负责TM*TN个元素，则需要申请TM*TN个寄存器保存累加值，额外的一个寄存器用于缓存；
+
+    // 计算ldg_a_num的所有参数必须全部是const，否则不能用来申明数组大小
+    float ldg_a_reg[4 * ldg_a_num] = {0.}; // 每个线程搬运ldg_a_num轮，寄存器缓存ldg_a_num个float4元素，用于转置As矩阵
+    float ldg_b_reg[4 * ldg_b_num] = {0.}; // 每个线程搬运ldg_a_num轮，寄存器缓存ldg_a_num个float4元素，用于转置As矩阵
+
+    float a_frag[2][TM];  // 缓存As共享内存,增加一倍寄存器大小用于缓存
+    float b_frag[2][TN];  // 缓存Bs共享内存,增加一倍寄存器大小用于缓存
+
+    // 移动到当前block
+    A = &A[by * BM * K];
+    B = &B[bx * BN];
+    C = &C[by * BM * N + bx * BN];
+
+    // first global to shared
+#pragma unroll
+    for (int i = 0; i < BM; i += a_tile_stride) {
+        int ldg_index = i / a_tile_stride * 4;  // 第ldg_index轮
+        FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
+                FETCH_FLOAT4(A[OFFSET(a_tile_row + i, a_tile_col, K)]);
+        // As转置存，其中ldg_a_reg做中间缓存，目的是读取时可以按FLOAT4读取
+        As[0][OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
+        As[0][OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
+        As[0][OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
+        As[0][OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
+    }
+#pragma unroll
+    for (int i = 0; i < BK; i += b_tile_stride) {
+        FETCH_FLOAT4(Bs[0][OFFSET(b_tile_row + i, b_tile_col, BN)]) =
+                FETCH_FLOAT4(B[OFFSET(b_tile_row + i, b_tile_col, N)]); // 不需要转置
+    }
+    __syncthreads();
+
+    // first shared to frag
+#pragma unroll
+    for (int m = 0; m < TM; m += 4) {
+        FETCH_FLOAT4(a_frag[0][m]) = FETCH_FLOAT4(As[0][OFFSET(0, ty + m, BM)]); // 偏移到当前thread tile
+    }
+#pragma unroll
+    for (int n = 0; n < TN; n += 4) {
+        FETCH_FLOAT4(b_frag[0][n]) = FETCH_FLOAT4(Bs[0][OFFSET(0, tx + n, BN)]); // 偏移到当前thread tile
+    }
+
+
+    int write_index = 1;
+    int load_index;
+    int k = 0;
+    do {
+        k += BK;
+        // load global to reg
+        if (k < K) {
+#pragma unroll
+            for (int i = 0; i < BM; i += a_tile_stride) {
+                int ldg_index = i / a_tile_stride * 4;  // 第ldg_index轮
+                FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
+                        FETCH_FLOAT4(A[OFFSET(a_tile_row + i, k + a_tile_col, K)]);
+            }
+#pragma unroll
+            for (int i = 0; i < BK; i += b_tile_stride) {
+                int ldg_index = i / b_tile_stride * 4;  // 第ldg_index轮
+                FETCH_FLOAT4(ldg_b_reg[ldg_index]) =
+                        FETCH_FLOAT4(B[OFFSET(k + b_tile_row + i, b_tile_col, N)]);
+            }
+        }
+
+        load_index = write_index ^ 1;
+#pragma unroll
+        for (int bk = 0; bk < BK - 1; bk++) {
+            for (int m = 0; m < TM; m += 4) {
+                FETCH_FLOAT4(a_frag[(bk + 1) % 2][m]) = FETCH_FLOAT4(
+                        As[load_index][OFFSET(bk + 1, ty + m, BM)]); // 偏移到当前thread tile
+            }
+#pragma unroll
+            for (int n = 0; n < TN; n += 4) {
+                FETCH_FLOAT4(b_frag[(bk + 1) % 2][n]) = FETCH_FLOAT4(
+                        Bs[load_index][OFFSET(bk + 1, tx + n, BN)]); // 偏移到当前thread tile
+            }
+#pragma unroll
+            for (int m = 0; m < TM; m++) {
+                for (int n = 0; n < TN; n++) {
+                    accum[m][n] += a_frag[bk % 2][m] * b_frag[bk % 2][n];
+                }
+            }
+        }
+        if (k < K) {
+#pragma unroll
+            for (int i = 0; i < BM; i += a_tile_stride) {
+                int ldg_index = i / a_tile_stride * 4;
+                As[write_index][OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
+                As[write_index][OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
+                As[write_index][OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
+                As[write_index][OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
+            }
+#pragma unroll
+            for (int i = 0; i < BK; i += b_tile_stride) {
+                int ldg_index = i / b_tile_stride * 4;
+                FETCH_FLOAT4(Bs[write_index][OFFSET(b_tile_row + i, b_tile_col, BN)]) =
+                        FETCH_FLOAT4(ldg_b_reg[ldg_index]);
+            }
+            __syncthreads();
+#pragma unroll
+            for (int m = 0; m < TM; m += 4) {
+                FETCH_FLOAT4(a_frag[0][m]) = FETCH_FLOAT4(
+                        As[write_index][OFFSET(0, ty + m, BM)]); // 偏移到当前thread tile
+            }
+#pragma unroll
+            for (int n = 0; n < TN; n += 4) {
+                FETCH_FLOAT4(b_frag[0][n]) = FETCH_FLOAT4(
+                        Bs[write_index][OFFSET(0, tx + n, BN)]); // 偏移到当前thread tile
+            }
+
+            write_index ^= 1;
+        }
+#pragma unroll
+        for (int m = 0; m < TM; m++) {
+#pragma unroll
+            for (int n = 0; n < TN; n++) {
+                accum[m][n] += a_frag[(BK - 1) % 2][m] * b_frag[(BK - 1) % 2][n];
+            }
+        }
+
+
+    } while (k < K);
+    
+    // C = alpha*AB+C
+#pragma unroll
+    for (int m = 0; m < TM; m++) {
+#pragma unroll
+        for (int n = 0; n < TN; n += 4) {
+            float4 ctmp = FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]);
+            ctmp.x = alpha * accum[m][n] + beta * ctmp.x;
+            ctmp.y = alpha * accum[m][n + 1] + beta * ctmp.y;
+            ctmp.z = alpha * accum[m][n + 2] + beta * ctmp.z;
+            ctmp.w = alpha * accum[m][n + 3] + beta * ctmp.w;
+            FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]) = ctmp;
+        }
+    }
+}
--- a/upstream_ref/sgemm_cuda/kernel_wangzyon.cuh
+++ b/upstream_ref/sgemm_cuda/kernel_wangzyon.cuh
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "kernel/kernel_1.cuh"
+#include "kernel/kernel_2.cuh"
+#include "kernel/kernel_3.cuh"
+#include "kernel/kernel_4.cuh"
+#include "kernel/kernel_5.cuh"
+#include "kernel/kernel_6.cuh"
+#include "kernel/kernel_7.cuh"
--- a/upstream_ref/sgemm_cuda/kernels.cuh
+++ b/upstream_ref/sgemm_cuda/kernels.cuh
@@ -0,0 +1,14 @@
+#pragma once
+
+#include "kernels/10_kernel_warptiling.cuh"
+#include "kernels/11_kernel_double_buffering.cuh"
+#include "kernels/12_kernel_double_buffering.cuh"
+#include "kernels/1_naive.cuh"
+#include "kernels/2_kernel_global_mem_coalesce.cuh"
+#include "kernels/3_kernel_shared_mem_blocking.cuh"
+#include "kernels/4_kernel_1D_blocktiling.cuh"
+#include "kernels/5_kernel_2D_blocktiling.cuh"
+#include "kernels/6_kernel_vectorize.cuh"
+#include "kernels/7_kernel_resolve_bank_conflicts.cuh"
+#include "kernels/8_kernel_bank_extra_col.cuh"
+#include "kernels/9_kernel_autotuned.cuh"
--- a/upstream_ref/sgemm_cuda/runner.cu
+++ b/upstream_ref/sgemm_cuda/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 (isnan(diff) || 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");
+  }
+}
--- a/upstream_ref/sgemm_cuda/runner.cuh
+++ b/upstream_ref/sgemm_cuda/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);
--- a/upstream_ref/sgemm_cuda/sgemm_siboehm.cu
+++ b/upstream_ref/sgemm_cuda/sgemm_siboehm.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;
+};
--- a/upstream_ref/sgemm_cuda/sgemm_wangzyon.cu
+++ b/upstream_ref/sgemm_cuda/sgemm_wangzyon.cu
@@ -0,0 +1,119 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <utils.cuh>
+
+#define cudaCheck(err) (cudaCheck(err, __FILE__, __LINE__))
+
+int main(int argc, char **argv) {
+    if (argc != 2) {
+        printf("Please select a kernel (range 0 - 11, here 0 is for NVIDIA cuBLAS).\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // cuda kernel num
+    int kernel_num = atoi(argv[1]);
+    if (kernel_num < 0 || kernel_num > 11) {
+        printf("Please enter a valid kernel number (0-11).\n");
+        exit(EXIT_FAILURE);
+    } else {
+        printf("Select kernel %d.\n", kernel_num);
+    };
+
+    // 申明句柄，创建句柄, cublasCreate会返回一个cublasStatus_t类型的值，用来判断句柄是否创建成功(值为0)
+    cublasHandle_t handle;
+    if (cublasCreate(&handle)) {
+        printf("Create cublas handle error.\n");
+        exit(EXIT_FAILURE);
+    };
+
+    // 采用cudaEvent进行gpu流计时，cudaEvent相当于在目标流中发布事件任务
+    float elapsed_time;
+    cudaEvent_t beg, end;
+    cudaEventCreate(&beg);
+    cudaEventCreate(&end);
+
+    // matrix size
+    int size_len = 24;
+    int SIZE[size_len];
+    for (int i = 0; i < size_len; i++)
+        SIZE[i] = 256 * (i + 1);
+
+    int m, n, k, max_size;
+    max_size = SIZE[size_len - 1];
+    printf("max_size=%d\n", max_size);
+
+    float alpha = 1.0, beta = 0.; //two arbitary input parameters，C=α*AB+β*C
+
+    float *A = NULL, *B = NULL, *C = NULL, *C_ref = NULL;     //host matrices
+    float *dA = NULL, *dB = NULL, *dC = NULL, *dC_ref = NULL; //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);
+    copy_matrix(C, C_ref, 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_ref, sizeof(float) * max_size * max_size, cudaMemcpyHostToDevice));
+
+    int repeat_times = 10;
+    for (int i = 0; i < size_len; i++) {
+        m = n = k = SIZE[i];
+
+        printf("m=n=k=%d\n", m);
+        // 验证计算正确性，同时在核函数计时前预先执行一次，避免冷启动误差
+        if (kernel_num != 0) {
+            test_kernel(0, m, n, k, alpha, dA, dB, beta, dC_ref, handle);      // cuBLAS
+            test_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC, handle); // user define
+            cudaDeviceSynchronize();
+            cudaMemcpy(C, dC, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
+            cudaMemcpy(C_ref, dC_ref, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
+            cudaDeviceSynchronize();
+
+            if (!verify_matrix(C_ref, C, m * n)) {
+                printf("Failed to pass the correctness verification against NVIDIA cuBLAS. Exited.\n");
+                exit(EXIT_FAILURE);
+            }
+        }
+        cudaDeviceSynchronize();
+
+        cudaEventRecord(beg);
+        for (int j = 0; j < repeat_times; j++) {
+            test_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.; //换算成秒
+
+        printf("Average elasped time: (%f) second, performance: (%f) GFLOPS. size: (%d).\n",
+               elapsed_time / repeat_times, 2. * 1e-9 * repeat_times * m * n * k / elapsed_time, m);
+        fflush(stdout);
+        copy_matrix(C_ref, C, m * n); //sync C with cuBLAS to prepare for the next run
+    }
+
+    // 释放CPU和GPU空间
+    free(A);
+    free(B);
+    free(C);
+    free(C_ref);
+    cudaFree(dA);
+    cudaFree(dB);
+    cudaFree(dC);
+    cudaFree(dC_ref);
+
+    return 0;
+};
--- a/upstream_ref/sgemm_cuda/simplest_kernel.cu
+++ b/upstream_ref/sgemm_cuda/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);
+}
--- a/upstream_ref/sgemm_cuda/utils.cu
+++ b/upstream_ref/sgemm_cuda/utils.cu
@@ -0,0 +1,199 @@
+#include <stdio.h>
+#include "utils.cuh"
+#include "kernel.cuh"
+
+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(line %d):\n%s\n", file, line, cudaGetErrorString(error));
+        exit(EXIT_FAILURE);
+    }
+    return;
+};
+
+void CudaDeviceInfo() {
+    int deviceId;
+
+    cudaGetDevice(&deviceId);
+
+    cudaDeviceProp props;
+    cudaGetDeviceProperties(&props, deviceId);
+
+    /*
+   * There should be no need to modify the output string below.
+   */
+
+    printf("Device ID: %d\n\
+       *Number of SMs: %d\n\
+       Compute Capability Major: %d\n\
+       Compute Capability Minor: %d\n\
+       memoryBusWidth: %d\n\
+       *maxThreadsPerBlock: %d\n\
+       maxThreadsPerMultiProcessor: %d\n\
+       *totalGlobalMem: %zuM\n\
+       sharedMemPerBlock: %zuKB\n\
+       *sharedMemPerMultiprocessor: %zuKB\n\
+       totalConstMem: %zuKB\n\
+       *multiProcessorCount: %d\n\
+       *Warp Size: %d\n",
+           deviceId,
+           props.multiProcessorCount,
+           props.major,
+           props.minor,
+           props.memoryBusWidth,
+           props.maxThreadsPerBlock,
+           props.maxThreadsPerMultiProcessor,
+           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: 使用gettimeofdays替代srand((unsigned)time(NULL));time精度过低，产生相同随机数
+    struct timeval time;
+    gettimeofday(&time, NULL);
+    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 copy_matrix(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) {
+    int i;
+    printf("[");
+    for (i = 0; i < M * N; i++) {
+        if ((i + 1) % N == 0)
+            printf("%5.2f ", A[i]);
+        else
+            printf("%5.2f, ", A[i]);
+        if ((i + 1) % N == 0) {
+            if (i + 1 < M * N)
+                printf(";\n");
+        }
+    }
+    printf("]\n");
+}
+
+bool verify_matrix(float *mat1, float *mat2, int N) {
+    double diff = 0.0;
+    int i;
+    for (i = 0; mat1 + i && mat2 + i && i < N; i++) {
+        diff = fabs((double) mat1[i] - (double) mat2[i]);
+        if (diff > 1e-2) {
+            printf("error. %5.2f,%5.2f,%d\n", mat1[i], mat2[i], i);
+            return false;
+        }
+    }
+    return true;
+}
+
+#define CEIL_DIV(M, N) ((M) + (N)-1) / (N)
+
+void test_cublas(cublasHandle_t handle, int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    //cublas列主序计算：https://www.cnblogs.com/cuancuancuanhao/p/7763256.html
+    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, N, A, K, &beta, C, N);
+}
+
+void test_mysgemm_v1(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(32, 32);
+    dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+    mysgemm_v1<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void test_mysgemm_v2(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(1024);
+    dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
+    mysgemm_v2<32><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void test_mysgemm_v3(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(512);
+    dim3 gridDim(CEIL_DIV(M, 64), CEIL_DIV(N, 64));
+    mysgemm_v3<64, 64, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void test_mysgemm_v4(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(256);
+    dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
+    mysgemm_v4<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void test_mysgemm_v5(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(256);
+    dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
+    mysgemm_v5<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+//void test_mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+//    dim3 blockDim(4);
+//    dim3 gridDim(CEIL_DIV(M, 8), CEIL_DIV(N, 8));
+//    mysgemm_v6<8, 8, 4, 4, 4><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+//}
+
+void test_mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(256);
+    dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
+    mysgemm_v6<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+void test_mysgemm_v7(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
+    dim3 blockDim(256);
+    dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
+    mysgemm_v7<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
+}
+
+
+
+void test_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:
+            test_cublas(handle, M, N, K, alpha, A, B, beta, C);
+            break;
+        case 1:
+            test_mysgemm_v1(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 2:
+            test_mysgemm_v2(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 3:
+            test_mysgemm_v3(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 4:
+            test_mysgemm_v4(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 5:
+            test_mysgemm_v5(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 6:
+            test_mysgemm_v6(M, N, K, alpha, A, B, beta, C);
+            break;
+        case 7:
+            test_mysgemm_v7(M, N, K, alpha, A, B, beta, C);
+            break;
+        default:
+            break;
+    }
+}
--- a/upstream_ref/sgemm_cuda/utils.cuh
+++ b/upstream_ref/sgemm_cuda/utils.cuh
@@ -0,0 +1,42 @@
+#pragma once
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <unistd.h>
+#include <sys/time.h>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+
+/*
+=====================================
+CUDA操作
+=====================================
+*/
+void cudaCheck(cudaError_t error, const char *file, int line); //CUDA错误检查
+void CudaDeviceInfo();                                         // 打印CUDA信息
+
+/*
+=====================================
+矩阵操作
+=====================================
+*/
+void randomize_matrix(float *mat, int N);            // 随机初始化矩阵
+void copy_matrix(float *src, float *dest, int N);    // 复制矩阵
+void print_matrix(const float *A, int M, int N);     // 打印矩阵
+bool verify_matrix(float *mat1, float *mat2, int N); // 验证矩阵
+
+/*
+=====================================
+计时操作
+=====================================
+*/
+float get_current_sec();                        // 获取当前时刻
+float cpu_elapsed_time(float &beg, float &end); // 计算时间差
+
+/*
+=====================================
+kernel操作
+=====================================
+*/
+//调用指定核函数计算矩阵乘法
+void test_kernel(int kernel_num, int m, int n, int k, float alpha, float *A, float *B, float beta, float *C, cublasHandle_t handle);