dylanyunlong/project_6
1
0
Fork 0
Code Issues Pull Requests Actions Projects Releases Wiki Activity

data: cat SGEMM files from 3 repos into cat_files/

Browse Source 
siboehm/SGEMM_CUDA (19 files):
  siboehm_sgemm.cu, siboehm_runner.cu, siboehm_runner.cuh, siboehm_kernels.cuh
  siboehm_cuBLAS_sgemm.cu, siboehm_simplest_kernel.cu, siboehm_CMakeLists.txt
  siboehm_{1_naive..12_kernel_double_buffering}.cuh

wangzyon/NVIDIA_SGEMM_PRACTICE (12 files):
  wangzyon_sgemm.cu, wangzyon_utils.cu, wangzyon_utils.cuh, wangzyon_kernel.cuh
  wangzyon_CMakeLists.txt, wangzyon_kernel_{1..7}.cuh

edtallison/sgemm-cuda (19 files):
  edtallison_sgemm.cu, edtallison_runner.cu, edtallison_runner.cuh
  edtallison_kernels.cuh, edtallison_cuBLAS_sgemm.cu, edtallison_simplest_kernel.cu
  edtallison_CMakeLists.txt, edtallison_{01_naive..12_kernel_double_buffering}.cuh

cat_files/ total: 25 → 75 files
This commit is contained in:
dylan
2026-08-15 06:59:18 +00:00
parent 284804ac53
commit 7cfa87b5ac
50 changed files with 5612 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
cat_files/edtallison_01_naive.cuh Normal file
View File

@@ -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];
}
}

25
cat_files/edtallison_02_kernel_global_mem_coalesce.cuh Normal file
View File

@@ -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];
}
}

54
cat_files/edtallison_03_kernel_shared_mem_blocking.cuh Normal file
View File

@@ -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];
}

80
cat_files/edtallison_04_kernel_1D_blocktiling.cuh Normal file
View File

@@ -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];
}
}

102
cat_files/edtallison_05_kernel_2D_blocktiling.cuh Normal file
View File

@@ -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];
}
}
}

98
cat_files/edtallison_06_kernel_vectorize.cuh Normal file
View File

@@ -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;
}
}
}

103
cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh Normal file
View File

@@ -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;
}
}
}

103
cat_files/edtallison_08_kernel_bank_extra_col.cuh Normal file
View File

@@ -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;
}
}
}

127
cat_files/edtallison_09_kernel_autotuned.cuh Normal file
View File

@@ -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;
}
}
}
}
}

187
cat_files/edtallison_10_kernel_warptiling.cuh Normal file
View File

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

220
cat_files/edtallison_11_kernel_double_buffering.cuh Normal file
View File

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

229
cat_files/edtallison_12_kernel_double_buffering.cuh Normal file
View File

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

36
cat_files/edtallison_CMakeLists.txt Normal file
View File

@@ -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})

108
cat_files/edtallison_cuBLAS_sgemm.cu Normal file
View File

@@ -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;
}

14
cat_files/edtallison_kernels.cuh Normal file
View File

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

549
cat_files/edtallison_runner.cu Normal file
View File

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

26
cat_files/edtallison_runner.cuh Normal file
View File

@@ -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);

168
cat_files/edtallison_sgemm.cu Normal file
View File

@@ -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;
};

46
cat_files/edtallison_simplest_kernel.cu Normal file
View File

@@ -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);
}

187
cat_files/siboehm_10_kernel_warptiling.cuh Normal file
View File

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

220
cat_files/siboehm_11_kernel_double_buffering.cuh Normal file
View File

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

229
cat_files/siboehm_12_kernel_double_buffering.cuh Normal file
View File

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

29
cat_files/siboehm_1_naive.cuh Normal file
View File

@@ -0,0 +1,29 @@
#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, float alpha, const float *A,
const float *B, float beta, float *C) {
const uint x = blockIdx.x * blockDim.x + threadIdx.x;
const uint y = blockIdx.y * blockDim.y + threadIdx.y;
// if statement is necessary to make things work under tile quantization
if (x < M && y < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) {
tmp += A[x * K + i] * B[i * N + y];
}
// C = α*(A@B)+β*C
C[x * N + y] = alpha * tmp + beta * C[x * N + y];
}
}

24
cat_files/siboehm_2_kernel_global_mem_coalesce.cuh Normal file
View File

@@ -0,0 +1,24 @@
#pragma once
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
template <const uint BLOCKSIZE>
__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);
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];
}
}

57
cat_files/siboehm_3_kernel_shared_mem_blocking.cuh Normal file
View File

@@ -0,0 +1,57 @@
#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) {
// the output block that we want to compute in this threadblock
const uint cRow = blockIdx.x;
const uint cCol = blockIdx.y;
// allocate buffer for current block in fast shared mem
// shared mem is shared between all threads in a block
__shared__ float As[BLOCKSIZE * BLOCKSIZE];
__shared__ float Bs[BLOCKSIZE * BLOCKSIZE];
// the inner row & col that we're accessing in this thread
const uint threadCol = threadIdx.x % BLOCKSIZE;
const uint threadRow = threadIdx.x / BLOCKSIZE;
// advance pointers to the starting positions
A += cRow * BLOCKSIZE * K; // row=cRow, col=0
B += cCol * BLOCKSIZE; // row=0, col=cCol
C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE; // row=cRow, col=cCol
float tmp = 0.0;
for (int bkIdx = 0; bkIdx < K; bkIdx += BLOCKSIZE) {
// Have each thread load one of the elements in A & B
// Make the threadCol (=threadIdx.x) the consecutive index
// to allow global memory access coalescing
As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
// block threads in this block until cache is fully populated
__syncthreads();
A += BLOCKSIZE;
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];
}
// need to sync again at the end, to avoid faster threads
// fetching the next block into the cache before slower threads are done
__syncthreads();
}
C[threadRow * N + threadCol] =
alpha * tmp + beta * C[threadRow * N + threadCol];
}

80
cat_files/siboehm_4_kernel_1D_blocktiling.cuh Normal file
View File

@@ -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];
}
}

102
cat_files/siboehm_5_kernel_2D_blocktiling.cuh Normal file
View File

@@ -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];
}
}
}

98
cat_files/siboehm_6_kernel_vectorize.cuh Normal file
View File

@@ -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;
}
}
}

103
cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh Normal file
View File

@@ -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;
}
}
}

103
cat_files/siboehm_8_kernel_bank_extra_col.cuh Normal file
View File

@@ -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;
}
}
}

127
cat_files/siboehm_9_kernel_autotuned.cuh Normal file
View File

@@ -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;
}
}
}
}
}

36
cat_files/siboehm_CMakeLists.txt Normal file
View File

@@ -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})

108
cat_files/siboehm_cuBLAS_sgemm.cu Normal file
View File

@@ -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;
}

14
cat_files/siboehm_kernels.cuh Normal file
View File

@@ -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"

549
cat_files/siboehm_runner.cu Normal file
View File

@@ -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");
}
}

26
cat_files/siboehm_runner.cuh Normal file
View File

@@ -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);

168
cat_files/siboehm_sgemm.cu Normal file
View File

@@ -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;
};

46
cat_files/siboehm_simplest_kernel.cu Normal file
View File

@@ -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);
}

36
cat_files/wangzyon_CMakeLists.txt Normal file
View File

@@ -0,0 +1,36 @@
cmake_minimum_required(VERSION 3.0)
project(NVIDIA_SGEMM_PRACTICE)
# gcc/g++编译参数说明：
# -O1~3编译器优化选项的4个级别，-O1默认，级别越大优化效果越好，但编译时间越长;
# -std=c++11，采用C++11标准编译
set(CMAKE_CXX_FLAGS "-O3 -std=c++11")
# nvcc编译参数说明：
# -g:主机代码添加调试信息;
# -G:设备代码产生调试信息，将会禁用大多数编译器优化，造成设备代码运行缓慢;
# -Xptxas -dlcm=ca启用L1缓存，-Xptxas -dlcm=cg关闭L1缓存
# set(CUDA_NVCC_FLAGS -g;-G;-Xptxas;-dlcm=ca)
# set(CUDA_NVCC_FLAGS -Xptxas;-dlcm=cg)
set(CUDA_NVCC_FLAGS -arch=compute_70;-code=compute_70)
# 若FIND CUDA ERROR，在~/.bashrc中添加配置环境变量和动态库路径
# CUDA_HOME=/usr/local/cuda
# export PATH=$CUDA_HOME/bin:$PATH
# export LD_LIBRARY_PATH=$CUDA_HOME/lib64:$LD_LIBRARY_PATH
find_package(CUDA REQUIRED)
# 配置头文件搜索路径
include_directories(${CUDA_INCLUDE_DIRS})
include_directories(${PROJECT_SOURCE_DIR}/src)
# 配置待编译的源文件路径
aux_source_directory(${PROJECT_SOURCE_DIR}/src SRC)
# 可执行文件输出路径
set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR})
# 生成可执行文件
CUDA_ADD_EXECUTABLE(sgemm sgemm.cu ${SRC})
# link cudart cublas
target_link_libraries(sgemm ${CUDA_LIBRARIES} ${CUDA_cublas_LIBRARY})

9
cat_files/wangzyon_kernel.cuh Normal file
View File

@@ -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"

19
cat_files/wangzyon_kernel_1.cuh Normal file
View File

@@ -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];
}

45
cat_files/wangzyon_kernel_2.cuh Normal file
View File

@@ -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];
}

71
cat_files/wangzyon_kernel_3.cuh Normal file
View File

@@ -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];
}
}

76
cat_files/wangzyon_kernel_4.cuh Normal file
View File

@@ -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];
}
}

88
cat_files/wangzyon_kernel_5.cuh Normal file
View File

@@ -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];
}
}

110
cat_files/wangzyon_kernel_6.cuh Normal file
View File

@@ -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;
}
}
}

180
cat_files/wangzyon_kernel_7.cuh Normal file
View File

@@ -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;
}
}
}

119
cat_files/wangzyon_sgemm.cu Normal file
View File

@@ -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;
};

199
cat_files/wangzyon_utils.cu Normal file
View File

@@ -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;
}
}

42
cat_files/wangzyon_utils.cuh Normal file
View File

@@ -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);