9ca33cf4d5
Sources (all CUDA 10.2 compatible, no CUTLASS/Triton dependency): - leimao/CUDA-GEMM-Optimization: v00-v07, fp16 WMMA variant, double buffered - siboehm/SGEMM_CUDA: kernel 1-12, warp tiling + double buffering - wangzyon/NVIDIA_SGEMM_PRACTICE: kernel 1-7 - edtallison/sgemm-cuda: kernel 1-12 (reimplementation with notes) Key porting issue: ALL kernels hardcode WARPSIZE=32. BI-V100 has warp_size=64. Need to: 1. Replace all 32U / WARPSIZE constants with 64 2. Adjust warp subtile decomposition (WMITER, WNITER, WSUBM, WSUBN) 3. Adjust shared memory bank conflict avoidance (may have different bank count) 4. Test __shfl_down_sync with mask=0xFFFFFFFFFFFFFFFF (64-bit)
98 lines
3.4 KiB
Plaintext
98 lines
3.4 KiB
Plaintext
#pragma once
|
|
|
|
#include <algorithm>
|
|
#include <cassert>
|
|
#include <cstdio>
|
|
#include <cstdlib>
|
|
#include <cublas_v2.h>
|
|
#include <cuda_runtime.h>
|
|
|
|
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
|
|
|
|
template <const int BM, const int BN, const int BK, const int TM, const int TN>
|
|
__global__ void sgemmVectorize(int M, int N, int K, float alpha, float *A,
|
|
float *B, float beta, float *C) {
|
|
const uint cRow = blockIdx.y;
|
|
const uint cCol = blockIdx.x;
|
|
|
|
// BN/TN are the number of threads to span a column
|
|
const int threadCol = threadIdx.x % (BN / TN);
|
|
const int threadRow = threadIdx.x / (BN / TN);
|
|
|
|
// allocate space for the current blocktile in smem
|
|
__shared__ float As[BM * BK];
|
|
__shared__ float Bs[BK * BN];
|
|
|
|
// Move blocktile to beginning of A's row and B's column
|
|
A += cRow * BM * K;
|
|
B += cCol * BN;
|
|
C += cRow * BM * N + cCol * BN;
|
|
|
|
// calculating the indices that this thread will load into SMEM
|
|
// we'll load 128bit / 32bit = 4 elements per thread at each step
|
|
const uint innerRowA = threadIdx.x / (BK / 4);
|
|
const uint innerColA = threadIdx.x % (BK / 4);
|
|
const uint innerRowB = threadIdx.x / (BN / 4);
|
|
const uint innerColB = threadIdx.x % (BN / 4);
|
|
|
|
// allocate thread-local cache for results in registerfile
|
|
float threadResults[TM * TN] = {0.0};
|
|
float regM[TM] = {0.0};
|
|
float regN[TN] = {0.0};
|
|
|
|
// outer-most loop over block tiles
|
|
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
|
|
// populate the SMEM caches
|
|
// transpose A while loading it
|
|
float4 tmp =
|
|
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
|
|
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
|
|
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
|
|
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
|
|
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
|
|
|
|
reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
|
|
reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
|
|
__syncthreads();
|
|
|
|
// advance blocktile
|
|
A += BK; // move BK columns to right
|
|
B += BK * N; // move BK rows down
|
|
|
|
// calculate per-thread results
|
|
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
|
|
// block into registers
|
|
for (uint i = 0; i < TM; ++i) {
|
|
regM[i] = As[dotIdx * BM + threadRow * TM + i];
|
|
}
|
|
for (uint i = 0; i < TN; ++i) {
|
|
regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
|
|
}
|
|
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
|
|
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
|
|
threadResults[resIdxM * TN + resIdxN] +=
|
|
regM[resIdxM] * regN[resIdxN];
|
|
}
|
|
}
|
|
}
|
|
__syncthreads();
|
|
}
|
|
|
|
// write out the results
|
|
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
|
|
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
|
|
// load C vector into registers
|
|
float4 tmp = reinterpret_cast<float4 *>(
|
|
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
|
|
// perform GEMM update in reg
|
|
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
|
|
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
|
|
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
|
|
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
|
|
// write back
|
|
reinterpret_cast<float4 *>(
|
|
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
|
|
tmp;
|
|
}
|
|
}
|
|
} |