From 7cfa87b5ac0b09be78cc91850efb1033b02f888b Mon Sep 17 00:00:00 2001 From: dylan Date: Sat, 15 Aug 2026 06:59:18 +0000 Subject: [PATCH] data: cat SGEMM files from 3 repos into cat_files/ MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- cat_files/edtallison_01_naive.cuh | 37 ++ ...tallison_02_kernel_global_mem_coalesce.cuh | 25 + ...tallison_03_kernel_shared_mem_blocking.cuh | 54 ++ .../edtallison_04_kernel_1D_blocktiling.cuh | 80 +++ .../edtallison_05_kernel_2D_blocktiling.cuh | 102 ++++ cat_files/edtallison_06_kernel_vectorize.cuh | 98 ++++ ...lison_07_kernel_resolve_bank_conflicts.cuh | 103 ++++ .../edtallison_08_kernel_bank_extra_col.cuh | 103 ++++ cat_files/edtallison_09_kernel_autotuned.cuh | 127 ++++ cat_files/edtallison_10_kernel_warptiling.cuh | 187 ++++++ .../edtallison_11_kernel_double_buffering.cuh | 220 +++++++ .../edtallison_12_kernel_double_buffering.cuh | 229 ++++++++ cat_files/edtallison_CMakeLists.txt | 36 ++ cat_files/edtallison_cuBLAS_sgemm.cu | 108 ++++ cat_files/edtallison_kernels.cuh | 14 + cat_files/edtallison_runner.cu | 549 ++++++++++++++++++ cat_files/edtallison_runner.cuh | 26 + cat_files/edtallison_sgemm.cu | 168 ++++++ cat_files/edtallison_simplest_kernel.cu | 46 ++ cat_files/siboehm_10_kernel_warptiling.cuh | 187 ++++++ .../siboehm_11_kernel_double_buffering.cuh | 220 +++++++ .../siboehm_12_kernel_double_buffering.cuh | 229 ++++++++ cat_files/siboehm_1_naive.cuh | 29 + .../siboehm_2_kernel_global_mem_coalesce.cuh | 24 + .../siboehm_3_kernel_shared_mem_blocking.cuh | 57 ++ cat_files/siboehm_4_kernel_1D_blocktiling.cuh | 80 +++ cat_files/siboehm_5_kernel_2D_blocktiling.cuh | 102 ++++ cat_files/siboehm_6_kernel_vectorize.cuh | 98 ++++ ...iboehm_7_kernel_resolve_bank_conflicts.cuh | 103 ++++ cat_files/siboehm_8_kernel_bank_extra_col.cuh | 103 ++++ cat_files/siboehm_9_kernel_autotuned.cuh | 127 ++++ cat_files/siboehm_CMakeLists.txt | 36 ++ cat_files/siboehm_cuBLAS_sgemm.cu | 108 ++++ cat_files/siboehm_kernels.cuh | 14 + cat_files/siboehm_runner.cu | 549 ++++++++++++++++++ cat_files/siboehm_runner.cuh | 26 + cat_files/siboehm_sgemm.cu | 168 ++++++ cat_files/siboehm_simplest_kernel.cu | 46 ++ cat_files/wangzyon_CMakeLists.txt | 36 ++ cat_files/wangzyon_kernel.cuh | 9 + cat_files/wangzyon_kernel_1.cuh | 19 + cat_files/wangzyon_kernel_2.cuh | 45 ++ cat_files/wangzyon_kernel_3.cuh | 71 +++ cat_files/wangzyon_kernel_4.cuh | 76 +++ cat_files/wangzyon_kernel_5.cuh | 88 +++ cat_files/wangzyon_kernel_6.cuh | 110 ++++ cat_files/wangzyon_kernel_7.cuh | 180 ++++++ cat_files/wangzyon_sgemm.cu | 119 ++++ cat_files/wangzyon_utils.cu | 199 +++++++ cat_files/wangzyon_utils.cuh | 42 ++ 50 files changed, 5612 insertions(+) create mode 100644 cat_files/edtallison_01_naive.cuh create mode 100644 cat_files/edtallison_02_kernel_global_mem_coalesce.cuh create mode 100644 cat_files/edtallison_03_kernel_shared_mem_blocking.cuh create mode 100644 cat_files/edtallison_04_kernel_1D_blocktiling.cuh create mode 100644 cat_files/edtallison_05_kernel_2D_blocktiling.cuh create mode 100644 cat_files/edtallison_06_kernel_vectorize.cuh create mode 100644 cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh create mode 100644 cat_files/edtallison_08_kernel_bank_extra_col.cuh create mode 100644 cat_files/edtallison_09_kernel_autotuned.cuh create mode 100644 cat_files/edtallison_10_kernel_warptiling.cuh create mode 100644 cat_files/edtallison_11_kernel_double_buffering.cuh create mode 100644 cat_files/edtallison_12_kernel_double_buffering.cuh create mode 100644 cat_files/edtallison_CMakeLists.txt create mode 100644 cat_files/edtallison_cuBLAS_sgemm.cu create mode 100644 cat_files/edtallison_kernels.cuh create mode 100644 cat_files/edtallison_runner.cu create mode 100644 cat_files/edtallison_runner.cuh create mode 100644 cat_files/edtallison_sgemm.cu create mode 100644 cat_files/edtallison_simplest_kernel.cu create mode 100644 cat_files/siboehm_10_kernel_warptiling.cuh create mode 100644 cat_files/siboehm_11_kernel_double_buffering.cuh create mode 100644 cat_files/siboehm_12_kernel_double_buffering.cuh create mode 100644 cat_files/siboehm_1_naive.cuh create mode 100644 cat_files/siboehm_2_kernel_global_mem_coalesce.cuh create mode 100644 cat_files/siboehm_3_kernel_shared_mem_blocking.cuh create mode 100644 cat_files/siboehm_4_kernel_1D_blocktiling.cuh create mode 100644 cat_files/siboehm_5_kernel_2D_blocktiling.cuh create mode 100644 cat_files/siboehm_6_kernel_vectorize.cuh create mode 100644 cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh create mode 100644 cat_files/siboehm_8_kernel_bank_extra_col.cuh create mode 100644 cat_files/siboehm_9_kernel_autotuned.cuh create mode 100644 cat_files/siboehm_CMakeLists.txt create mode 100644 cat_files/siboehm_cuBLAS_sgemm.cu create mode 100644 cat_files/siboehm_kernels.cuh create mode 100644 cat_files/siboehm_runner.cu create mode 100644 cat_files/siboehm_runner.cuh create mode 100644 cat_files/siboehm_sgemm.cu create mode 100644 cat_files/siboehm_simplest_kernel.cu create mode 100644 cat_files/wangzyon_CMakeLists.txt create mode 100644 cat_files/wangzyon_kernel.cuh create mode 100644 cat_files/wangzyon_kernel_1.cuh create mode 100644 cat_files/wangzyon_kernel_2.cuh create mode 100644 cat_files/wangzyon_kernel_3.cuh create mode 100644 cat_files/wangzyon_kernel_4.cuh create mode 100644 cat_files/wangzyon_kernel_5.cuh create mode 100644 cat_files/wangzyon_kernel_6.cuh create mode 100644 cat_files/wangzyon_kernel_7.cuh create mode 100644 cat_files/wangzyon_sgemm.cu create mode 100644 cat_files/wangzyon_utils.cu create mode 100644 cat_files/wangzyon_utils.cuh diff --git a/cat_files/edtallison_01_naive.cuh b/cat_files/edtallison_01_naive.cuh new file mode 100644 index 00000000..04ed8e8c --- /dev/null +++ b/cat_files/edtallison_01_naive.cuh @@ -0,0 +1,37 @@ +# pragma once + +#include +#include +#include +#include + +/* + +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]; + } +} diff --git a/cat_files/edtallison_02_kernel_global_mem_coalesce.cuh b/cat_files/edtallison_02_kernel_global_mem_coalesce.cuh new file mode 100644 index 00000000..2f6f4504 --- /dev/null +++ b/cat_files/edtallison_02_kernel_global_mem_coalesce.cuh @@ -0,0 +1,25 @@ +#pragma once + +#include +#include +#include +#include +#include + +template +// __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]; + } +} diff --git a/cat_files/edtallison_03_kernel_shared_mem_blocking.cuh b/cat_files/edtallison_03_kernel_shared_mem_blocking.cuh new file mode 100644 index 00000000..975cc47f --- /dev/null +++ b/cat_files/edtallison_03_kernel_shared_mem_blocking.cuh @@ -0,0 +1,54 @@ +#pragma once + +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; +} \ No newline at end of file diff --git a/cat_files/edtallison_04_kernel_1D_blocktiling.cuh b/cat_files/edtallison_04_kernel_1D_blocktiling.cuh new file mode 100644 index 00000000..12c9c5d7 --- /dev/null +++ b/cat_files/edtallison_04_kernel_1D_blocktiling.cuh @@ -0,0 +1,80 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/edtallison_05_kernel_2D_blocktiling.cuh b/cat_files/edtallison_05_kernel_2D_blocktiling.cuh new file mode 100644 index 00000000..0b361f7c --- /dev/null +++ b/cat_files/edtallison_05_kernel_2D_blocktiling.cuh @@ -0,0 +1,102 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_06_kernel_vectorize.cuh b/cat_files/edtallison_06_kernel_vectorize.cuh new file mode 100644 index 00000000..665d4e50 --- /dev/null +++ b/cat_files/edtallison_06_kernel_vectorize.cuh @@ -0,0 +1,98 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&Bs[innerRowB * BN + innerColB * 4])[0] = + reinterpret_cast(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh b/cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh new file mode 100644 index 00000000..d571008e --- /dev/null +++ b/cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh @@ -0,0 +1,103 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_08_kernel_bank_extra_col.cuh b/cat_files/edtallison_08_kernel_bank_extra_col.cuh new file mode 100644 index 00000000..c362cba8 --- /dev/null +++ b/cat_files/edtallison_08_kernel_bank_extra_col.cuh @@ -0,0 +1,103 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_09_kernel_autotuned.cuh b/cat_files/edtallison_09_kernel_autotuned.cuh new file mode 100644 index 00000000..6841952b --- /dev/null +++ b/cat_files/edtallison_09_kernel_autotuned.cuh @@ -0,0 +1,127 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) +const int K9_NUM_THREADS = 256; + +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &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( + &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(&C_interim[(threadRow * TM + resIdxM) * N + + threadCol * TN + resIdxN])[0] = + tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_10_kernel_warptiling.cuh b/cat_files/edtallison_10_kernel_warptiling.cuh new file mode 100644 index 00000000..2cc66f36 --- /dev/null +++ b/cat_files/edtallison_10_kernel_warptiling.cuh @@ -0,0 +1,187 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) +const int WARPSIZE = 32; // warpSize is not constexpr + +namespace wt { +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &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 +__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 +__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( + N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB); + __syncthreads(); + wt::processFromSmem(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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_11_kernel_double_buffering.cuh b/cat_files/edtallison_11_kernel_double_buffering.cuh new file mode 100644 index 00000000..f54b61cf --- /dev/null +++ b/cat_files/edtallison_11_kernel_double_buffering.cuh @@ -0,0 +1,220 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +namespace db { + +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &B[(innerRowB + offset) * N + innerColB * 4])[0]; + } +} + +template +__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 +__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( + 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(regM, regN, threadResults, As, Bs, warpRow, + warpCol, threadRowInWarp, threadColInWarp); + __syncthreads(); + + // process current+1 (B1) + if (bkIdx + BK < K) { + db::processFromSmem(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( + 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( + N, K, A + BK, B + BK * N, As + (BM * BK), Bs + (BK * BN), innerRowA, + innerColA, innerRowB, innerColB); + } + __syncthreads(); + + // process current (B0) + db::processFromSmem(regM, regN, threadResults, As, Bs, warpRow, + warpCol, threadRowInWarp, threadColInWarp); + __syncthreads(); + + // process current+1 (B1) + if (bkIdx + BK < K) { + db::processFromSmem(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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_12_kernel_double_buffering.cuh b/cat_files/edtallison_12_kernel_double_buffering.cuh new file mode 100644 index 00000000..c0b0a40f --- /dev/null +++ b/cat_files/edtallison_12_kernel_double_buffering.cuh @@ -0,0 +1,229 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +namespace { +template +__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)), + barrier); + cuda::memcpy_async(&As[(innerColA * 4 + 1) * BM + innerRowA + offset], + &A[(innerRowA + offset) * K + innerColA * 4 + 1], + cuda::aligned_size_t(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)), + barrier); + cuda::memcpy_async(&As[(innerColA * 4 + 3) * BM + innerRowA + offset], + &A[(innerRowA + offset) * K + innerColA * 4 + 3], + cuda::aligned_size_t(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)), + barrier); + } +} + +template +__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 +__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 frontBarrier; + __shared__ cuda::barrier 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( + 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( + 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( + 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( + 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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/edtallison_CMakeLists.txt b/cat_files/edtallison_CMakeLists.txt new file mode 100644 index 00000000..fde93413 --- /dev/null +++ b/cat_files/edtallison_CMakeLists.txt @@ -0,0 +1,36 @@ +cmake_minimum_required(VERSION 3.19) +project(NVIDIA_SGEMM_PRACTICE LANGUAGES CXX CUDA) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +find_package(CUDA REQUIRED) + +# ensure cuda is available +include(CheckLanguage) +check_language(CUDA) + +set(CMAKE_CXX_STANDARD 20) +set(CUDA_COMPUTE_CAPABILITY 75) + +# in debug mode, add debug symbols to device code +# this disables most optimizations and kills performance +add_compile_options("$<$,$>:-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}) \ No newline at end of file diff --git a/cat_files/edtallison_cuBLAS_sgemm.cu b/cat_files/edtallison_cuBLAS_sgemm.cu new file mode 100644 index 00000000..c6062184 --- /dev/null +++ b/cat_files/edtallison_cuBLAS_sgemm.cu @@ -0,0 +1,108 @@ +#include +#include +#include + +/* + * 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; +} \ No newline at end of file diff --git a/cat_files/edtallison_kernels.cuh b/cat_files/edtallison_kernels.cuh new file mode 100644 index 00000000..b027b017 --- /dev/null +++ b/cat_files/edtallison_kernels.cuh @@ -0,0 +1,14 @@ +#pragma once + +#include "kernels/01_naive.cuh" +#include "kernels/02_kernel_global_mem_coalesce.cuh" +#include "kernels/03_kernel_shared_mem_blocking.cuh" +#include "kernels/04_kernel_1D_blocktiling.cuh" +#include "kernels/05_kernel_2D_blocktiling.cuh" +#include "kernels/06_kernel_vectorize.cuh" +#include "kernels/07_kernel_resolve_bank_conflicts.cuh" +#include "kernels/08_kernel_bank_extra_col.cuh" +#include "kernels/09_kernel_autotuned.cuh" +#include "kernels/10_kernel_warptiling.cuh" +#include "kernels/11_kernel_double_buffering.cuh" +#include "kernels/12_kernel_double_buffering.cuh" \ No newline at end of file diff --git a/cat_files/edtallison_runner.cu b/cat_files/edtallison_runner.cu new file mode 100644 index 00000000..4c0142dc --- /dev/null +++ b/cat_files/edtallison_runner.cu @@ -0,0 +1,549 @@ +#include "kernels.cuh" +#include "runner.cuh" +#include +#include +#include +#include + +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<<>>(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> + <<>>(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> + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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"); + } +} \ No newline at end of file diff --git a/cat_files/edtallison_runner.cuh b/cat_files/edtallison_runner.cuh new file mode 100644 index 00000000..3ac967e6 --- /dev/null +++ b/cat_files/edtallison_runner.cuh @@ -0,0 +1,26 @@ +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include + +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); \ No newline at end of file diff --git a/cat_files/edtallison_sgemm.cu b/cat_files/edtallison_sgemm.cu new file mode 100644 index 00000000..bb6dd38c --- /dev/null +++ b/cat_files/edtallison_sgemm.cu @@ -0,0 +1,168 @@ +#include +#include +#include +#include +#include +#include +#include + +#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 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; +}; \ No newline at end of file diff --git a/cat_files/edtallison_simplest_kernel.cu b/cat_files/edtallison_simplest_kernel.cu new file mode 100644 index 00000000..2fc15ccc --- /dev/null +++ b/cat_files/edtallison_simplest_kernel.cu @@ -0,0 +1,46 @@ +#include +#include +#include + +__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<<>>(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); +} diff --git a/cat_files/siboehm_10_kernel_warptiling.cuh b/cat_files/siboehm_10_kernel_warptiling.cuh new file mode 100644 index 00000000..2cc66f36 --- /dev/null +++ b/cat_files/siboehm_10_kernel_warptiling.cuh @@ -0,0 +1,187 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) +const int WARPSIZE = 32; // warpSize is not constexpr + +namespace wt { +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &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 +__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 +__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( + N, K, A, B, As, Bs, innerRowA, innerColA, innerRowB, innerColB); + __syncthreads(); + wt::processFromSmem(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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_11_kernel_double_buffering.cuh b/cat_files/siboehm_11_kernel_double_buffering.cuh new file mode 100644 index 00000000..f54b61cf --- /dev/null +++ b/cat_files/siboehm_11_kernel_double_buffering.cuh @@ -0,0 +1,220 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +namespace db { + +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &B[(innerRowB + offset) * N + innerColB * 4])[0]; + } +} + +template +__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 +__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( + 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(regM, regN, threadResults, As, Bs, warpRow, + warpCol, threadRowInWarp, threadColInWarp); + __syncthreads(); + + // process current+1 (B1) + if (bkIdx + BK < K) { + db::processFromSmem(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( + 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( + N, K, A + BK, B + BK * N, As + (BM * BK), Bs + (BK * BN), innerRowA, + innerColA, innerRowB, innerColB); + } + __syncthreads(); + + // process current (B0) + db::processFromSmem(regM, regN, threadResults, As, Bs, warpRow, + warpCol, threadRowInWarp, threadColInWarp); + __syncthreads(); + + // process current+1 (B1) + if (bkIdx + BK < K) { + db::processFromSmem(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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_12_kernel_double_buffering.cuh b/cat_files/siboehm_12_kernel_double_buffering.cuh new file mode 100644 index 00000000..c0b0a40f --- /dev/null +++ b/cat_files/siboehm_12_kernel_double_buffering.cuh @@ -0,0 +1,229 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +namespace { +template +__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)), + barrier); + cuda::memcpy_async(&As[(innerColA * 4 + 1) * BM + innerRowA + offset], + &A[(innerRowA + offset) * K + innerColA * 4 + 1], + cuda::aligned_size_t(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)), + barrier); + cuda::memcpy_async(&As[(innerColA * 4 + 3) * BM + innerRowA + offset], + &A[(innerRowA + offset) * K + innerColA * 4 + 3], + cuda::aligned_size_t(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)), + barrier); + } +} + +template +__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 +__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 frontBarrier; + __shared__ cuda::barrier 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( + 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( + 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( + 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( + 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( + &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( + &C_interim[(threadRowInWarp * TM + resIdxM) * N + + threadColInWarp * TN + resIdxN])[0] = tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_1_naive.cuh b/cat_files/siboehm_1_naive.cuh new file mode 100644 index 00000000..47038a37 --- /dev/null +++ b/cat_files/siboehm_1_naive.cuh @@ -0,0 +1,29 @@ +#pragma once + +#include +#include +#include +#include + +/* + +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]; + } +} \ No newline at end of file diff --git a/cat_files/siboehm_2_kernel_global_mem_coalesce.cuh b/cat_files/siboehm_2_kernel_global_mem_coalesce.cuh new file mode 100644 index 00000000..ef0a8b94 --- /dev/null +++ b/cat_files/siboehm_2_kernel_global_mem_coalesce.cuh @@ -0,0 +1,24 @@ +#pragma once + +#include +#include +#include +#include +#include + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/siboehm_3_kernel_shared_mem_blocking.cuh b/cat_files/siboehm_3_kernel_shared_mem_blocking.cuh new file mode 100644 index 00000000..806b6982 --- /dev/null +++ b/cat_files/siboehm_3_kernel_shared_mem_blocking.cuh @@ -0,0 +1,57 @@ +#pragma once + +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; +} \ No newline at end of file diff --git a/cat_files/siboehm_4_kernel_1D_blocktiling.cuh b/cat_files/siboehm_4_kernel_1D_blocktiling.cuh new file mode 100644 index 00000000..12c9c5d7 --- /dev/null +++ b/cat_files/siboehm_4_kernel_1D_blocktiling.cuh @@ -0,0 +1,80 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/siboehm_5_kernel_2D_blocktiling.cuh b/cat_files/siboehm_5_kernel_2D_blocktiling.cuh new file mode 100644 index 00000000..0b361f7c --- /dev/null +++ b/cat_files/siboehm_5_kernel_2D_blocktiling.cuh @@ -0,0 +1,102 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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]; + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_6_kernel_vectorize.cuh b/cat_files/siboehm_6_kernel_vectorize.cuh new file mode 100644 index 00000000..665d4e50 --- /dev/null +++ b/cat_files/siboehm_6_kernel_vectorize.cuh @@ -0,0 +1,98 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&Bs[innerRowB * BN + innerColB * 4])[0] = + reinterpret_cast(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh b/cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh new file mode 100644 index 00000000..d571008e --- /dev/null +++ b/cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh @@ -0,0 +1,103 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_8_kernel_bank_extra_col.cuh b/cat_files/siboehm_8_kernel_bank_extra_col.cuh new file mode 100644 index 00000000..c362cba8 --- /dev/null +++ b/cat_files/siboehm_8_kernel_bank_extra_col.cuh @@ -0,0 +1,103 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) + +template +__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(&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(&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( + &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( + &C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] = + tmp; + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_9_kernel_autotuned.cuh b/cat_files/siboehm_9_kernel_autotuned.cuh new file mode 100644 index 00000000..6841952b --- /dev/null +++ b/cat_files/siboehm_9_kernel_autotuned.cuh @@ -0,0 +1,127 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#define CEIL_DIV(M, N) (((M) + (N)-1) / (N)) +const int K9_NUM_THREADS = 256; + +template +__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( + &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( + &Bs[(innerRowB + offset) * BN + innerColB * 4])[0] = + reinterpret_cast( + &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( + &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(&C_interim[(threadRow * TM + resIdxM) * N + + threadCol * TN + resIdxN])[0] = + tmp; + } + } + } + } +} \ No newline at end of file diff --git a/cat_files/siboehm_CMakeLists.txt b/cat_files/siboehm_CMakeLists.txt new file mode 100644 index 00000000..7e344af8 --- /dev/null +++ b/cat_files/siboehm_CMakeLists.txt @@ -0,0 +1,36 @@ +cmake_minimum_required(VERSION 3.19) +project(NVIDIA_SGEMM_PRACTICE LANGUAGES CXX CUDA) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +find_package(CUDA REQUIRED) + +# ensure cuda is available +include(CheckLanguage) +check_language(CUDA) + +set(CMAKE_CXX_STANDARD 20) +set(CUDA_COMPUTE_CAPABILITY 86) + +# in debug mode, add debug symbols to device code +# this disables most optimizations and kills performance +add_compile_options("$<$,$>:-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}) \ No newline at end of file diff --git a/cat_files/siboehm_cuBLAS_sgemm.cu b/cat_files/siboehm_cuBLAS_sgemm.cu new file mode 100644 index 00000000..c6062184 --- /dev/null +++ b/cat_files/siboehm_cuBLAS_sgemm.cu @@ -0,0 +1,108 @@ +#include +#include +#include + +/* + * 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; +} \ No newline at end of file diff --git a/cat_files/siboehm_kernels.cuh b/cat_files/siboehm_kernels.cuh new file mode 100644 index 00000000..7691b297 --- /dev/null +++ b/cat_files/siboehm_kernels.cuh @@ -0,0 +1,14 @@ +#pragma once + +#include "kernels/10_kernel_warptiling.cuh" +#include "kernels/11_kernel_double_buffering.cuh" +#include "kernels/12_kernel_double_buffering.cuh" +#include "kernels/1_naive.cuh" +#include "kernels/2_kernel_global_mem_coalesce.cuh" +#include "kernels/3_kernel_shared_mem_blocking.cuh" +#include "kernels/4_kernel_1D_blocktiling.cuh" +#include "kernels/5_kernel_2D_blocktiling.cuh" +#include "kernels/6_kernel_vectorize.cuh" +#include "kernels/7_kernel_resolve_bank_conflicts.cuh" +#include "kernels/8_kernel_bank_extra_col.cuh" +#include "kernels/9_kernel_autotuned.cuh" \ No newline at end of file diff --git a/cat_files/siboehm_runner.cu b/cat_files/siboehm_runner.cu new file mode 100644 index 00000000..b3b731a2 --- /dev/null +++ b/cat_files/siboehm_runner.cu @@ -0,0 +1,549 @@ +#include "kernels.cuh" +#include "runner.cuh" +#include +#include +#include +#include + +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<<>>(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> + <<>>(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> + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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 + <<>>(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"); + } +} \ No newline at end of file diff --git a/cat_files/siboehm_runner.cuh b/cat_files/siboehm_runner.cuh new file mode 100644 index 00000000..3ac967e6 --- /dev/null +++ b/cat_files/siboehm_runner.cuh @@ -0,0 +1,26 @@ +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include + +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); \ No newline at end of file diff --git a/cat_files/siboehm_sgemm.cu b/cat_files/siboehm_sgemm.cu new file mode 100644 index 00000000..bb6dd38c --- /dev/null +++ b/cat_files/siboehm_sgemm.cu @@ -0,0 +1,168 @@ +#include +#include +#include +#include +#include +#include +#include + +#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 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; +}; \ No newline at end of file diff --git a/cat_files/siboehm_simplest_kernel.cu b/cat_files/siboehm_simplest_kernel.cu new file mode 100644 index 00000000..2fc15ccc --- /dev/null +++ b/cat_files/siboehm_simplest_kernel.cu @@ -0,0 +1,46 @@ +#include +#include +#include + +__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<<>>(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); +} diff --git a/cat_files/wangzyon_CMakeLists.txt b/cat_files/wangzyon_CMakeLists.txt new file mode 100644 index 00000000..64682336 --- /dev/null +++ b/cat_files/wangzyon_CMakeLists.txt @@ -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}) + diff --git a/cat_files/wangzyon_kernel.cuh b/cat_files/wangzyon_kernel.cuh new file mode 100644 index 00000000..cde4381c --- /dev/null +++ b/cat_files/wangzyon_kernel.cuh @@ -0,0 +1,9 @@ +#pragma once + +#include "kernel/kernel_1.cuh" +#include "kernel/kernel_2.cuh" +#include "kernel/kernel_3.cuh" +#include "kernel/kernel_4.cuh" +#include "kernel/kernel_5.cuh" +#include "kernel/kernel_6.cuh" +#include "kernel/kernel_7.cuh" \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_1.cuh b/cat_files/wangzyon_kernel_1.cuh new file mode 100644 index 00000000..94ce9a77 --- /dev/null +++ b/cat_files/wangzyon_kernel_1.cuh @@ -0,0 +1,19 @@ +#pragma once + +#include +#include +#include +#include + +__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]; +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_2.cuh b/cat_files/wangzyon_kernel_2.cuh new file mode 100644 index 00000000..73bde1a1 --- /dev/null +++ b/cat_files/wangzyon_kernel_2.cuh @@ -0,0 +1,45 @@ +#pragma once + +#include +#include +#include +#include + +template +__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]; +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_3.cuh b/cat_files/wangzyon_kernel_3.cuh new file mode 100644 index 00000000..f55a32ce --- /dev/null +++ b/cat_files/wangzyon_kernel_3.cuh @@ -0,0 +1,71 @@ +#pragma once + +#include +#include +#include +#include + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_4.cuh b/cat_files/wangzyon_kernel_4.cuh new file mode 100644 index 00000000..bb2214fa --- /dev/null +++ b/cat_files/wangzyon_kernel_4.cuh @@ -0,0 +1,76 @@ +#pragma once + +#include +#include +#include +#include + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_5.cuh b/cat_files/wangzyon_kernel_5.cuh new file mode 100644 index 00000000..2fcd0587 --- /dev/null +++ b/cat_files/wangzyon_kernel_5.cuh @@ -0,0 +1,88 @@ +#pragma once + +#include +#include +#include +#include + +template +__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]; + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_6.cuh b/cat_files/wangzyon_kernel_6.cuh new file mode 100644 index 00000000..f40c2501 --- /dev/null +++ b/cat_files/wangzyon_kernel_6.cuh @@ -0,0 +1,110 @@ +#pragma once + +#include +#include +#include +#include + +#define OFFSET(row, col, ld) ((row)*(ld)+(col)) +#define FETCH_FLOAT4(pointer) (reinterpret_cast(&(pointer))[0]) + +template +__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; + } + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_kernel_7.cuh b/cat_files/wangzyon_kernel_7.cuh new file mode 100644 index 00000000..53639644 --- /dev/null +++ b/cat_files/wangzyon_kernel_7.cuh @@ -0,0 +1,180 @@ +#pragma once + +#include +#include +#include +#include + +#define OFFSET(row, col, ld) ((row)*(ld)+(col)) +#define FETCH_FLOAT4(pointer) (reinterpret_cast(&(pointer))[0]) + +template +__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; + } + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_sgemm.cu b/cat_files/wangzyon_sgemm.cu new file mode 100644 index 00000000..b93c59d3 --- /dev/null +++ b/cat_files/wangzyon_sgemm.cu @@ -0,0 +1,119 @@ +#include +#include +#include +#include + +#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; +}; diff --git a/cat_files/wangzyon_utils.cu b/cat_files/wangzyon_utils.cu new file mode 100644 index 00000000..bd94ef8b --- /dev/null +++ b/cat_files/wangzyon_utils.cu @@ -0,0 +1,199 @@ +#include +#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<<>>(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><<>>(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><<>>(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><<>>(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><<>>(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><<>>(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><<>>(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><<>>(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; + } +} \ No newline at end of file diff --git a/cat_files/wangzyon_utils.cuh b/cat_files/wangzyon_utils.cuh new file mode 100644 index 00000000..fa423199 --- /dev/null +++ b/cat_files/wangzyon_utils.cuh @@ -0,0 +1,42 @@ +#pragma once +#include +#include +#include +#include +#include +#include +#include + +/* +===================================== +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); \ No newline at end of file