diff --git a/cat_files/edtallison_01_naive.cuh b/cat_files/edtallison_01_naive.cuh deleted file mode 100644 index 04ed8e8c..00000000 --- a/cat_files/edtallison_01_naive.cuh +++ /dev/null @@ -1,37 +0,0 @@ -# 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 deleted file mode 100644 index 2f6f4504..00000000 --- a/cat_files/edtallison_02_kernel_global_mem_coalesce.cuh +++ /dev/null @@ -1,25 +0,0 @@ -#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 deleted file mode 100644 index 975cc47f..00000000 --- a/cat_files/edtallison_03_kernel_shared_mem_blocking.cuh +++ /dev/null @@ -1,54 +0,0 @@ -#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 deleted file mode 100644 index 12c9c5d7..00000000 --- a/cat_files/edtallison_04_kernel_1D_blocktiling.cuh +++ /dev/null @@ -1,80 +0,0 @@ -#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 deleted file mode 100644 index 0b361f7c..00000000 --- a/cat_files/edtallison_05_kernel_2D_blocktiling.cuh +++ /dev/null @@ -1,102 +0,0 @@ -#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 deleted file mode 100644 index 665d4e50..00000000 --- a/cat_files/edtallison_06_kernel_vectorize.cuh +++ /dev/null @@ -1,98 +0,0 @@ -#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 deleted file mode 100644 index d571008e..00000000 --- a/cat_files/edtallison_07_kernel_resolve_bank_conflicts.cuh +++ /dev/null @@ -1,103 +0,0 @@ -#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 deleted file mode 100644 index c362cba8..00000000 --- a/cat_files/edtallison_08_kernel_bank_extra_col.cuh +++ /dev/null @@ -1,103 +0,0 @@ -#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 deleted file mode 100644 index 6841952b..00000000 --- a/cat_files/edtallison_09_kernel_autotuned.cuh +++ /dev/null @@ -1,127 +0,0 @@ -#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 deleted file mode 100644 index 2cc66f36..00000000 --- a/cat_files/edtallison_10_kernel_warptiling.cuh +++ /dev/null @@ -1,187 +0,0 @@ -#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 deleted file mode 100644 index f54b61cf..00000000 --- a/cat_files/edtallison_11_kernel_double_buffering.cuh +++ /dev/null @@ -1,220 +0,0 @@ -#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 deleted file mode 100644 index c0b0a40f..00000000 --- a/cat_files/edtallison_12_kernel_double_buffering.cuh +++ /dev/null @@ -1,229 +0,0 @@ -#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 deleted file mode 100644 index fde93413..00000000 --- a/cat_files/edtallison_CMakeLists.txt +++ /dev/null @@ -1,36 +0,0 @@ -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 deleted file mode 100644 index c6062184..00000000 --- a/cat_files/edtallison_cuBLAS_sgemm.cu +++ /dev/null @@ -1,108 +0,0 @@ -#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 deleted file mode 100644 index b027b017..00000000 --- a/cat_files/edtallison_kernels.cuh +++ /dev/null @@ -1,14 +0,0 @@ -#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 deleted file mode 100644 index 4c0142dc..00000000 --- a/cat_files/edtallison_runner.cu +++ /dev/null @@ -1,549 +0,0 @@ -#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 deleted file mode 100644 index 3ac967e6..00000000 --- a/cat_files/edtallison_runner.cuh +++ /dev/null @@ -1,26 +0,0 @@ -#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 deleted file mode 100644 index bb6dd38c..00000000 --- a/cat_files/edtallison_sgemm.cu +++ /dev/null @@ -1,168 +0,0 @@ -#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 deleted file mode 100644 index 2fc15ccc..00000000 --- a/cat_files/edtallison_simplest_kernel.cu +++ /dev/null @@ -1,46 +0,0 @@ -#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 deleted file mode 100644 index 2cc66f36..00000000 --- a/cat_files/siboehm_10_kernel_warptiling.cuh +++ /dev/null @@ -1,187 +0,0 @@ -#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 deleted file mode 100644 index f54b61cf..00000000 --- a/cat_files/siboehm_11_kernel_double_buffering.cuh +++ /dev/null @@ -1,220 +0,0 @@ -#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 deleted file mode 100644 index c0b0a40f..00000000 --- a/cat_files/siboehm_12_kernel_double_buffering.cuh +++ /dev/null @@ -1,229 +0,0 @@ -#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 deleted file mode 100644 index 47038a37..00000000 --- a/cat_files/siboehm_1_naive.cuh +++ /dev/null @@ -1,29 +0,0 @@ -#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 deleted file mode 100644 index ef0a8b94..00000000 --- a/cat_files/siboehm_2_kernel_global_mem_coalesce.cuh +++ /dev/null @@ -1,24 +0,0 @@ -#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 deleted file mode 100644 index 806b6982..00000000 --- a/cat_files/siboehm_3_kernel_shared_mem_blocking.cuh +++ /dev/null @@ -1,57 +0,0 @@ -#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 deleted file mode 100644 index 12c9c5d7..00000000 --- a/cat_files/siboehm_4_kernel_1D_blocktiling.cuh +++ /dev/null @@ -1,80 +0,0 @@ -#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 deleted file mode 100644 index 0b361f7c..00000000 --- a/cat_files/siboehm_5_kernel_2D_blocktiling.cuh +++ /dev/null @@ -1,102 +0,0 @@ -#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 deleted file mode 100644 index 665d4e50..00000000 --- a/cat_files/siboehm_6_kernel_vectorize.cuh +++ /dev/null @@ -1,98 +0,0 @@ -#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 deleted file mode 100644 index d571008e..00000000 --- a/cat_files/siboehm_7_kernel_resolve_bank_conflicts.cuh +++ /dev/null @@ -1,103 +0,0 @@ -#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 deleted file mode 100644 index c362cba8..00000000 --- a/cat_files/siboehm_8_kernel_bank_extra_col.cuh +++ /dev/null @@ -1,103 +0,0 @@ -#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 deleted file mode 100644 index 6841952b..00000000 --- a/cat_files/siboehm_9_kernel_autotuned.cuh +++ /dev/null @@ -1,127 +0,0 @@ -#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 deleted file mode 100644 index 7e344af8..00000000 --- a/cat_files/siboehm_CMakeLists.txt +++ /dev/null @@ -1,36 +0,0 @@ -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 deleted file mode 100644 index c6062184..00000000 --- a/cat_files/siboehm_cuBLAS_sgemm.cu +++ /dev/null @@ -1,108 +0,0 @@ -#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 deleted file mode 100644 index 7691b297..00000000 --- a/cat_files/siboehm_kernels.cuh +++ /dev/null @@ -1,14 +0,0 @@ -#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 deleted file mode 100644 index b3b731a2..00000000 --- a/cat_files/siboehm_runner.cu +++ /dev/null @@ -1,549 +0,0 @@ -#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 deleted file mode 100644 index 3ac967e6..00000000 --- a/cat_files/siboehm_runner.cuh +++ /dev/null @@ -1,26 +0,0 @@ -#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 deleted file mode 100644 index bb6dd38c..00000000 --- a/cat_files/siboehm_sgemm.cu +++ /dev/null @@ -1,168 +0,0 @@ -#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 deleted file mode 100644 index 2fc15ccc..00000000 --- a/cat_files/siboehm_simplest_kernel.cu +++ /dev/null @@ -1,46 +0,0 @@ -#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 deleted file mode 100644 index 64682336..00000000 --- a/cat_files/wangzyon_CMakeLists.txt +++ /dev/null @@ -1,36 +0,0 @@ -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 deleted file mode 100644 index cde4381c..00000000 --- a/cat_files/wangzyon_kernel.cuh +++ /dev/null @@ -1,9 +0,0 @@ -#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 deleted file mode 100644 index 94ce9a77..00000000 --- a/cat_files/wangzyon_kernel_1.cuh +++ /dev/null @@ -1,19 +0,0 @@ -#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 deleted file mode 100644 index 73bde1a1..00000000 --- a/cat_files/wangzyon_kernel_2.cuh +++ /dev/null @@ -1,45 +0,0 @@ -#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 deleted file mode 100644 index f55a32ce..00000000 --- a/cat_files/wangzyon_kernel_3.cuh +++ /dev/null @@ -1,71 +0,0 @@ -#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 deleted file mode 100644 index bb2214fa..00000000 --- a/cat_files/wangzyon_kernel_4.cuh +++ /dev/null @@ -1,76 +0,0 @@ -#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 deleted file mode 100644 index 2fcd0587..00000000 --- a/cat_files/wangzyon_kernel_5.cuh +++ /dev/null @@ -1,88 +0,0 @@ -#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 deleted file mode 100644 index f40c2501..00000000 --- a/cat_files/wangzyon_kernel_6.cuh +++ /dev/null @@ -1,110 +0,0 @@ -#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 deleted file mode 100644 index 53639644..00000000 --- a/cat_files/wangzyon_kernel_7.cuh +++ /dev/null @@ -1,180 +0,0 @@ -#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 deleted file mode 100644 index b93c59d3..00000000 --- a/cat_files/wangzyon_sgemm.cu +++ /dev/null @@ -1,119 +0,0 @@ -#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 deleted file mode 100644 index bd94ef8b..00000000 --- a/cat_files/wangzyon_utils.cu +++ /dev/null @@ -1,199 +0,0 @@ -#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 deleted file mode 100644 index fa423199..00000000 --- a/cat_files/wangzyon_utils.cuh +++ /dev/null @@ -1,42 +0,0 @@ -#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