data: cat 3 SGEMM repos — siboehm, wangzyon, edtallison (full clone, no --depth)

Sources:
  siboehm/SGEMM_CUDA        → upstream_ref/sgemm_siboehm/     (25 files)
  wangzyon/NVIDIA_SGEMM_PRACTICE → upstream_ref/nvidia_sgemm_practice/ (23 files, filled gaps)
  edtallison/sgemm-cuda      → upstream_ref/sgemm_edtallison/  (41 files)

All files cat'd one by one from git clone (no --depth).
These are the 3 public SGEMM repos that can compile on CUDA 10.2 + CoreX ivcore10.

Key files for BI-V100 porting:
  kernel 10 (warp tiling) — already proven on device with WARPSIZE=64
  kernel 11/12 (double buffering) — next optimization target
  sgemm.cu + runner.cu — complete build+benchmark harness
  CMakeLists.txt — build system reference
This commit is contained in:
dylan
2026-08-15 06:58:07 +00:00
parent 854fb93a8e
commit 284804ac53
65 changed files with 6963 additions and 0 deletions

View File

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

View File

@@ -0,0 +1,431 @@
![](images/head.png)
![](https://img.shields.io/badge/build-passing-brightgreen) ![](https://img.shields.io/badge/ubuntu-18.04-blue) ![](https://img.shields.io/badge/cuda-10.2-blue) ![](https://img.shields.io/badge/nvidia-RTX3090-blue) ![](https://img.shields.io/badge/cmake-3.21-blue)
# 概述
面向NVIDIA GPU使用CUDA编程逐步优化矩阵乘法运算性能
| 核函数 | 描述 | GFLOPS | 自定义核函数/CUBLAS% |
| -------- | ----------------------- | -------- | ------------------------ |
| CUBLAS | 官方库函数 | 14448.69 | 基准 |
| kernel_1 | 朴素实现 | 2262.168 | 15.65657 |
| kernel_2 | 共享内存缓存 | 4216.536 | 29.18283 |
| kernel_3 | 一维Thread Tile并行优化 | 7809.629 | 54.05078 |
| kernel_4 | 二维Thread Tile并行优化 | 12251.3 | 84.79179 |
| kernel_5 | 寄存器缓存 | 12177.95 | 84.28412 |
| kernel_6 | FLOAT4向量访存 | 13161.49 | 91.09125 |
| kernel_7 | 双缓存预取 | 13634.98 | 94.36832 |
> NVIDIA GeForce RTX 3090矩阵尺寸5120
# 配置
- 编译采用 `gcc 7.5.0` under Ubuntu 18.04.5 LTS
- NVIDIA CUDA version: `CUDA 10.2`
# 目录
```
NVIDIA_SGEMM_PRACTICE # 根目录
├── images # 图片结果
│ ├── describe_kernel_1.png
│ ├── describe_kernel_x.png
│ └── kernel_x_vs_y.png
├── test # 测试结果
│ ├── test_kernel_0.txt
│ ├── test_kernel_1.txt
│ └── test_kernel_x.txt
└── src # 源文件
│ ├── kernel
│ │ ├── kernel_1.cuh # 声明和定义
│ │ ├── kernel_2.cuh
│ │ └── kernel_x.cuh
│ ├── kernel.cuh
│ ├── utils.cuh # 辅助函数
│ └── utils.cu
├── plot.py # 根据test结果绘图
├── run.sh # 运行编译后可执行文件
├── sgemm.cu # 主程序
└── CMakeLists.txt # 编译相关
```
# 运行
1. 配置NVCC编译参数
> 在CMakeLists.txt中修改`set(CUDA_NVCC_FLAGS -arch=compute_70;-code=compute_70)`
2. 配置矩阵计算最大尺寸
> 在`sgemm.cu:16`中修改`size_len`建议初次运行设置为16过大尺寸可能导致电源超负荷主机重启
3. 编译
`cd build && cmake .. && make`
4. 运行run.sh统计各个核函数计算效率结果保存在test目录
5. 计算效率折线绘图
> `python plot.py 0 1`表示绘制CUBLAS和kernel_1计算效率对比图
# 逐步优化
## kernel 1
**Naive基础版矩阵乘法实现**
将每个逻辑线程与矩阵C的每一个元素相对应每个线程负责C中一个元素的计算
![](./images/describe_kernel_1.png)
```cpp
__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];
}
```
![](./images/kernel_culas_vs_1.png)
未经过优化的矩阵乘法性能不足CUBLAS的1/10具体分析如下
- 计算访存比每次迭代需要进行一次FMA乘累加和两次全局内存读取计算访存比1/2
- 访存量访问全局内存C矩阵每个元素计算需要访问`2K`个单精度浮点数,完成全部计算需要` 2*K*M*N`
全局内存访问延迟高几百cycle同时相同位置元素被重复读取C中同一行元素计算共享A中同一行元素C中同一列元素计算共享B中同一列元素另一方面较低的计算访存比无法有效隐藏访存延迟因此访存延迟和计算访存比是导致kernel 1效率低下的原因。
## kernel 2
**利用共享内存缓存减少全局内存访存量和访存延迟**
访存延迟来自于全局内存的高延迟和全局内存的重复访问。共享内存是片上内存具有较低的访存延迟几十cycle,使用共享内存进行缓存可降低访存延迟;
![](./images/describe_kernel_2.png)
> BM和BN表示block tile的高和宽BK表示待缓存的全局内存的步长即一个block的计算需要缓存K/BK次
共享内存缓存全局内存A tile和B tile完成C block中所有元素的FMA计算不断滑动缓存区域更新block
```cpp
/*
dim3 blockDim(1024);
dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
mysgemm_v2<32><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
*/
template<const int BLOCK_SIZE>
__global__ void mysgemm_v2(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
const int BM = BLOCK_SIZE;
const int BN = BLOCK_SIZE;
const int BK = BLOCK_SIZE;
int tx = threadIdx.x % BN;
int ty = threadIdx.x / BN;
// 申请共享内存空间
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
float tmp = 0.;
for (int k = 0; k < K; k += BK) {
// 缓存A_tile和B_tile
As[ty * BK + tx] = A[ty * K + tx];
Bs[ty * BN + tx] = B[ty * N + tx];
// 同步所有线程缓存完成
__syncthreads();
A += BK;
B += BK * N;
for (int i = 0; i < BK; i++) {
tmp += As[ty * BK + i] * Bs[i * BN + tx];
}
// FMA计算需要读取缓存数据在新一轮写入缓存前进行同步确保所有线程计算完成
__syncthreads();
}
C[ty * N + tx] = alpha * tmp + beta * C[ty * N + tx];
}
```
![](./images/kernel_1_vs_2.png)
- 访存量每个block需要从global memory中读取`(K/BK)*(BM*BK+BK*BN)`个单精度浮点数整个C存在`(M/BM)*(N/BN)`个block因此完成C中所有元素计算需要读取`(M/BM)*(N/BN)*(K/BK)*(BM*BK+BK*BN)`个单精度浮点数
kernel 1受限于全局内存的访存延迟和重复访问优化前全局访存量为`2*K*M*N`,共享内存缓存优化后,访存量减少为原来的`1/2*(1/BN)*(1/BM)`,当`BN=BM=32`访存减少至1/32另一方面shared memory访存延迟远低于全局内存因此计算效率得到了一定程度的提升。
## kernel 3
**利用一维thread tile优化**
已知可以通过增加block大小BMBN进一步降低全局内存的访问量因此将BM和BN从32提升至64
> **是否能通过无限增加block size降低全局访存**
>
> 不能一方面block分块矩阵尺寸过大block数量减少这样会造成大量 SMStreaming Multiprocessor的闲置浪费另一方面BN和BM的增加需要申请更多的共享内存单线程内共享内存占用越多活跃线程束越少不利于隐藏指令延迟
因此在增加BM和BN值的同时为了减少共享内存占用一方面减小BK值降低为8
> 当增加block size时应尤其注意共享内存的消耗限制共享内存尺寸和block中线程的数量避免因资源不足无法启动核函数
![](./images/describe_kernel_3_1.png)
另一方面通过共享内存缓存减少了全局内存访存量和FMA乘累加的访存延迟但计算访存比没有得到改善每次迭代计算都需要两个访存指令和一个计算指令因此引入thread tile即一个线程负责block中多个元素的计算TM和TN分别表示thread tile的高和宽。
![](./images/describe_kernel_3_2.png)
```cpp
/*
dim3 blockDim(512);
dim3 gridDim(CEIL_DIV(M, 64), CEIL_DIV(N, 64));
mysgemm_v3<64, 64, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
*/
template<const int BM,
const int BN,
const int BK,
const int TM>
__global__ void mysgemm_v3(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
int thread_num = BM * BN / TM; // 一个线程负责block中计算TM个元素
int tx = threadIdx.x % BN;
int ty = threadIdx.x / BN * TM;
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
/*
当前线程负责搬运全局内存中第a_tile_row行第a_tile_col列元素至共享内存第a_tile_row行第a_tile_col列
a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存
若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM表示每个线程搬运一轮即可完成所需元素的搬运;
若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
*/
int a_tile_row = threadIdx.x / BK;
int a_tile_col = threadIdx.x % BK;
int a_tile_stride = thread_num / BK;
int b_tile_row = threadIdx.x / BN;
int b_tile_col = threadIdx.x % BN;
int b_tile_stride = thread_num / BN;
float tmp[TM + 1] = {0.}; // 每个线程负责TM个元素则需要申请TM个寄存器保存累加值额外的一个寄存器用于缓存
#pragma unroll
for (int k = 0; k < K; k += BK) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
}
__syncthreads();
A += BK;
B += BK * N;
#pragma unroll
for (int i = 0; i < BK; i++) {
tmp[TM] = Bs[tx + i * BN]; // 额外的一个寄存器避免反复从共享内存中读取Bs[tx + i * BN]
#pragma unroll // 循环展开,增加指令并行度
for (int j = 0; j < TM; j++) {
tmp[j] += As[(ty + j) * BK + i] * tmp[TM];
}
}
__syncthreads();
}
#pragma unroll
for (int j = 0; j < TM; j++) {
C[(ty + j) * N + tx] = alpha * tmp[j] + beta * C[(ty + j) * N + tx];
}
}
```
![](./images/kernel_2_vs_3.png)
本例从两方面进行优化:
- 全局内存访存量:相比于初始版本,通过对`64*64`block size进行缓存访存量降至1/64
- 计算访存比引入thread tile利用单个线程负责多个元素计算增加计算访存比当TM=8时每执行共享内存As的8个次访存指令和共享内存Bs的1个访存指令可执行8次计算指令相比初始版本的计算访存比1:2提高至8:9有效隐藏访存延迟
通过本例的两方面优化,矩阵乘法计算效率显著提高近一倍;
## kernel 4
**利用二维thread tile优化**
将thread tile设置为二维即一个线程负责一小块元素的计算从而进一步增加block尺寸减少全局访存数量
> 增加thread tile尺寸可以在相同的线程数量或更少的线程数量下计算更大的block size;
更重要的是单线程负责计算更多的C元素区域可以增加指令级并行程度
> 为什么可以提高指令并行程度?
>
> 单线程处理的指令数量越多,流水线级越长,由于单线程流水线可并行处理多条指令,虽然单条指令执行变慢,但单位时间内处理的指令数量变多,提高了吞吐量,隐藏指令延迟;指令级并发相比与线程级并发更具优势。
![](./images/describe_kernel_4.png)
设置一个线程负责8×8区域内元素计算即thread tile=8×8TM=8,TN=8
```cpp
// BM=BN=128BK=8TM=TN=8共享内存大小128*8
dim3 blockDim(256);
dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
mysgemm_v4<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
int a_tile_row = threadIdx.x / BK;
int a_tile_col = threadIdx.x % BK;
int a_tile_stride = thread_num / BK; // 128*8/256=4需要所有线程搬运4轮可将全局内存中128*8大小区域搬运至共享内存
int b_tile_row = threadIdx.x / BN;
int b_tile_col = threadIdx.x % BN;
int b_tile_stride = thread_num / BN;
// 每个线程负责TM*TN个元素则需要申请TM*TN个寄存器保存累加值
float tmp[TM][TN] = {0.};
// 单个线程循环TMTN完成thread tile内元素的乘累加
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];
}
```
全局访存量:相比未引入共享内存缓存版本,全局内存访存量减少至`1/2*(1/BM+1/BN)=1/128`,访存量显著降低。
![](./images/kernel_3_vs_4.png)
实际测试发现相比与一维thread tile由于二维thread tile进一步降低了全局访存量、提升计算访存比矩阵乘法效率显著提升一倍。
## kernel 5
**寄存器缓存共享内存**
![](./images/describe_kernel_5.png)
由下方代码可知单个线程计算thread tile元素乘累加时共享内存会被重复访问。
```cpp
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]; //内层循环中 As[(ty + j) * BK + i] 重复访问TN次
}
```
共享内存相比全局内存能够大大减少访存延迟但共享内存延迟几十cycle相比于计算延迟几cycle仍然较大因此采用寄存器对共享内存As、Bs进行缓存避免共享内存的重复访问
```cpp
float a_frag[TM] = {0.};
float b_frag[TN] = {0.};
for (int i = 0; i < BK; i++) {
for (int j = 0; j < TM; j++) {
a_frag[j] = As[(ty + j) * BK + i]; // 采用a_frag寄存器数组缓存thread tile所需的As共享内存数据
}
for (int l = 0; l < TN; l++) {
b_frag[l] = Bs[tx + l + i * BN]; // 采用b_frag寄存器数组缓存thread tile所需的Bs共享内存数据
}
for (int j = 0; j < TM; j++) {
for (int l = 0; l < TN; l++)
tmp[j][l] += a_frag[j] * b_frag[l];
}
}
```
当TM=TN=8时经过寄存器缓存每个thread tile需要执行8个As共享内存访存指令和8个Bs共享内存访存指令可进行8×8=64个计算指令计算访存比相比于初始版本的1/2提升至64:16可有效隐藏访存延迟
![](./images/kernel_4_vs_5.png)
实际测试发现,经寄存器缓存实际性能并未发生明显变化,原因可能是当前性能瓶颈并非共享内存的重复访问;
## kernel 6
**向量内存指令FLOAT4优化**
- 计算指令GPU是以4维向量为基本单位进行计算的4个浮点数组成的float4向量是GPU最基本的类型使用GPU对两个float4进行向量计算与对两个整数或两个浮点数进行计算一样只需要一个指令即可完成
- 内存指令:与发出单个指令生成单独的内存事务获取相同数量的字节相比,通过向量内存指令所需的内存事务更少,减少了内存控制器的争用;另一方面,使用矢量加载每个字节需要更少的索引计算;
![](./images/describe_kernel_6.png)
例如BM=128BK=8线程数量为256若每个线程每次取1个浮点数每个线程需要消耗4次内存指令才能将全局内存搬运至共享内存若采用float4向量内存指令每个线程每次可以搬运4个浮点数则每个线程仅需要执行一次内存指令即可完成搬运。
关键代码示例如下:
```cpp
#define OFFSET(row, col, ld) ((row)*(ld)+(col))
#define FETCH_FLOAT4(pointer) (reinterpret_cast<float4*>(&(pointer))[0])
float ldg_a_reg[4 * ldg_a_num] = {0.}; // 每个线程搬运ldg_a_num轮寄存器缓存ldg_a_num个float4元素用于转置As矩阵
// 共享内存缓存全局内存
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];
}
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)]); // 不需要转置
}
// 寄存器缓存共享内存
// ty,tx为当前线程对应thread tile的左上角元素在block中的位置
#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
}
```
全局内存无法直接写入共享内存需要寄存器做中介其中As写入将全局内存->将寄存器->共享内存过程显示的描述出来而Bs写入并不是不需要寄存器参与只是编译器隐藏了这段代码As缓存显示运用寄存器的目的在于将As进行转置转置前的一列在转置后变成一行内存连续便于float4读取
![kernel_1](./images/kernel_5_vs_6.png)
实际测试,整体计算效率增加;
## kernel 7
**数据预取**
单缓存是指申请单块共享内存,缓存全局数据,申请单块寄存器内存,缓存共享数据,单块缓存不能实现读取和存储并行进行,因为数据之间存在依赖。例如单缓存场景,计算依赖共享内存数据,为保证计算前全局内存完全存入共享内存,需要进行一次同步;同样因为计算依赖共享内存数据,所以在存新一轮全局内存到共享内存前也需要进行一次同步,保证上一轮计算完成。
双缓存通过申请双倍存储空间,将读和写分开,计算数据读取一块存储空间同时,可以同时向另一块内存写入下一轮依赖的数据,因此,只需要保证计算前待读取共享内存完成写入,即一次同步即可。
> 双缓存使读写同步进行,实现数据预取,隐藏内存延迟。
![](./images/describe_kernel_7.png)
![](./images/kernel_6_vs_7.png)
采用双缓存技术实现数据预取,计算效率得到了进一步提升;
![](./images/kernel_culas_vs_7.png)
基本可以接近CUBLAS官方矩阵乘法的计算效率

View File

@@ -0,0 +1,66 @@
import os
import re
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
import argparse
def parse_file(file):
with open(file, 'r') as f:
lines = [line.strip() for line in f.readlines()]
data = []
pattern = "Average elasped time: \((.*?)\) second, performance: \((.*?)\) GFLOPS. size: \((.*?)\)."
for line in lines:
r = re.match(pattern, line)
if r:
gflops = float(r.group(2))
data.append(gflops)
return data
def plot(num1, num2, y1, y2, save_dir):
x = [(i + 1) * 256 for i in range(len(y1))]
fig = plt.figure(figsize=(12, 10))
if num1 == 0:
num1 = "culas"
plt.plot(x, y1, c='k', linewidth=2, label=f"kernel_{num1}")
plt.plot(x, y2, c='b', linewidth=2, label=f"kernel_{num2}")
plt.legend()
plt.scatter(x, y1, marker="s", s=60, c='', edgecolors='k', linewidth=2)
plt.scatter(x, y2, marker="^", s=60, c='', edgecolors='b', linewidth=2)
plt.tick_params(labelsize=10)
plt.xlabel("Matrix size (M=N=K)", fontsize=12, fontweight='bold')
plt.ylabel("Performance (GFLOPS)", fontsize=12, fontweight='bold')
plt.title(f"Comparison bewteen: kernel_{num1} and kernel_{num2}", fontsize=16, fontweight='bold')
x_major_locator = MultipleLocator(256)
plt.gca().xaxis.set_major_locator(x_major_locator)
plt.savefig(f"{save_dir}/kernel_{num1}_vs_{num2}.png")
def main(args):
root = os.path.dirname(os.path.abspath(__file__))
data1 = parse_file(os.path.join(root, f'test/test_kernel_{args.one}.txt'))
data2 = parse_file(os.path.join(root, f'test/test_kernel_{args.another}.txt'))
plot(args.one, args.another, data1, data2, args.save_dir)
def parse_args():
parser = argparse.ArgumentParser(description='plot kernel performance')
parser.add_argument('one', type=int, help='one kernel num')
parser.add_argument('another', type=int, help='another kernel num')
parser.add_argument('--save_dir', default='images')
return parser.parse_args()
if __name__ == "__main__":
args = parse_args()
main(args)
# python plot.py 0 1

View File

@@ -0,0 +1,16 @@
#!/bin/bash
# 全部kernel运行
rm ./test/test_kernel*
echo -n "test_kernel:"
for((i=0;i<=7;i++))
do
echo -n "${i}..."
file_name="./test/test_kernel_${i}.txt"
./sgemm ${i} >> ${file_name}
done
# 单个kernel运行
# kernel_num=$1
# file_name="test_kernel_${kernel_num}.txt"
# ./sgemm ${kernel_num} | tee ./test/${file_name}

View File

@@ -0,0 +1,9 @@
#pragma once
#include "kernel/kernel_1.cuh"
#include "kernel/kernel_2.cuh"
#include "kernel/kernel_3.cuh"
#include "kernel/kernel_4.cuh"
#include "kernel/kernel_5.cuh"
#include "kernel/kernel_6.cuh"
#include "kernel/kernel_7.cuh"

View File

@@ -0,0 +1,19 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
__global__ __launch_bounds__(1024) void
mysgemm_v1(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int gx = blockIdx.x * blockDim.x + threadIdx.x; // 全局x
int gy = blockIdx.y * blockDim.y + threadIdx.y; // 全局y
float tmp = 0.;
for (int i = 0; i < K; i++) {
tmp += A[gy * K + i] * B[i * N + gx]; // 两次全局内存访问和一次FMA累加乘
}
C[gy * N + gx] = alpha * tmp + beta * C[gy * N + gx];
}

View File

@@ -0,0 +1,45 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
template<const int BLOCK_SIZE>
__global__ void mysgemm_v2(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
const int BM = BLOCK_SIZE;
const int BN = BLOCK_SIZE;
const int BK = BLOCK_SIZE;
int tx = threadIdx.x % BN;
int ty = threadIdx.x / BN;
// 申请共享内存空间
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
float tmp = 0.;
for (int k = 0; k < K; k += BK) {
// 缓存A_tile和B_tile
As[ty * BK + tx] = A[ty * K + tx];
Bs[ty * BN + tx] = B[ty * N + tx];
// 同步所有线程缓存完成
__syncthreads();
A += BK;
B += BK * N;
for (int i = 0; i < BK; i++) {
tmp += As[ty * BK + i] * Bs[i * BN + tx];
}
// FMA计算需要读取缓存数据在新一轮写入缓存前进行同步确保所有线程计算完成
__syncthreads();
}
C[ty * N + tx] = alpha * tmp + beta * C[ty * N + tx];
}

View File

@@ -0,0 +1,71 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
template<const int BM,
const int BN,
const int BK,
const int TM>
__global__ void mysgemm_v3(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
int thread_num = BM * BN / TM; // 一个线程负责block中计算TM个元素
int tx = threadIdx.x % BN;
int ty = threadIdx.x / BN * TM;
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
/*
当前线程负责搬运全局内存中第a_tile_row行第a_tile_col列元素至共享内存第a_tile_row行第a_tile_col列
a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存
若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM表示每个线程搬运一轮即可完成所需元素的搬运;
若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
*/
int a_tile_row = threadIdx.x / BK;
int a_tile_col = threadIdx.x % BK;
int a_tile_stride = thread_num / BK;
int b_tile_row = threadIdx.x / BN;
int b_tile_col = threadIdx.x % BN;
int b_tile_stride = thread_num / BN;
float tmp[TM + 1] = {0.}; // 每个线程负责TM个元素则需要申请TM个寄存器保存累加值额外的一个寄存器用于缓存
#pragma unroll
for (int k = 0; k < K; k += BK) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
}
__syncthreads();
A += BK;
B += BK * N;
#pragma unroll
for (int i = 0; i < BK; i++) {
tmp[TM] = Bs[tx + i * BN]; // 额外的一个寄存器避免反复从共享内存中读取Bs[tx + i * BN]
#pragma unroll // 循环展开,增加指令并行度
for (int j = 0; j < TM; j++) {
tmp[j] += As[(ty + j) * BK + i] * tmp[TM];
}
}
__syncthreads();
}
#pragma unroll
for (int j = 0; j < TM; j++) {
C[(ty + j) * N + tx] = alpha * tmp[j] + beta * C[(ty + j) * N + tx];
}
}

View File

@@ -0,0 +1,76 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
template<const int BM,
const int BN,
const int BK,
const int TM,
const int TN>
__global__ void mysgemm_v4(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
int block_row_thread = BN / TN;
int block_col_thread = BM / TM;
int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
int tx = (threadIdx.x % block_row_thread) * TN;
int ty = (threadIdx.x / block_row_thread) * TM;
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
/*
当前线程负责搬运全局内存中第a_tile_row行第a_tile_col列元素至共享内存第a_tile_row行第a_tile_col列
a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存
若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM表示每个线程搬运一轮即可完成所需元素的搬运;
若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
*/
int a_tile_row = threadIdx.x / BK;
int a_tile_col = threadIdx.x % BK;
int a_tile_stride = thread_num / BK;
int b_tile_row = threadIdx.x / BN;
int b_tile_col = threadIdx.x % BN;
int b_tile_stride = thread_num / BN;
float tmp[TM][TN] = {0.}; // 每个线程负责TM*TN个元素则需要申请TM*TN个寄存器保存累加值额外的一个寄存器用于缓存
#pragma unroll
for (int k = 0; k < K; k += BK) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
}
__syncthreads();
A += BK;
B += BK * N;
#pragma unroll
for (int i = 0; i < BK; i++) {
#pragma unroll // 循环展开,增加指令并行度
for (int j = 0; j < TM; j++) {
for (int l = 0; l < TN; l++)
tmp[j][l] += As[(ty + j) * BK + i] * Bs[tx + l + i * BN];
}
}
__syncthreads();
}
#pragma unroll
for (int j = 0; j < TM; j++) {
for (int l = 0; l < TN; l++)
C[(ty + j) * N + tx + l] = alpha * tmp[j][l] + beta * C[(ty + j) * N + tx + l];
}
}

View File

@@ -0,0 +1,88 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
template<const int BM,
const int BN,
const int BK,
const int TM,
const int TN>
__global__ void mysgemm_v5(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
int block_row_thread = BN / TN;
int block_col_thread = BM / TM;
int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
int tx = (threadIdx.x % block_row_thread) * TN;
int ty = (threadIdx.x / block_row_thread) * TM;
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
/*
当前线程负责搬运全局内存中第a_tile_row行第a_tile_col列元素至共享内存第a_tile_row行第a_tile_col列
a_tile_stride表示block中线程可搬运a_tile_stride行至共享内存
若BM=64,BK=8,thread_num=512,则a_tile_stride=64,a_tile_stride=BM表示每个线程搬运一轮即可完成所需元素的搬运;
若BM=128,BK=8,thread_num=512,则a_tile_stride=64,表示每个线程搬运两轮即可完成所需元素的搬运;
*/
int a_tile_row = threadIdx.x / BK;
int a_tile_col = threadIdx.x % BK;
int a_tile_stride = thread_num / BK;
int b_tile_row = threadIdx.x / BN;
int b_tile_col = threadIdx.x % BN;
int b_tile_stride = thread_num / BN;
float tmp[TM][TN] = {0.}; // 每个线程负责TM*TN个元素则需要申请TM*TN个寄存器保存累加值额外的一个寄存器用于缓存
float a_frag[TM] = {0.};
float b_frag[TN] = {0.};
#pragma unroll
for (int k = 0; k < K; k += BK) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
As[(a_tile_row + i) * BK + a_tile_col] = A[(a_tile_row + i) * K + a_tile_col];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
Bs[(b_tile_row + i) * BN + b_tile_col] = B[(b_tile_row + i) * N + b_tile_col];
}
__syncthreads();
A += BK;
B += BK * N;
#pragma unroll
for (int i = 0; i < BK; i++) {
#pragma unroll
for (int j = 0; j < TM; j++) {
a_frag[j] = As[(ty + j) * BK + i];
}
#pragma unroll
for (int l = 0; l < TN; l++) {
b_frag[l] = Bs[tx + l + i * BN];
}
#pragma unroll
for (int j = 0; j < TM; j++) {
#pragma unroll
for (int l = 0; l < TN; l++)
tmp[j][l] += a_frag[j] * b_frag[l];
}
}
__syncthreads();
}
#pragma unroll
for (int j = 0; j < TM; j++) {
for (int l = 0; l < TN; l++)
C[(ty + j) * N + tx + l] = alpha * tmp[j][l] + beta * C[(ty + j) * N + tx + l];
}
}

View File

@@ -0,0 +1,110 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
#define OFFSET(row, col, ld) ((row)*(ld)+(col))
#define FETCH_FLOAT4(pointer) (reinterpret_cast<float4*>(&(pointer))[0])
template<const int BM,
const int BN,
const int BK,
const int TM,
const int TN>
__global__ void mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
const int block_row_thread = BN / TN;
const int block_col_thread = BM / TM;
const int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
// 当前线程对应thread tile的左上角元素在block中的位置
int tx = (threadIdx.x % block_row_thread) * TN;
int ty = (threadIdx.x / block_row_thread) * TM;
__shared__ float As[BK * BM];
__shared__ float Bs[BK * BN];
const int ldg_a_num = BK * BM / thread_num / 4; // 每个线程搬运4个浮点数完成搬运至As需要所有线程搬运ldg_a_num轮
const int ldg_b_num = BK * BN / thread_num / 4; // 每个线程搬运4个浮点数完成搬运至Bs需要所有线程搬运ldg_b_num轮
int a_tile_row = threadIdx.x / (BK / 4); // 每行4个字节作为一个内存块当前线程负责第a_tile_row行的第a_tile_col个内存块的搬运
int a_tile_col = threadIdx.x % (BK / 4) * 4;
int a_tile_stride = BM / ldg_a_num; // 一共BM行搬运ldg_a_num轮每论搬运a_tile_stride行
int b_tile_row = threadIdx.x / (BN / 4); // 每行4个字节作为一个内存块当前线程负责第b_tile_row行的第b_tile_col个内存块的搬运
int b_tile_col = threadIdx.x % (BN / 4) * 4;
int b_tile_stride = BK / ldg_b_num; // 一共BK行搬运ldg_b_num轮每论搬运b_tile_stride行
float accum[TM][TN] = {0.}; // 每个线程负责TM*TN个元素则需要申请TM*TN个寄存器保存累加值额外的一个寄存器用于缓存
// 计算ldg_a_num的所有参数必须全部是const否则不能用来申明数组大小
float ldg_a_reg[4 * ldg_a_num] = {0.}; // 每个线程搬运ldg_a_num轮寄存器缓存ldg_a_num个float4元素用于转置As矩阵
float a_frag[TM]; // 缓存As共享内存
float b_frag[TN]; // 缓存Bs共享内存
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
#pragma unroll
for (int k = 0; k < K; k += BK) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
int ldg_index = i / a_tile_stride * 4; // 第ldg_index轮
FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
FETCH_FLOAT4(A[OFFSET(a_tile_row + i, a_tile_col, K)]);
// As转置存其中ldg_a_reg做中间缓存目的是读取时可以按FLOAT4读取
As[OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
As[OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
As[OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
As[OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
FETCH_FLOAT4(Bs[OFFSET(b_tile_row + i, b_tile_col, BN)]) =
FETCH_FLOAT4(B[OFFSET(b_tile_row + i, b_tile_col, N)]); // 不需要转置
}
__syncthreads();
A += BK;
B += BK * N;
#pragma unroll
for (int i = 0; i < BK; i++) {
#pragma unroll
for (int m = 0; m < TM; m += 4) {
FETCH_FLOAT4(a_frag[m]) = FETCH_FLOAT4(As[OFFSET(i, ty + m, BM)]); // 偏移到当前thread tile
}
#pragma unroll
for (int n = 0; n < TN; n += 4) {
FETCH_FLOAT4(b_frag[n]) = FETCH_FLOAT4(Bs[OFFSET(i, tx + n, BN)]); // 偏移到当前thread tile
}
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n++) {
accum[m][n] += a_frag[m] * b_frag[n];
}
}
}
__syncthreads();
}
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n += 4) {
float4 ctmp = FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]);
//float4 atmp = FETCH_FLOAT4(accum[m][n]);
ctmp.x = alpha * accum[m][n] + beta * ctmp.x;
ctmp.y = alpha * accum[m][n + 1] + beta * ctmp.y;
ctmp.z = alpha * accum[m][n + 2] + beta * ctmp.z;
ctmp.w = alpha * accum[m][n + 3] + beta * ctmp.w;
FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]) = ctmp;
}
}
}

View File

@@ -0,0 +1,180 @@
#pragma once
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
#define OFFSET(row, col, ld) ((row)*(ld)+(col))
#define FETCH_FLOAT4(pointer) (reinterpret_cast<float4*>(&(pointer))[0])
template<const int BM,
const int BN,
const int BK,
const int TM,
const int TN>
__global__ void mysgemm_v7(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
int bx = blockIdx.x;
int by = blockIdx.y;
const int block_row_thread = BN / TN;
const int block_col_thread = BM / TM;
const int thread_num = block_row_thread * block_col_thread; // 一个线程负责计算block中TM*TN个元素
// 当前线程对应thread tile的左上角元素在block中的位置
int tx = (threadIdx.x % block_row_thread) * TN;
int ty = (threadIdx.x / block_row_thread) * TM;
__shared__ float As[2][BK * BM]; // 增加一倍共享内存大小用于缓存
__shared__ float Bs[2][BK * BN];
const int ldg_a_num = BK * BM / thread_num / 4; // 每个线程搬运4个浮点数完成搬运至As需要所有线程搬运ldg_a_num轮
const int ldg_b_num = BK * BN / thread_num / 4; // 每个线程搬运4个浮点数完成搬运至Bs需要所有线程搬运ldg_b_num轮
int a_tile_row = threadIdx.x / (BK / 4); // 每行4个字节作为一个内存块当前线程负责第a_tile_row行的第a_tile_col个内存块的搬运
int a_tile_col = threadIdx.x % (BK / 4) * 4;
int a_tile_stride = BM / ldg_a_num; // 一共BM行搬运ldg_a_num轮每论搬运a_tile_stride行
int b_tile_row = threadIdx.x / (BN / 4); // 每行4个字节作为一个内存块当前线程负责第b_tile_row行的第b_tile_col个内存块的搬运
int b_tile_col = threadIdx.x % (BN / 4) * 4;
int b_tile_stride = BK / ldg_b_num; // 一共BK行搬运ldg_b_num轮每论搬运b_tile_stride行
float accum[TM][TN] = {0.}; // 每个线程负责TM*TN个元素则需要申请TM*TN个寄存器保存累加值额外的一个寄存器用于缓存
// 计算ldg_a_num的所有参数必须全部是const否则不能用来申明数组大小
float ldg_a_reg[4 * ldg_a_num] = {0.}; // 每个线程搬运ldg_a_num轮寄存器缓存ldg_a_num个float4元素用于转置As矩阵
float ldg_b_reg[4 * ldg_b_num] = {0.}; // 每个线程搬运ldg_a_num轮寄存器缓存ldg_a_num个float4元素用于转置As矩阵
float a_frag[2][TM]; // 缓存As共享内存,增加一倍寄存器大小用于缓存
float b_frag[2][TN]; // 缓存Bs共享内存,增加一倍寄存器大小用于缓存
// 移动到当前block
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * N + bx * BN];
// first global to shared
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
int ldg_index = i / a_tile_stride * 4; // 第ldg_index轮
FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
FETCH_FLOAT4(A[OFFSET(a_tile_row + i, a_tile_col, K)]);
// As转置存其中ldg_a_reg做中间缓存目的是读取时可以按FLOAT4读取
As[0][OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
As[0][OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
As[0][OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
As[0][OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
FETCH_FLOAT4(Bs[0][OFFSET(b_tile_row + i, b_tile_col, BN)]) =
FETCH_FLOAT4(B[OFFSET(b_tile_row + i, b_tile_col, N)]); // 不需要转置
}
__syncthreads();
// first shared to frag
#pragma unroll
for (int m = 0; m < TM; m += 4) {
FETCH_FLOAT4(a_frag[0][m]) = FETCH_FLOAT4(As[0][OFFSET(0, ty + m, BM)]); // 偏移到当前thread tile
}
#pragma unroll
for (int n = 0; n < TN; n += 4) {
FETCH_FLOAT4(b_frag[0][n]) = FETCH_FLOAT4(Bs[0][OFFSET(0, tx + n, BN)]); // 偏移到当前thread tile
}
int write_index = 1;
int load_index;
int k = 0;
do {
k += BK;
// load global to reg
if (k < K) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
int ldg_index = i / a_tile_stride * 4; // 第ldg_index轮
FETCH_FLOAT4(ldg_a_reg[ldg_index]) =
FETCH_FLOAT4(A[OFFSET(a_tile_row + i, k + a_tile_col, K)]);
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
int ldg_index = i / b_tile_stride * 4; // 第ldg_index轮
FETCH_FLOAT4(ldg_b_reg[ldg_index]) =
FETCH_FLOAT4(B[OFFSET(k + b_tile_row + i, b_tile_col, N)]);
}
}
load_index = write_index ^ 1;
#pragma unroll
for (int bk = 0; bk < BK - 1; bk++) {
for (int m = 0; m < TM; m += 4) {
FETCH_FLOAT4(a_frag[(bk + 1) % 2][m]) = FETCH_FLOAT4(
As[load_index][OFFSET(bk + 1, ty + m, BM)]); // 偏移到当前thread tile
}
#pragma unroll
for (int n = 0; n < TN; n += 4) {
FETCH_FLOAT4(b_frag[(bk + 1) % 2][n]) = FETCH_FLOAT4(
Bs[load_index][OFFSET(bk + 1, tx + n, BN)]); // 偏移到当前thread tile
}
#pragma unroll
for (int m = 0; m < TM; m++) {
for (int n = 0; n < TN; n++) {
accum[m][n] += a_frag[bk % 2][m] * b_frag[bk % 2][n];
}
}
}
if (k < K) {
#pragma unroll
for (int i = 0; i < BM; i += a_tile_stride) {
int ldg_index = i / a_tile_stride * 4;
As[write_index][OFFSET(a_tile_col, i + a_tile_row, BM)] = ldg_a_reg[ldg_index];
As[write_index][OFFSET(a_tile_col + 1, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 1];
As[write_index][OFFSET(a_tile_col + 2, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 2];
As[write_index][OFFSET(a_tile_col + 3, i + a_tile_row, BM)] = ldg_a_reg[ldg_index + 3];
}
#pragma unroll
for (int i = 0; i < BK; i += b_tile_stride) {
int ldg_index = i / b_tile_stride * 4;
FETCH_FLOAT4(Bs[write_index][OFFSET(b_tile_row + i, b_tile_col, BN)]) =
FETCH_FLOAT4(ldg_b_reg[ldg_index]);
}
__syncthreads();
#pragma unroll
for (int m = 0; m < TM; m += 4) {
FETCH_FLOAT4(a_frag[0][m]) = FETCH_FLOAT4(
As[write_index][OFFSET(0, ty + m, BM)]); // 偏移到当前thread tile
}
#pragma unroll
for (int n = 0; n < TN; n += 4) {
FETCH_FLOAT4(b_frag[0][n]) = FETCH_FLOAT4(
Bs[write_index][OFFSET(0, tx + n, BN)]); // 偏移到当前thread tile
}
write_index ^= 1;
}
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n++) {
accum[m][n] += a_frag[(BK - 1) % 2][m] * b_frag[(BK - 1) % 2][n];
}
}
} while (k < K);
// C = alpha*AB+C
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n += 4) {
float4 ctmp = FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]);
ctmp.x = alpha * accum[m][n] + beta * ctmp.x;
ctmp.y = alpha * accum[m][n + 1] + beta * ctmp.y;
ctmp.z = alpha * accum[m][n + 2] + beta * ctmp.z;
ctmp.w = alpha * accum[m][n + 3] + beta * ctmp.w;
FETCH_FLOAT4(C[OFFSET(ty + m, tx + n, N)]) = ctmp;
}
}
}

View File

@@ -0,0 +1,199 @@
#include <stdio.h>
#include "utils.cuh"
#include "kernel.cuh"
float get_sec() {
struct timeval time;
gettimeofday(&time, NULL);
return (1e6 * time.tv_sec + time.tv_usec);
}
float cpu_elapsed_time(float &beg, float &end) {
return 1.0e-6 * (end - beg);
}
void cudaCheck(cudaError_t error, const char *file, int line) {
if (error != cudaSuccess) {
printf("[CUDA ERROR] at file %s(line %d):\n%s\n", file, line, cudaGetErrorString(error));
exit(EXIT_FAILURE);
}
return;
};
void CudaDeviceInfo() {
int deviceId;
cudaGetDevice(&deviceId);
cudaDeviceProp props;
cudaGetDeviceProperties(&props, deviceId);
/*
* There should be no need to modify the output string below.
*/
printf("Device ID: %d\n\
*Number of SMs: %d\n\
Compute Capability Major: %d\n\
Compute Capability Minor: %d\n\
memoryBusWidth: %d\n\
*maxThreadsPerBlock: %d\n\
maxThreadsPerMultiProcessor: %d\n\
*totalGlobalMem: %zuM\n\
sharedMemPerBlock: %zuKB\n\
*sharedMemPerMultiprocessor: %zuKB\n\
totalConstMem: %zuKB\n\
*multiProcessorCount: %d\n\
*Warp Size: %d\n",
deviceId,
props.multiProcessorCount,
props.major,
props.minor,
props.memoryBusWidth,
props.maxThreadsPerBlock,
props.maxThreadsPerMultiProcessor,
props.totalGlobalMem / 1024 / 1024,
props.sharedMemPerBlock / 1024,
props.sharedMemPerMultiprocessor / 1024,
props.totalConstMem / 1024,
props.multiProcessorCount,
props.warpSize);
};
void randomize_matrix(float *mat, int N) {
// NOTICE: 使用gettimeofdays替代srand((unsigned)time(NULL));time精度过低产生相同随机数
struct timeval time;
gettimeofday(&time, NULL);
srand(time.tv_usec);
for (int i = 0; i < N; i++) {
float tmp = (float) (rand() % 5) + 0.01 * (rand() % 5);
tmp = (rand() % 2 == 0) ? tmp : tmp * (-1.);
mat[i] = tmp;
}
}
void copy_matrix(float *src, float *dest, int N) {
int i;
for (i = 0; src + i && dest + i && i < N; i++)
*(dest + i) = *(src + i);
if (i != N)
printf("copy failed at %d while there are %d elements in total.\n", i, N);
}
void print_matrix(const float *A, int M, int N) {
int i;
printf("[");
for (i = 0; i < M * N; i++) {
if ((i + 1) % N == 0)
printf("%5.2f ", A[i]);
else
printf("%5.2f, ", A[i]);
if ((i + 1) % N == 0) {
if (i + 1 < M * N)
printf(";\n");
}
}
printf("]\n");
}
bool verify_matrix(float *mat1, float *mat2, int N) {
double diff = 0.0;
int i;
for (i = 0; mat1 + i && mat2 + i && i < N; i++) {
diff = fabs((double) mat1[i] - (double) mat2[i]);
if (diff > 1e-2) {
printf("error. %5.2f,%5.2f,%d\n", mat1[i], mat2[i], i);
return false;
}
}
return true;
}
#define CEIL_DIV(M, N) ((M) + (N)-1) / (N)
void test_cublas(cublasHandle_t handle, int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
//cublas列主序计算https://www.cnblogs.com/cuancuancuanhao/p/7763256.html
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B, N, A, K, &beta, C, N);
}
void test_mysgemm_v1(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(32, 32);
dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
mysgemm_v1<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_mysgemm_v2(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(1024);
dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));
mysgemm_v2<32><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_mysgemm_v3(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(512);
dim3 gridDim(CEIL_DIV(M, 64), CEIL_DIV(N, 64));
mysgemm_v3<64, 64, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_mysgemm_v4(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(256);
dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
mysgemm_v4<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_mysgemm_v5(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(256);
dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
mysgemm_v5<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
//void test_mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
// dim3 blockDim(4);
// dim3 gridDim(CEIL_DIV(M, 8), CEIL_DIV(N, 8));
// mysgemm_v6<8, 8, 4, 4, 4><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
//}
void test_mysgemm_v6(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(256);
dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
mysgemm_v6<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_mysgemm_v7(int M, int N, int K, float alpha, float *A, float *B, float beta, float *C) {
dim3 blockDim(256);
dim3 gridDim(CEIL_DIV(M, 128), CEIL_DIV(N, 128));
mysgemm_v7<128, 128, 8, 8, 8><<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
}
void test_kernel(int kernel_num, int M, int N, int K, float alpha, float *A, float *B, float beta, float *C,
cublasHandle_t handle) {
switch (kernel_num) {
case 0:
test_cublas(handle, M, N, K, alpha, A, B, beta, C);
break;
case 1:
test_mysgemm_v1(M, N, K, alpha, A, B, beta, C);
break;
case 2:
test_mysgemm_v2(M, N, K, alpha, A, B, beta, C);
break;
case 3:
test_mysgemm_v3(M, N, K, alpha, A, B, beta, C);
break;
case 4:
test_mysgemm_v4(M, N, K, alpha, A, B, beta, C);
break;
case 5:
test_mysgemm_v5(M, N, K, alpha, A, B, beta, C);
break;
case 6:
test_mysgemm_v6(M, N, K, alpha, A, B, beta, C);
break;
case 7:
test_mysgemm_v7(M, N, K, alpha, A, B, beta, C);
break;
default:
break;
}
}

View File

@@ -0,0 +1,42 @@
#pragma once
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#include <sys/time.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
/*
=====================================
CUDA操作
=====================================
*/
void cudaCheck(cudaError_t error, const char *file, int line); //CUDA错误检查
void CudaDeviceInfo(); // 打印CUDA信息
/*
=====================================
矩阵操作
=====================================
*/
void randomize_matrix(float *mat, int N); // 随机初始化矩阵
void copy_matrix(float *src, float *dest, int N); // 复制矩阵
void print_matrix(const float *A, int M, int N); // 打印矩阵
bool verify_matrix(float *mat1, float *mat2, int N); // 验证矩阵
/*
=====================================
计时操作
=====================================
*/
float get_current_sec(); // 获取当前时刻
float cpu_elapsed_time(float &beg, float &end); // 计算时间差
/*
=====================================
kernel操作
=====================================
*/
//调用指定核函数计算矩阵乘法
void test_kernel(int kernel_num, int m, int n, int k, float alpha, float *A, float *B, float beta, float *C, cublasHandle_t handle);

View File

@@ -0,0 +1,36 @@
cmake_minimum_required(VERSION 3.19)
project(NVIDIA_SGEMM_PRACTICE LANGUAGES CXX CUDA)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
find_package(CUDA REQUIRED)
# ensure cuda is available
include(CheckLanguage)
check_language(CUDA)
set(CMAKE_CXX_STANDARD 20)
set(CUDA_COMPUTE_CAPABILITY 75)
# in debug mode, add debug symbols to device code
# this disables most optimizations and kills performance
add_compile_options("$<$<AND:$<CONFIG:Debug>,$<COMPILE_LANGUAGE:CUDA>>:-G;-src-in-ptx>")
# add_compile_options("--ptxas-options=-v")
# Configure header file search paths
include_directories(${CUDA_INCLUDE_DIRS})
include_directories(${PROJECT_SOURCE_DIR}/src)
# Configure the source file path to be compiled
aux_source_directory(${PROJECT_SOURCE_DIR}/src SRC)
# generate executable
add_executable(sgemm sgemm.cu ${SRC})
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
add_executable(cuBLAS_sgemm cuBLAS_sgemm.cu )
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(cuBLAS_sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
add_executable(simplest_kernel simplest_kernel.cu)
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(simplest_kernel ${CUDA_LIBRARIES})

View File

@@ -0,0 +1,34 @@
.PHONY: all build debug clean profile bench cuobjdump
CMAKE := cmake
BUILD_DIR := build
BENCHMARK_DIR := benchmark_results
all: build
build:
@mkdir -p $(BUILD_DIR)
@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Release ..
@$(MAKE) -C $(BUILD_DIR)
debug:
@mkdir -p $(BUILD_DIR)
@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Debug ..
@$(MAKE) -C $(BUILD_DIR)
clean:
@rm -rf $(BUILD_DIR)
FUNCTION := $$(cuobjdump -symbols build/sgemm | grep -i Warptiling | awk '{print $$NF}')
cuobjdump: build
@cuobjdump -arch sm_86 -sass -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.sass
@cuobjdump -arch sm_86 -ptx -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.ptx
# Usage: make profile KERNEL=<integer> PREFIX=<optional string>
profile: build
@ncu --set full --export $(BENCHMARK_DIR)/$(PREFIX)kernel_$(KERNEL) --force-overwrite $(BUILD_DIR)/sgemm $(KERNEL)
bench: build
@bash gen_benchmark_results.sh

View File

@@ -0,0 +1,78 @@
Reimplementation of Simon Boehm's [CUDA SGEMM](https://github.com/siboehm/SGEMM_CUDA) kernels.
Following the [article](https://siboehm.com/articles/22/CUDA-MMM), for my learning :).
## Run on Google Colab
[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/edtallison/sgemm-cuda/blob/master/run_on_colab.ipynb)
Click the link above to open and run the project in a GPU-enabled Google Colab environment. No additional setup required.
# Notes
(also scattered throughout kernel code)
## 1. Naive
- **Three-level hierarchy of computation**
- Grid, block, thread. Assume grid and block are 2D, thread is the atomic unit of computation.
- Blocks can have up to 1024 threads.
- Threads within the same block share memory (SMEM).
- **Grid and block indexing**
- `gridDim` specifies dimensions of the grid i.e. rows and columns of blocks.
- `blockDim` specifies dimensions of the block i.e. rows and columns of threads.
- `blockIdx.x/y/z` specifies the block's position in the grid.
- `threadIdx.x/y/z `specifies the thread's position in the block.
- When used within a kernel, these vars are automatically assigned by the CUDA runtime.
- **Matrix Multiplication**
- Matrix multiplication: element ij of C is the dot product of row i of A and column j of B.
- In this kernel, each thread computes one element of C. This can obviously be done in parallel so no synchronisation is required.
- **Kernel Launch**
- When the kernel is launched, we make the grid as big as necessary to cover all of C, depending on the block size.
- The kernel execution is launched asynchronously i.e. the function call on the host (CPU) returns immediately.
- **Memory Access Pattern**
- Threads within the same block e.g. `threadIds` (0, 0) and (0, 1) use the same column of B.
- They each load the whole column from global memory. Hmmm this seems inefficient...
## 2. Global Memory Coalescing
- **Warps**
- In execution, within a block, threads are grouped into "warps" of 32 threads.
- Each streaming multiprocessor (SM) has four warp schedulers - physical cores that execute instructions.
- Each warp is assigned to a warp scheduler, based on a consecutive `threadId` (x, y, z).
- Threads with neighbouring `threadId` become part of the same warp.
- **Global Memory Coalescing**
- Sequential memory acceses by threads in the same warp can be grouped and executed as one.
- Important to keep in mind when optimising GMEM memory access.
- For coalescing, the memory addresses need to be consecutive, but the within-warp accesses don't need to be consecutive.
- GPU supports 32B, 64B, and 128B memory accesses.
- **Memory Access Pattern** (this part took me some time to get my head around)
- In naive kernel, iterating threads with `threadIdx.x` (which aligns with consecutive `threadId`) actually leads to consecutive threads operating on consecutive rows of A, and the same row of B
- If, instead, the threads operated on the same row of A but consecutive columns of B, this accessing of the B values could be coalesced.
- This is achieved simply by changing the x and y position indices of the C element computed by each thread.
- Note that in either case, we can use within-warp broadcasting as the same row of A or col of B is being accessed by the threads.
## 3. Shared Memory Cache-Blocking
- **SMEM in GPU Memory Architecture**
- GPU has global memory GMEM.
- Each Streaming Multiprocessor (SM) has a much smaller memory called shared memory (SMEM).
- This SMEM is partitioned among the blocks.
- Each block of threads runs on a single SM. Multiple blocks can be assigned to the same SM.
- A thread can communicate with the other threads in its block via the SMEM chunk.
- SMEM, being located on-chip, has much lower latency and higher bandwidth than GMEM.
- **Kernel Memory Access**
- Load a chunk of A and a chunk of B from GMEM into SMEM.
- Perform as much work as possible on the chunks.
- Perform partial sums on C, moving the chunks along the columns of A (same row) and rows of B (same col) until result fully computed.
- I.e. in this kernel, each block of threads computes one `BLOCKSIZE*BLOCKSIZE` tile of C.
- **Improvement**
- For this kernel, resources mostly spent in waiting for SMEM accesses to return.
- Need to make the kernel issue less SMEM instructions to improve efficiency.

View File

@@ -0,0 +1,108 @@
#include <cstdio>
#include <cublas_v2.h>
#include <cuda_runtime.h>
/*
* A stand-alone script to invoke & benchmark standard cuBLAS SGEMM performance
*/
int main(int argc, char *argv[]) {
int m = 2;
int k = 3;
int n = 4;
int print = 1;
cudaError_t cudaStat; // cudaMalloc status
cublasStatus_t stat; // cuBLAS functions status
cublasHandle_t handle; // cuBLAS context
int i, j;
float *a, *b, *c;
// malloc for a,b,c...
a = (float *)malloc(m * k * sizeof(float));
b = (float *)malloc(k * n * sizeof(float));
c = (float *)malloc(m * n * sizeof(float));
int ind = 11;
for (j = 0; j < m * k; j++) {
a[j] = (float)ind++;
}
ind = 11;
for (j = 0; j < k * n; j++) {
b[j] = (float)ind++;
}
ind = 11;
for (j = 0; j < m * n; j++) {
c[j] = (float)ind++;
}
// DEVICE
float *d_a, *d_b, *d_c;
// cudaMalloc for d_a, d_b, d_c...
cudaMalloc((void **)&d_a, m * k * sizeof(float));
cudaMalloc((void **)&d_b, k * n * sizeof(float));
cudaMalloc((void **)&d_c, m * n * sizeof(float));
stat = cublasCreate(&handle); // initialize CUBLAS context
cudaMemcpy(d_a, a, m * k * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_b, b, k * n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_c, c, m * n * sizeof(float), cudaMemcpyHostToDevice);
float alpha = 1.0f;
float beta = 0.5f;
if (print == 1) {
printf("alpha = %4.0f, beta = %4.0f\n", alpha, beta);
printf("A = (mxk: %d x %d)\n", m, k);
for (i = 0; i < m; i++) {
for (j = 0; j < k; j++) {
printf("%4.1f ", a[i * m + j]);
}
printf("\n");
}
printf("B = (kxn: %d x %d)\n", k, n);
for (i = 0; i < k; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", b[i * n + j]);
}
printf("\n");
}
printf("C = (mxn: %d x %d)\n", m, n);
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", c[i * n + j]);
}
printf("\n");
}
}
stat = cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, d_b, n,
d_a, k, &beta, d_c, n);
cudaMemcpy(c, d_c, m * n * sizeof(float), cudaMemcpyDeviceToHost);
if (print == 1) {
printf("\nC after SGEMM = \n");
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", c[i * n + j]);
}
printf("\n");
}
}
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
cublasDestroy(handle); // destroy CUBLAS context
free(a);
free(b);
free(c);
return EXIT_SUCCESS;
}

View File

@@ -0,0 +1,71 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "b9326784",
"metadata": {},
"source": [
"This notebook installs dependencies, builds the project, and runs a selected kernel on Colab's GPU.\n",
"\n",
"**Note:** Ensure Colab runtime type is set to GPU (Runtime → Change runtime type → GPU)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35ef02c6",
"metadata": {},
"outputs": [],
"source": [
"# install system dependencies\n",
"!apt-get update && apt-get install -y cmake ninja-build"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "37a03120",
"metadata": {},
"outputs": [],
"source": [
"# clone\n",
"!git clone https://github.com/edtallison/sgemm-cuda.git\n",
"%cd sgemm-cuda"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cd1855de",
"metadata": {},
"outputs": [],
"source": [
"# build\n",
"!mkdir -p build && cd build && cmake -G Ninja .. && ninja"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "093a0b6d",
"metadata": {},
"outputs": [],
"source": [
"# run for kernel 1\n",
"!cd build && ./sgemm 1"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

View File

@@ -0,0 +1,33 @@
banks_naive = lambda r, c: (r * 32 + c) % 32
banks_one_extra = lambda r, c: (r * 33 + c) % 32
ITEMS_PER_WARP = 8
def printBankConflicts(bank_fun):
for c in range(1):
banks = []
for i in range(32):
row = (i * ITEMS_PER_WARP) // 16
col = (i * ITEMS_PER_WARP + c) % 16
banks.append((i, row, col, bank_fun(row, col)))
print("Step", c, "\n", "\n".join(["(" + ",".join(str(x) for x in i) + ")" for i in banks]))
d = {k: 0 for k in range(32)}
for i in banks:
d[i[-1]] += 1
count = 0
for key, val in d.items():
if val > 0:
count += 1
print(
f"Bank conflicts (Step {c}): {sorted(d.items(), key=lambda item: item[1], reverse=True)[0][1]}, banks accessed: {count}/32\n"
)
print("---NAIVE---")
printBankConflicts(banks_naive, 32)
print("\n---EXTRA COL---")
printBankConflicts(banks_one_extra, 33)

View File

@@ -0,0 +1,115 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
WM_VALUES=(32 64 128 256)
WN_VALUES=(32 64 128 256)
WNITER_VALUES=(1 2 4 8)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
NUM_THREADS_VALUES=(128 256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
OUTPUT="../benchmark_results/kernel_10_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="0"
WARPSIZE=32
TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for BK in "${BK_VALUES[@]}"; do
for BM in "${BM_VALUES[@]}"; do
for BN in "${BN_VALUES[@]}"; do
for WM in "${WM_VALUES[@]}"; do
for WN in "${WN_VALUES[@]}"; do
for WN_ITER in "${WNITER_VALUES[@]}"; do
for TM in "${TM_VALUES[@]}"; do
for TN in "${TN_VALUES[@]}"; do
for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
echo ""
CONFIG_NUM=$(( CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
NUM_WARPS=$(( NUM_THREADS / 32 ))
if ! (( BN % WN == 0 && BM % WM == 0 )); then
echo "Error: BN % WN must be 0 and BM % WM must be 0."
continue
fi
if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
continue
fi
if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
continue
fi
WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BK == 0 )); then
echo "Error: (NUM_THREADS * 4) % BK must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BN == 0 )); then
echo "Error: (NUM_THREADS * 4) % BN must be 0."
continue
fi
if ! (( BN % (16 * TN) == 0 )); then
echo "Error: BN must be a multiple of 16 * TN."
continue
fi
if ! (( BM % (16 * TM) == 0 )); then
echo "Error: BM must be a multiple of 16 * TM."
continue
fi
if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K10_NUM_THREADS = .*/const uint K10_NUM_THREADS = $NUM_THREADS;/" $RUNNER
sed -i "s/const uint K10_BN = .*/const uint K10_BN = $BN;/" $RUNNER
sed -i "s/const uint K10_BM = .*/const uint K10_BM = $BM;/" $RUNNER
sed -i "s/const uint K10_BK = .*/const uint K10_BK = $BK;/" $RUNNER
sed -i "s/const uint K10_WM = .*/const uint K10_WM = $WM;/" $RUNNER
sed -i "s/const uint K10_WN = .*/const uint K10_WN = $WN;/" $RUNNER
sed -i "s/const uint K10_WNITER = .*/const uint K10_WNITER = $WN_ITER;/" $RUNNER
sed -i "s/const uint K10_TM = .*/const uint K10_TM = $TM;/" $RUNNER
sed -i "s/const uint K10_TN = .*/const uint K10_TN = $TN;/" $RUNNER
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 4 seconds if it doesn't finish
timeout -v 8 ./sgemm 10 | tee -a $OUTPUT
done
done
done
done
done
done
done
done
done

View File

@@ -0,0 +1,115 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
WM_VALUES=(32 64 128 256)
WN_VALUES=(32 64 128 256)
WNITER_VALUES=(1 2 4 8)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
NUM_THREADS_VALUES=(128 256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
OUTPUT="../benchmark_results/kernel_11_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="0"
WARPSIZE=32
TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for BK in "${BK_VALUES[@]}"; do
for BM in "${BM_VALUES[@]}"; do
for BN in "${BN_VALUES[@]}"; do
for WM in "${WM_VALUES[@]}"; do
for WN in "${WN_VALUES[@]}"; do
for WN_ITER in "${WNITER_VALUES[@]}"; do
for TM in "${TM_VALUES[@]}"; do
for TN in "${TN_VALUES[@]}"; do
for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
echo ""
CONFIG_NUM=$(( CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
NUM_WARPS=$(( NUM_THREADS / 32 ))
if ! (( BN % WN == 0 && BM % WM == 0 )); then
echo "Error: BN % WN must be 0 and BM % WM must be 0."
continue
fi
if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
continue
fi
if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
continue
fi
WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BK == 0 )); then
echo "Error: (NUM_THREADS * 4) % BK must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BN == 0 )); then
echo "Error: (NUM_THREADS * 4) % BN must be 0."
continue
fi
if ! (( BN % (16 * TN) == 0 )); then
echo "Error: BN must be a multiple of 16 * TN."
continue
fi
if ! (( BM % (16 * TM) == 0 )); then
echo "Error: BM must be a multiple of 16 * TM."
continue
fi
if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K11_NUM_THREADS = .*/const uint K11_NUM_THREADS = $NUM_THREADS;/" $RUNNER
sed -i "s/const uint K11_BN = .*/const uint K11_BN = $BN;/" $RUNNER
sed -i "s/const uint K11_BM = .*/const uint K11_BM = $BM;/" $RUNNER
sed -i "s/const uint K11_BK = .*/const uint K11_BK = $BK;/" $RUNNER
sed -i "s/const uint K11_WM = .*/const uint K11_WM = $WM;/" $RUNNER
sed -i "s/const uint K11_WN = .*/const uint K11_WN = $WN;/" $RUNNER
sed -i "s/const uint K11_WNITER = .*/const uint K11_WNITER = $WN_ITER;/" $RUNNER
sed -i "s/const uint K11_TM = .*/const uint K11_TM = $TM;/" $RUNNER
sed -i "s/const uint K11_TN = .*/const uint K11_TN = $TN;/" $RUNNER
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 8 seconds if it doesn't finish
timeout -v 8 ./sgemm 11 | tee -a $OUTPUT
done
done
done
done
done
done
done
done
done

View File

@@ -0,0 +1,86 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
NUM_THREADS_VALUES=(256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
KERNEL="../src/kernels/9_kernel_autotuned.cuh"
OUTPUT="../benchmark_results/kernel_9_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="2"
TOTAL_CONFIGS="$(( ${#NUM_THREADS_VALUES[@]} * ${#BK_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for bk in ${BK_VALUES[@]}; do
for tm in ${TM_VALUES[@]}; do
for tn in ${TN_VALUES[@]}; do
for bm in ${BM_VALUES[@]}; do
for bn in ${BN_VALUES[@]}; do
for nt in ${NUM_THREADS_VALUES[@]}; do
echo ""
CONFIG_NUM=$(( $CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
config="BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NT=$nt"
if [[ $(( ($nt * 4) % bk )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BK = $(( ($nt * 4) % bk )) != 0))"
continue
fi
if [[ $(( ($nt * 4) % bn )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BN = $(( ($nt * 4) % bn )) != 0))"
continue
fi
if [[ $(( $bn % (16 * $tn ) )) -ne 0 ]]; then
echo "QUANTIZATION: Skipping $config because BN % (16 * TN) = $(( $bn % (16 * $tn ) )) != 0))"
continue
fi
if [[ $(( $bm % (16 * $tm ) )) -ne 0 ]]; then
echo "QUANTIZATION: Skipping $config because BM % (16 * TM) = $(( $bm % (16 * $tm ) )) != 0))"
continue
fi
if [[ $(( ($bm * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (BM * BK) % (4 * NUM_THREADS) = $(( ($bm * $bk) % ( 4 * 256 ) )) != 0))"
continue
fi
if [[ $(( ($bn * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (BN * BK) % (4 * NUM_THREADS) = $(( ($bn * $bk) % ( 4 * 256 ) )) != 0))"
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K9_BK = .*/const uint K9_BK = $bk;/" $RUNNER
sed -i "s/const uint K9_TM = .*/const uint K9_TM = $tm;/" $RUNNER
sed -i "s/const uint K9_TN = .*/const uint K9_TN = $tn;/" $RUNNER
sed -i "s/const uint K9_BM = .*/const uint K9_BM = $bm;/" $RUNNER
sed -i "s/const uint K9_BN = .*/const uint K9_BN = $bn;/" $RUNNER
sed -i "s/const int K9_NUM_THREADS = .*/const int K9_NUM_THREADS = $nt;/" $KERNEL
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NUM_THREADS=$nt" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 4 seconds if it doesn't finish
timeout -v 4 ./sgemm 9 | tee -a $OUTPUT
done
done
done
done
done
done

View File

@@ -0,0 +1,168 @@
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <runner.cuh>
#include <vector>
#define cudaCheck(err) (cudaCheck(err, __FILE__, __LINE__))
const std::string errLogFile = "matrixValidationFailure.txt";
int main(int argc, char **argv) {
if (argc != 2) {
std::cerr << "Please select a kernel (range 0 - 12, 0 for NVIDIA cuBLAS)"
<< std::endl;
exit(EXIT_FAILURE);
}
// get kernel number
int kernel_num = std::stoi(argv[1]);
if (kernel_num < 0 || kernel_num > 12) {
std::cerr << "Please enter a valid kernel number (0-12)" << std::endl;
exit(EXIT_FAILURE);
}
// get environment variable for device
int deviceIdx = 0;
if (getenv("DEVICE") != NULL) {
deviceIdx = atoi(getenv("DEVICE"));
}
cudaCheck(cudaSetDevice(deviceIdx));
printf("Running kernel %d on device %d.\n", kernel_num, deviceIdx);
// print some device info
// CudaDeviceInfo();
// Declare the handle, create the handle, cublasCreate will return a value of
// type cublasStatus_t to determine whether the handle was created
// successfully (the value is 0)
cublasHandle_t handle;
if (cublasCreate(&handle)) {
std::cerr << "Create cublas handle error." << std::endl;
exit(EXIT_FAILURE);
};
// Using cudaEvent for gpu stream timing, cudaEvent is equivalent to
// publishing event tasks in the target stream
float elapsed_time;
cudaEvent_t beg, end;
cudaEventCreate(&beg);
cudaEventCreate(&end);
// cuBLAS FLOPs ceiling is reached at 8192
std::vector<int> SIZE = {128, 256, 512, 1024, 2048, 4096};
long m, n, k, max_size;
max_size = SIZE[SIZE.size() - 1];
std::cout << "Max size: " << max_size << std::endl;
float alpha = 0.5, beta = 3.0; // GEMM input parameters, C=α*AB+β*C
float *A = nullptr, *B = nullptr, *C = nullptr,
*C_ref = nullptr; // host matrices
float *dA = nullptr, *dB = nullptr, *dC = nullptr,
*dC_ref = nullptr; // device matrices
A = (float *)malloc(sizeof(float) * max_size * max_size);
B = (float *)malloc(sizeof(float) * max_size * max_size);
C = (float *)malloc(sizeof(float) * max_size * max_size);
C_ref = (float *)malloc(sizeof(float) * max_size * max_size);
randomize_matrix(A, max_size * max_size);
randomize_matrix(B, max_size * max_size);
randomize_matrix(C, max_size * max_size);
cudaCheck(cudaMalloc((void **)&dA, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dB, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dC, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dC_ref, sizeof(float) * max_size * max_size));
cudaCheck(cudaMemcpy(dA, A, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dB, B, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dC, C, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dC_ref, C, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
int repeat_times = 50;
for (int size : SIZE) {
m = n = k = size;
std::cout << "dimensions(m=n=k) " << m << ", alpha: " << alpha
<< ", beta: " << beta << std::endl;
// Verify the correctness of the calculation, and execute it once before the
// kernel function timing to avoid cold start errors
if (kernel_num != 0) {
run_kernel(0, m, n, k, alpha, dA, dB, beta, dC_ref,
handle); // cuBLAS
run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC,
handle); // Executes the kernel, modifies the result matrix
cudaCheck(cudaDeviceSynchronize());
cudaCheck(cudaGetLastError()); // Check for async errors during kernel run
cudaMemcpy(C, dC, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
cudaMemcpy(C_ref, dC_ref, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
if (!verify_matrix(C_ref, C, m * n)) {
std::cout
<< "Failed to pass the correctness verification against NVIDIA "
"cuBLAS."
<< std::endl;
if (m <= 128) {
std::cout << " Logging faulty output into " << errLogFile << "\n";
std::ofstream fs;
fs.open(errLogFile);
fs << "A:\n";
print_matrix(A, m, n, fs);
fs << "B:\n";
print_matrix(B, m, n, fs);
fs << "C:\n";
print_matrix(C, m, n, fs);
fs << "Should:\n";
print_matrix(C_ref, m, n, fs);
}
exit(EXIT_FAILURE);
}
}
cudaEventRecord(beg);
for (int j = 0; j < repeat_times; j++) {
// We don't reset dC between runs to save time
run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC, handle);
}
cudaEventRecord(end);
cudaEventSynchronize(beg);
cudaEventSynchronize(end);
cudaEventElapsedTime(&elapsed_time, beg, end);
elapsed_time /= 1000.; // Convert to seconds
long flops = 2 * m * n * k;
printf(
"Average elapsed time: (%7.6f) s, performance: (%7.1f) GFLOPS. size: "
"(%ld).\n",
elapsed_time / repeat_times,
(repeat_times * flops * 1e-9) / elapsed_time, m);
fflush(stdout);
// make dC and dC_ref equal again (we modified dC while calling our kernel
// for benchmarking)
cudaCheck(cudaMemcpy(dC, dC_ref, sizeof(float) * m * n,
cudaMemcpyDeviceToDevice));
}
// Free up CPU and GPU space
free(A);
free(B);
free(C);
free(C_ref);
cudaFree(dA);
cudaFree(dB);
cudaFree(dC);
cudaFree(dC_ref);
cublasDestroy(handle);
return 0;
};

View File

@@ -0,0 +1,46 @@
#include <cuda_runtime.h>
#include <iostream>
#include <vector>
__global__ void kernel(uint *A, uint *B, int row) {
auto x = threadIdx.x / 4;
auto y = threadIdx.x % 4;
A[x * row + y] = x;
B[x * row + y] = y;
}
int main(int argc, char **argv) {
uint *Xs, *Ys;
uint *Xs_d, *Ys_d;
uint SIZE = 4;
Xs = (uint *)malloc(SIZE * SIZE * sizeof(uint));
Ys = (uint *)malloc(SIZE * SIZE * sizeof(uint));
cudaMalloc((void **)&Xs_d, SIZE * SIZE * sizeof(uint));
cudaMalloc((void **)&Ys_d, SIZE * SIZE * sizeof(uint));
dim3 grid_size(1, 1, 1);
dim3 block_size(4 * 4);
kernel<<<grid_size, block_size>>>(Xs_d, Ys_d, 4);
cudaMemcpy(Xs, Xs_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
cudaMemcpy(Ys, Ys_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
cudaDeviceSynchronize();
for (int row = 0; row < SIZE; ++row) {
for (int col = 0; col < SIZE; ++col) {
std::cout << "[" << Xs[row * SIZE + col] << "|" << Ys[row * SIZE + col]
<< "] ";
}
std::cout << "\n";
}
cudaFree(Xs_d);
cudaFree(Ys_d);
free(Xs);
free(Ys);
}

View File

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

View File

@@ -0,0 +1,37 @@
# pragma once
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
/*
Matrix sizes:
MxK * KxN = MxN
*/
__global__ void sgemm_naive(
int M, int N, int K, // sizes
float alpha, const float *A, const float *B, float beta, float *C // pointers used to point to matrices
) {
// compute position in C that this thread is responsible for
// "which block" * "width of block" to get to start of block + "which thread"
const uint x = blockIdx.x * blockDim.x + threadIdx.x; // "which row?" (inverted from graphical intuition, confusingly)
const uint y = blockIdx.y * blockDim.y + threadIdx.y; // "which column?"
// if M or N are not multiples of 32, there will be "extra"/"remainder" threads on the last block in x/y.
// we don't want those leftover threads to do anything (tile quantisation)
if (x < M && y < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) { // K is the size of the row in A, col in B i.e. the dot product
// A: x * K gives the start of relevant row, i enumerates across the row (col by col)
// B: y gives the relevant column, i * N enumerates down the column, (row by row)
tmp += A[x * K + i] * B[i * N + y];
}
// C = alpha*(A@B) + beta*C
// x * N takes to start of relevant row, y moves across to the relevant column
C[x * N + y] = alpha * tmp + beta * C[x * N + y];
}
}

View File

@@ -0,0 +1,25 @@
#pragma once
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
template <const uint BLOCKSIZE>
// __global__ is used to specify that the function is run on GPU, called by host (CPU)
__global__ void sgemm_global_mem_coalesce(int M, int N, int K, float alpha,
const float *A, const float *B,
float beta, float *C) {
const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); // note that blockDim is now 1-dimensional
const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);
// if statement is necessary to make things work under tile quantization
if (cRow < M && cCol < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) {
tmp += A[cRow * K + i] * B[i * N + cCol];
}
C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol];
}
}

View File

@@ -0,0 +1,54 @@
#pragma once
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BLOCKSIZE>
__global__ void sgemm_shared_mem_block(int M, int N, int K, float alpha,
const float *A, const float *B,
float beta, float *C) {
// output C block we want to compute with this threadBlock
const uint cRow = blockIdx.x;
const uint cCol = blockIdx.y;
// allocate buffer for current block in fast SMEM (shared between all threads in block)
__shared__ float As[BLOCKSIZE * BLOCKSIZE];
__shared__ float Bs[BLOCKSIZE * BLOCKSIZE];
// the inner row and col that we are accessing in this specific thread
const uint threadRow = threadIdx.x / BLOCKSIZE; // note similarity to previous kernel
const uint threadCol = threadIdx.x % BLOCKSIZE;
// advance pointers to the starting positions (they are input as pointers to first elements in the matrices)
A += cRow * BLOCKSIZE * K; // row=cRow, col=0 (the start of the relevant row)
B += cCol * BLOCKSIZE; // row=0, col=cCol (top of relevant col)
C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE; // row=cRow, col=cCol
float tmp = 0.0;
for (int bkIdx=0; bkIdx < K; bkIdx+=BLOCKSIZE) { // shifting the whole block along the row of A and col of B
// have each thread load one of the elements in A and B
// make the threadCol (=threadIdx.x) the consecutive index
// to allow GMEM access coalescing
As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
// ensure cache is fully populated
__syncthreads();
A += BLOCKSIZE; // for next iteration
B += BLOCKSIZE * N;
// execute the dotproduct on the currently cached block
for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx) {
tmp += As[threadRow * BLOCKSIZE + dotIdx] *
Bs[dotIdx * BLOCKSIZE + threadCol];
}
// sync so faster threads don't fetch the next block into cache
_syncthreads();
}
C[threadRow * N + threadCol] = alpha * tmp + beta * C[threadRow * N + threadCol];
}

View File

@@ -0,0 +1,80 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM>
__global__ void sgemm1DBlocktiling(int M, int N, int K, float alpha,
const float *A, const float *B, float beta,
float *C) {
// If we flip x and y here we get ~30% less performance for large matrices.
// The current, 30% faster configuration ensures that blocks with sequential
// blockIDs access columns of B sequentially, while sharing the same row of A.
// The slower configuration would share columns of A, but access into B would
// be non-sequential. So the faster configuration has better spatial locality
// and hence a greater L2 hit rate.
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// each warp will calculate 32*TM elements, with 32 being the columnar dim.
const int threadCol = threadIdx.x % BN;
const int threadRow = threadIdx.x / BN;
// allocate space for the current blocktile in SMEM
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// todo: adjust this to each thread to load multiple entries and
// better exploit the cache sizes
assert(BM * BK == blockDim.x);
assert(BN * BK == blockDim.x);
const uint innerColA = threadIdx.x % BK; // warp-level GMEM coalescing
const uint innerRowA = threadIdx.x / BK;
const uint innerColB = threadIdx.x % BN; // warp-level GMEM coalescing
const uint innerRowB = threadIdx.x / BN;
// allocate thread-local cache for results in registerfile
float threadResults[TM] = {0.0};
// outer loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
__syncthreads();
// advance blocktile
A += BK;
B += BK * N;
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// we make the dotproduct loop the outside loop, which facilitates
// reuse of the Bs entry, which we can cache in a tmp var.
float tmpB = Bs[dotIdx * BN + threadCol];
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
threadResults[resIdx] +=
As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
}
}
__syncthreads();
}
// write out the results
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
C[(threadRow * TM + resIdx) * N + threadCol] =
alpha * threadResults[resIdx] +
beta * C[(threadRow * TM + resIdx) * N + threadCol];
}
}

View File

@@ -0,0 +1,102 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void __launch_bounds__((BM * BN) / (TM * TN), 1)
sgemm2DBlocktiling(int M, int N, int K, float alpha, const float *A,
const float *B, float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
const uint totalResultsBlocktile = BM * BN;
// A thread is responsible for calculating TM*TN elements in the blocktile
const uint numThreadsBlocktile = totalResultsBlocktile / (TM * TN);
// ResultsPerBlock / ResultsPerThread == ThreadsPerBlock
assert(numThreadsBlocktile == blockDim.x);
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
const uint innerRowA = threadIdx.x / BK;
const uint innerColA = threadIdx.x % BK;
// calculates the number of rows of As that are being loaded in a single step
// by a single block
const uint strideA = numThreadsBlocktile / BK;
const uint innerRowB = threadIdx.x / BN;
const uint innerColB = threadIdx.x % BN;
// for both As and Bs we want each load to span the full column-width, for
// better GMEM coalescing (as opposed to spanning full row-width and iterating
// across columns)
const uint strideB = numThreadsBlocktile / BN;
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
// register caches for As and Bs
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
As[(innerRowA + loadOffset) * BK + innerColA] =
A[(innerRowA + loadOffset) * K + innerColA];
}
for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
Bs[(innerRowB + loadOffset) * BN + innerColB] =
B[(innerRowB + loadOffset) * N + innerColB];
}
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN] =
alpha * threadResults[resIdxM * TN + resIdxN] +
beta * C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN];
}
}
}

View File

@@ -0,0 +1,98 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmVectorize(int M, int N, int K, float alpha, float *A,
float *B, float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,103 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmResolveBankConflicts(int M, int N, int K, float alpha,
float *A, float *B, float beta,
float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
// "linearize" Bs while storing it
tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 0) * 16 + innerColB / 2] = tmp.x;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 1) * 16 + innerColB / 2] = tmp.y;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 2) * 16 + innerColB / 2] = tmp.z;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 3) * 16 + innerColB / 2] = tmp.w;
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[(dotIdx * 8 + i) * 16 + threadCol];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,103 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmResolveBankExtraCol(int M, int N, int K, float alpha,
float *A, float *B, float beta,
float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
const int extraCols = 5;
__shared__ float Bs[BK * (BN + extraCols)];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 0] = tmp.x;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 1] = tmp.y;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 2] = tmp.z;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 3] = tmp.w;
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * (BN + extraCols) + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,127 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
const int K9_NUM_THREADS = 256;
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void __launch_bounds__(K9_NUM_THREADS)
sgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// size of warptile
constexpr int WM = TM * 16;
constexpr int WN = TN * 16;
// iterations of warptile
constexpr int WMITER = CEIL_DIV(BM, WM);
constexpr int WNITER = CEIL_DIV(BN, WN);
// Placement of the thread in the warptile
const int threadCol = threadIdx.x % (WN / TN);
const int threadRow = threadIdx.x / (WN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
constexpr uint rowStrideA = (K9_NUM_THREADS * 4) / BK;
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
constexpr uint rowStrideB = K9_NUM_THREADS / (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[WMITER * WNITER * TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
float4 tmp = reinterpret_cast<float4 *>(
&A[(innerRowA + offset) * K + innerColA * 4])[0];
// transpose A while storing it
As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
}
for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
reinterpret_cast<float4 *>(
&Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
reinterpret_cast<float4 *>(
&B[(innerRowB + offset) * N + innerColB * 4])[0];
}
__syncthreads();
for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + (wmIdx * WM) + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + (wnIdx * WN) + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[(wmIdx * TM + resIdxM) * (WNITER * TN) +
wnIdx * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
}
}
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
}
// write out the results
for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
float *C_interim = C + (wmIdx * WM * N) + (wnIdx * WN);
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C_interim[(threadRow * TM + resIdxM) * N + threadCol * TN +
resIdxN])[0];
// perform GEMM update in reg
const int i =
(wmIdx * TM + resIdxM) * (WNITER * TN) + wnIdx * TN + resIdxN;
tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(&C_interim[(threadRow * TM + resIdxM) * N +
threadCol * TN + resIdxN])[0] =
tmp;
}
}
}
}
}

View File

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

View File

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

View File

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

View File

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

View File

@@ -0,0 +1,26 @@
#pragma once
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
void cudaCheck(cudaError_t error, const char *file,
int line); // CUDA error check
void CudaDeviceInfo(); // print CUDA information
void range_init_matrix(float *mat, int N);
void randomize_matrix(float *mat, int N);
void zero_init_matrix(float *mat, int N);
void copy_matrix(const float *src, float *dest, int N);
void print_matrix(const float *A, int M, int N, std::ofstream &fs);
bool verify_matrix(float *mat1, float *mat2, int N);
float get_current_sec(); // Get the current moment
float cpu_elapsed_time(float &beg, float &end); // Calculate time difference
void run_kernel(int kernel_num, int m, int n, int k, float alpha, float *A,
float *B, float beta, float *C, cublasHandle_t handle);

View File

@@ -0,0 +1,36 @@
cmake_minimum_required(VERSION 3.19)
project(NVIDIA_SGEMM_PRACTICE LANGUAGES CXX CUDA)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
find_package(CUDA REQUIRED)
# ensure cuda is available
include(CheckLanguage)
check_language(CUDA)
set(CMAKE_CXX_STANDARD 20)
set(CUDA_COMPUTE_CAPABILITY 86)
# in debug mode, add debug symbols to device code
# this disables most optimizations and kills performance
add_compile_options("$<$<AND:$<CONFIG:Debug>,$<COMPILE_LANGUAGE:CUDA>>:-G;-src-in-ptx>")
# add_compile_options("--ptxas-options=-v")
# Configure header file search paths
include_directories(${CUDA_INCLUDE_DIRS})
include_directories(${PROJECT_SOURCE_DIR}/src)
# Configure the source file path to be compiled
aux_source_directory(${PROJECT_SOURCE_DIR}/src SRC)
# generate executable
add_executable(sgemm sgemm.cu ${SRC})
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
add_executable(cuBLAS_sgemm cuBLAS_sgemm.cu )
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(cuBLAS_sgemm ${CUDA_LIBRARIES} ${CUDA_CUBLAS_LIBRARIES})
add_executable(simplest_kernel simplest_kernel.cu)
set_target_properties(sgemm PROPERTIES CUDA_ARCHITECTURES ${CUDA_COMPUTE_CAPABILITY})
target_link_libraries(simplest_kernel ${CUDA_LIBRARIES})

View File

@@ -0,0 +1,34 @@
.PHONY: all build debug clean profile bench cuobjdump
CMAKE := cmake
BUILD_DIR := build
BENCHMARK_DIR := benchmark_results
all: build
build:
@mkdir -p $(BUILD_DIR)
@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Release ..
@$(MAKE) -C $(BUILD_DIR)
debug:
@mkdir -p $(BUILD_DIR)
@cd $(BUILD_DIR) && $(CMAKE) -DCMAKE_BUILD_TYPE=Debug ..
@$(MAKE) -C $(BUILD_DIR)
clean:
@rm -rf $(BUILD_DIR)
FUNCTION := $$(cuobjdump -symbols build/sgemm | grep -i Warptiling | awk '{print $$NF}')
cuobjdump: build
@cuobjdump -arch sm_86 -sass -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.sass
@cuobjdump -arch sm_86 -ptx -fun $(FUNCTION) build/sgemm | c++filt > build/cuobjdump.ptx
# Usage: make profile KERNEL=<integer> PREFIX=<optional string>
profile: build
@ncu --set full --export $(BENCHMARK_DIR)/$(PREFIX)kernel_$(KERNEL) --force-overwrite $(BUILD_DIR)/sgemm $(KERNEL)
bench: build
@bash gen_benchmark_results.sh

View File

@@ -0,0 +1,42 @@
# Fast CUDA SGEMM from Scratch
Step-by-step optimization of matrix multiplication, implemented in CUDA.
For an explanation of each kernel, see [siboehm.com/CUDA-MMM](https://siboehm.com/articles/22/CUDA-MMM).
## Overview
Running the kernels on a NVIDIA A6000 (Ampere):
![](benchmark_results.png)
GFLOPs at matrix size 4096x4096:
<!-- benchmark_results -->
| Kernel | GFLOPs/s | Performance relative to cuBLAS |
|:------------------------------------|----------:|:-------------------------------|
| 1: Naive | `309.0` | 1.3% |
| 2: GMEM Coalescing | `1986.5` | 8.5% |
| 3: SMEM Caching | `2980.3` | 12.8% |
| 4: 1D Blocktiling | `8474.7` | 36.5% |
| 5: 2D Blocktiling | `15971.7` | 68.7% |
| 7: Avoid Bank Conflicts (Linearize) | `16213.4` | 69.7% |
| 8: Avoid Bank Conflicts (Offset) | `16459.2` | 70.8% |
| 11: Double Buffering | `17278.3` | 74.3% |
| 6: Vectorized Mem Access | `18237.3` | 78.4% |
| 9: Autotuning | `19721.0` | 84.8% |
| 10: Warptiling | `21779.3` | 93.7% |
| 0: cuBLAS | `23249.6` | 100.0% |
<!-- benchmark_results -->
## Setup
1. Install dependencies: CUDA toolkit 12, Python (+ Seaborn), CMake, Ninja. See [environment.yml](environment.yml).
1. Configure NVCC compilation parameters. Look up your GPUs compute
capability [here](https://developer.nvidia.com/cuda-gpus). Then configure the `CMakeLists.txt` and change:
```cmake
set(CUDA_COMPUTE_CAPABILITY 80)
```
1. Build: `mkdir build && cd build && cmake .. && cmake --build .`
1. Run one of the kernels: `DEVICE=<device_id> ./sgemm <kernel number>`
1. Profiling via [NVIDIA Nsight Compute](https://developer.nvidia.com/nsight-compute) (ncu): `make profile KERNEL=<kernel number>`
Credit goes to [wangzyon/NVIDIA_SGEMM_PRACTICE](https://github.com/wangzyon/NVIDIA_SGEMM_PRACTICE) for the benchmarking setup.

View File

@@ -0,0 +1,108 @@
#include <cstdio>
#include <cublas_v2.h>
#include <cuda_runtime.h>
/*
* A stand-alone script to invoke & benchmark standard cuBLAS SGEMM performance
*/
int main(int argc, char *argv[]) {
int m = 2;
int k = 3;
int n = 4;
int print = 1;
cudaError_t cudaStat; // cudaMalloc status
cublasStatus_t stat; // cuBLAS functions status
cublasHandle_t handle; // cuBLAS context
int i, j;
float *a, *b, *c;
// malloc for a,b,c...
a = (float *)malloc(m * k * sizeof(float));
b = (float *)malloc(k * n * sizeof(float));
c = (float *)malloc(m * n * sizeof(float));
int ind = 11;
for (j = 0; j < m * k; j++) {
a[j] = (float)ind++;
}
ind = 11;
for (j = 0; j < k * n; j++) {
b[j] = (float)ind++;
}
ind = 11;
for (j = 0; j < m * n; j++) {
c[j] = (float)ind++;
}
// DEVICE
float *d_a, *d_b, *d_c;
// cudaMalloc for d_a, d_b, d_c...
cudaMalloc((void **)&d_a, m * k * sizeof(float));
cudaMalloc((void **)&d_b, k * n * sizeof(float));
cudaMalloc((void **)&d_c, m * n * sizeof(float));
stat = cublasCreate(&handle); // initialize CUBLAS context
cudaMemcpy(d_a, a, m * k * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_b, b, k * n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_c, c, m * n * sizeof(float), cudaMemcpyHostToDevice);
float alpha = 1.0f;
float beta = 0.5f;
if (print == 1) {
printf("alpha = %4.0f, beta = %4.0f\n", alpha, beta);
printf("A = (mxk: %d x %d)\n", m, k);
for (i = 0; i < m; i++) {
for (j = 0; j < k; j++) {
printf("%4.1f ", a[i * m + j]);
}
printf("\n");
}
printf("B = (kxn: %d x %d)\n", k, n);
for (i = 0; i < k; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", b[i * n + j]);
}
printf("\n");
}
printf("C = (mxn: %d x %d)\n", m, n);
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", c[i * n + j]);
}
printf("\n");
}
}
stat = cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, d_b, n,
d_a, k, &beta, d_c, n);
cudaMemcpy(c, d_c, m * n * sizeof(float), cudaMemcpyDeviceToHost);
if (print == 1) {
printf("\nC after SGEMM = \n");
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
printf("%4.1f ", c[i * n + j]);
}
printf("\n");
}
}
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
cublasDestroy(handle); // destroy CUBLAS context
free(a);
free(b);
free(c);
return EXIT_SUCCESS;
}

View File

@@ -0,0 +1,33 @@
banks_naive = lambda r, c: (r * 32 + c) % 32
banks_one_extra = lambda r, c: (r * 33 + c) % 32
ITEMS_PER_WARP = 8
def printBankConflicts(bank_fun):
for c in range(1):
banks = []
for i in range(32):
row = (i * ITEMS_PER_WARP) // 16
col = (i * ITEMS_PER_WARP + c) % 16
banks.append((i, row, col, bank_fun(row, col)))
print("Step", c, "\n", "\n".join(["(" + ",".join(str(x) for x in i) + ")" for i in banks]))
d = {k: 0 for k in range(32)}
for i in banks:
d[i[-1]] += 1
count = 0
for key, val in d.items():
if val > 0:
count += 1
print(
f"Bank conflicts (Step {c}): {sorted(d.items(), key=lambda item: item[1], reverse=True)[0][1]}, banks accessed: {count}/32\n"
)
print("---NAIVE---")
printBankConflicts(banks_naive, 32)
print("\n---EXTRA COL---")
printBankConflicts(banks_one_extra, 33)

View File

@@ -0,0 +1,115 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
WM_VALUES=(32 64 128 256)
WN_VALUES=(32 64 128 256)
WNITER_VALUES=(1 2 4 8)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
NUM_THREADS_VALUES=(128 256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
OUTPUT="../benchmark_results/kernel_10_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="0"
WARPSIZE=32
TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for BK in "${BK_VALUES[@]}"; do
for BM in "${BM_VALUES[@]}"; do
for BN in "${BN_VALUES[@]}"; do
for WM in "${WM_VALUES[@]}"; do
for WN in "${WN_VALUES[@]}"; do
for WN_ITER in "${WNITER_VALUES[@]}"; do
for TM in "${TM_VALUES[@]}"; do
for TN in "${TN_VALUES[@]}"; do
for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
echo ""
CONFIG_NUM=$(( CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
NUM_WARPS=$(( NUM_THREADS / 32 ))
if ! (( BN % WN == 0 && BM % WM == 0 )); then
echo "Error: BN % WN must be 0 and BM % WM must be 0."
continue
fi
if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
continue
fi
if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
continue
fi
WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BK == 0 )); then
echo "Error: (NUM_THREADS * 4) % BK must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BN == 0 )); then
echo "Error: (NUM_THREADS * 4) % BN must be 0."
continue
fi
if ! (( BN % (16 * TN) == 0 )); then
echo "Error: BN must be a multiple of 16 * TN."
continue
fi
if ! (( BM % (16 * TM) == 0 )); then
echo "Error: BM must be a multiple of 16 * TM."
continue
fi
if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K10_NUM_THREADS = .*/const uint K10_NUM_THREADS = $NUM_THREADS;/" $RUNNER
sed -i "s/const uint K10_BN = .*/const uint K10_BN = $BN;/" $RUNNER
sed -i "s/const uint K10_BM = .*/const uint K10_BM = $BM;/" $RUNNER
sed -i "s/const uint K10_BK = .*/const uint K10_BK = $BK;/" $RUNNER
sed -i "s/const uint K10_WM = .*/const uint K10_WM = $WM;/" $RUNNER
sed -i "s/const uint K10_WN = .*/const uint K10_WN = $WN;/" $RUNNER
sed -i "s/const uint K10_WNITER = .*/const uint K10_WNITER = $WN_ITER;/" $RUNNER
sed -i "s/const uint K10_TM = .*/const uint K10_TM = $TM;/" $RUNNER
sed -i "s/const uint K10_TN = .*/const uint K10_TN = $TN;/" $RUNNER
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 4 seconds if it doesn't finish
timeout -v 8 ./sgemm 10 | tee -a $OUTPUT
done
done
done
done
done
done
done
done
done

View File

@@ -0,0 +1,115 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
WM_VALUES=(32 64 128 256)
WN_VALUES=(32 64 128 256)
WNITER_VALUES=(1 2 4 8)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
NUM_THREADS_VALUES=(128 256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
OUTPUT="../benchmark_results/kernel_11_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="0"
WARPSIZE=32
TOTAL_CONFIGS="$(( ${#BK_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} * ${#WM_VALUES[@]} * ${#WN_VALUES[@]} * ${#WNITER_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#NUM_THREADS_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for BK in "${BK_VALUES[@]}"; do
for BM in "${BM_VALUES[@]}"; do
for BN in "${BN_VALUES[@]}"; do
for WM in "${WM_VALUES[@]}"; do
for WN in "${WN_VALUES[@]}"; do
for WN_ITER in "${WNITER_VALUES[@]}"; do
for TM in "${TM_VALUES[@]}"; do
for TN in "${TN_VALUES[@]}"; do
for NUM_THREADS in "${NUM_THREADS_VALUES[@]}"; do
echo ""
CONFIG_NUM=$(( CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
NUM_WARPS=$(( NUM_THREADS / 32 ))
if ! (( BN % WN == 0 && BM % WM == 0 )); then
echo "Error: BN % WN must be 0 and BM % WM must be 0."
continue
fi
if ! (( (BN / WN) * (BM / WM) == NUM_WARPS )); then
echo "Error: (BN / WN) * (BM / WM) must be equal to NUM_WARPS."
continue
fi
if ! (( (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) == 0 )); then
echo "Error: (WM * WN) % (WARPSIZE * TM * TN * WN_ITER) must be 0."
continue
fi
WM_ITER=$(( (WM * WN) / (WARPSIZE * TM * TN * WN_ITER) ))
if ! (( WM % WM_ITER == 0 && WN % WN_ITER == 0 )); then
echo "Error: WM % WM_ITER must be 0 and WN % WN_ITER must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BK == 0 )); then
echo "Error: (NUM_THREADS * 4) % BK must be 0."
continue
fi
if ! (( (NUM_THREADS * 4) % BN == 0 )); then
echo "Error: (NUM_THREADS * 4) % BN must be 0."
continue
fi
if ! (( BN % (16 * TN) == 0 )); then
echo "Error: BN must be a multiple of 16 * TN."
continue
fi
if ! (( BM % (16 * TM) == 0 )); then
echo "Error: BM must be a multiple of 16 * TM."
continue
fi
if ! (( (BM * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BM * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
if ! (( (BN * BK) % (4 * NUM_THREADS) == 0 )); then
echo "Error: (BN * BK) % (4 * NUM_THREADS) must be 0."
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K11_NUM_THREADS = .*/const uint K11_NUM_THREADS = $NUM_THREADS;/" $RUNNER
sed -i "s/const uint K11_BN = .*/const uint K11_BN = $BN;/" $RUNNER
sed -i "s/const uint K11_BM = .*/const uint K11_BM = $BM;/" $RUNNER
sed -i "s/const uint K11_BK = .*/const uint K11_BK = $BK;/" $RUNNER
sed -i "s/const uint K11_WM = .*/const uint K11_WM = $WM;/" $RUNNER
sed -i "s/const uint K11_WN = .*/const uint K11_WN = $WN;/" $RUNNER
sed -i "s/const uint K11_WNITER = .*/const uint K11_WNITER = $WN_ITER;/" $RUNNER
sed -i "s/const uint K11_TM = .*/const uint K11_TM = $TM;/" $RUNNER
sed -i "s/const uint K11_TN = .*/const uint K11_TN = $TN;/" $RUNNER
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$BK BM=$BM BN=$BN WM=$WM WN=$WN WN_ITER=$WN_ITER TM=$TM TN=$TN NUM_THREADS=$NUM_THREADS" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 8 seconds if it doesn't finish
timeout -v 8 ./sgemm 11 | tee -a $OUTPUT
done
done
done
done
done
done
done
done
done

View File

@@ -0,0 +1,86 @@
#!/usr/bin/env bash
set -u
# Define the range of values for each parameter
BK_VALUES=(8 16 32 64)
TM_VALUES=(4 8 16 32)
TN_VALUES=(4 8 16 32)
BM_VALUES=(64 128 256)
BN_VALUES=(64 128 256)
NUM_THREADS_VALUES=(256)
cd "$(dirname "$0")"
cd "../build"
RUNNER="../src/runner.cu"
KERNEL="../src/kernels/9_kernel_autotuned.cuh"
OUTPUT="../benchmark_results/kernel_9_autotune_results.txt"
# Clear the output file
echo "" > $OUTPUT
# Set GPU to use
export DEVICE="2"
TOTAL_CONFIGS="$(( ${#NUM_THREADS_VALUES[@]} * ${#BK_VALUES[@]} * ${#TM_VALUES[@]} * ${#TN_VALUES[@]} * ${#BM_VALUES[@]} * ${#BN_VALUES[@]} ))"
CONFIG_NUM=0
# Loop through all combinations of parameters
for bk in ${BK_VALUES[@]}; do
for tm in ${TM_VALUES[@]}; do
for tn in ${TN_VALUES[@]}; do
for bm in ${BM_VALUES[@]}; do
for bn in ${BN_VALUES[@]}; do
for nt in ${NUM_THREADS_VALUES[@]}; do
echo ""
CONFIG_NUM=$(( $CONFIG_NUM + 1 ))
# skip configurations that don't fullfil preconditions
config="BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NT=$nt"
if [[ $(( ($nt * 4) % bk )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BK = $(( ($nt * 4) % bk )) != 0))"
continue
fi
if [[ $(( ($nt * 4) % bn )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (NUM_THREADS * 4) % BN = $(( ($nt * 4) % bn )) != 0))"
continue
fi
if [[ $(( $bn % (16 * $tn ) )) -ne 0 ]]; then
echo "QUANTIZATION: Skipping $config because BN % (16 * TN) = $(( $bn % (16 * $tn ) )) != 0))"
continue
fi
if [[ $(( $bm % (16 * $tm ) )) -ne 0 ]]; then
echo "QUANTIZATION: Skipping $config because BM % (16 * TM) = $(( $bm % (16 * $tm ) )) != 0))"
continue
fi
if [[ $(( ($bm * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (BM * BK) % (4 * NUM_THREADS) = $(( ($bm * $bk) % ( 4 * 256 ) )) != 0))"
continue
fi
if [[ $(( ($bn * $bk) % ( 4 * $nt ) )) -ne 0 ]]; then
echo "VECTORIZE: Skipping $config because (BN * BK) % (4 * NUM_THREADS) = $(( ($bn * $bk) % ( 4 * 256 ) )) != 0))"
continue
fi
# Update the parameters in the source code
sed -i "s/const uint K9_BK = .*/const uint K9_BK = $bk;/" $RUNNER
sed -i "s/const uint K9_TM = .*/const uint K9_TM = $tm;/" $RUNNER
sed -i "s/const uint K9_TN = .*/const uint K9_TN = $tn;/" $RUNNER
sed -i "s/const uint K9_BM = .*/const uint K9_BM = $bm;/" $RUNNER
sed -i "s/const uint K9_BN = .*/const uint K9_BN = $bn;/" $RUNNER
sed -i "s/const int K9_NUM_THREADS = .*/const int K9_NUM_THREADS = $nt;/" $KERNEL
# Rebuild the program
make
echo "($CONFIG_NUM/$TOTAL_CONFIGS): BK=$bk TM=$tm TN=$tn BM=$bm BN=$bn NUM_THREADS=$nt" |& tee -a $OUTPUT
# Run the benchmark and get the result
# Kill the program after 4 seconds if it doesn't finish
timeout -v 4 ./sgemm 9 | tee -a $OUTPUT
done
done
done
done
done
done

View File

@@ -0,0 +1,168 @@
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <runner.cuh>
#include <vector>
#define cudaCheck(err) (cudaCheck(err, __FILE__, __LINE__))
const std::string errLogFile = "matrixValidationFailure.txt";
int main(int argc, char **argv) {
if (argc != 2) {
std::cerr << "Please select a kernel (range 0 - 12, 0 for NVIDIA cuBLAS)"
<< std::endl;
exit(EXIT_FAILURE);
}
// get kernel number
int kernel_num = std::stoi(argv[1]);
if (kernel_num < 0 || kernel_num > 12) {
std::cerr << "Please enter a valid kernel number (0-12)" << std::endl;
exit(EXIT_FAILURE);
}
// get environment variable for device
int deviceIdx = 0;
if (getenv("DEVICE") != NULL) {
deviceIdx = atoi(getenv("DEVICE"));
}
cudaCheck(cudaSetDevice(deviceIdx));
printf("Running kernel %d on device %d.\n", kernel_num, deviceIdx);
// print some device info
// CudaDeviceInfo();
// Declare the handle, create the handle, cublasCreate will return a value of
// type cublasStatus_t to determine whether the handle was created
// successfully (the value is 0)
cublasHandle_t handle;
if (cublasCreate(&handle)) {
std::cerr << "Create cublas handle error." << std::endl;
exit(EXIT_FAILURE);
};
// Using cudaEvent for gpu stream timing, cudaEvent is equivalent to
// publishing event tasks in the target stream
float elapsed_time;
cudaEvent_t beg, end;
cudaEventCreate(&beg);
cudaEventCreate(&end);
// cuBLAS FLOPs ceiling is reached at 8192
std::vector<int> SIZE = {128, 256, 512, 1024, 2048, 4096};
long m, n, k, max_size;
max_size = SIZE[SIZE.size() - 1];
std::cout << "Max size: " << max_size << std::endl;
float alpha = 0.5, beta = 3.0; // GEMM input parameters, C=α*AB+β*C
float *A = nullptr, *B = nullptr, *C = nullptr,
*C_ref = nullptr; // host matrices
float *dA = nullptr, *dB = nullptr, *dC = nullptr,
*dC_ref = nullptr; // device matrices
A = (float *)malloc(sizeof(float) * max_size * max_size);
B = (float *)malloc(sizeof(float) * max_size * max_size);
C = (float *)malloc(sizeof(float) * max_size * max_size);
C_ref = (float *)malloc(sizeof(float) * max_size * max_size);
randomize_matrix(A, max_size * max_size);
randomize_matrix(B, max_size * max_size);
randomize_matrix(C, max_size * max_size);
cudaCheck(cudaMalloc((void **)&dA, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dB, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dC, sizeof(float) * max_size * max_size));
cudaCheck(cudaMalloc((void **)&dC_ref, sizeof(float) * max_size * max_size));
cudaCheck(cudaMemcpy(dA, A, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dB, B, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dC, C, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
cudaCheck(cudaMemcpy(dC_ref, C, sizeof(float) * max_size * max_size,
cudaMemcpyHostToDevice));
int repeat_times = 50;
for (int size : SIZE) {
m = n = k = size;
std::cout << "dimensions(m=n=k) " << m << ", alpha: " << alpha
<< ", beta: " << beta << std::endl;
// Verify the correctness of the calculation, and execute it once before the
// kernel function timing to avoid cold start errors
if (kernel_num != 0) {
run_kernel(0, m, n, k, alpha, dA, dB, beta, dC_ref,
handle); // cuBLAS
run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC,
handle); // Executes the kernel, modifies the result matrix
cudaCheck(cudaDeviceSynchronize());
cudaCheck(cudaGetLastError()); // Check for async errors during kernel run
cudaMemcpy(C, dC, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
cudaMemcpy(C_ref, dC_ref, sizeof(float) * m * n, cudaMemcpyDeviceToHost);
if (!verify_matrix(C_ref, C, m * n)) {
std::cout
<< "Failed to pass the correctness verification against NVIDIA "
"cuBLAS."
<< std::endl;
if (m <= 128) {
std::cout << " Logging faulty output into " << errLogFile << "\n";
std::ofstream fs;
fs.open(errLogFile);
fs << "A:\n";
print_matrix(A, m, n, fs);
fs << "B:\n";
print_matrix(B, m, n, fs);
fs << "C:\n";
print_matrix(C, m, n, fs);
fs << "Should:\n";
print_matrix(C_ref, m, n, fs);
}
exit(EXIT_FAILURE);
}
}
cudaEventRecord(beg);
for (int j = 0; j < repeat_times; j++) {
// We don't reset dC between runs to save time
run_kernel(kernel_num, m, n, k, alpha, dA, dB, beta, dC, handle);
}
cudaEventRecord(end);
cudaEventSynchronize(beg);
cudaEventSynchronize(end);
cudaEventElapsedTime(&elapsed_time, beg, end);
elapsed_time /= 1000.; // Convert to seconds
long flops = 2 * m * n * k;
printf(
"Average elapsed time: (%7.6f) s, performance: (%7.1f) GFLOPS. size: "
"(%ld).\n",
elapsed_time / repeat_times,
(repeat_times * flops * 1e-9) / elapsed_time, m);
fflush(stdout);
// make dC and dC_ref equal again (we modified dC while calling our kernel
// for benchmarking)
cudaCheck(cudaMemcpy(dC, dC_ref, sizeof(float) * m * n,
cudaMemcpyDeviceToDevice));
}
// Free up CPU and GPU space
free(A);
free(B);
free(C);
free(C_ref);
cudaFree(dA);
cudaFree(dB);
cudaFree(dC);
cudaFree(dC_ref);
cublasDestroy(handle);
return 0;
};

View File

@@ -0,0 +1,46 @@
#include <cuda_runtime.h>
#include <iostream>
#include <vector>
__global__ void kernel(uint *A, uint *B, int row) {
auto x = threadIdx.x / 4;
auto y = threadIdx.x % 4;
A[x * row + y] = x;
B[x * row + y] = y;
}
int main(int argc, char **argv) {
uint *Xs, *Ys;
uint *Xs_d, *Ys_d;
uint SIZE = 4;
Xs = (uint *)malloc(SIZE * SIZE * sizeof(uint));
Ys = (uint *)malloc(SIZE * SIZE * sizeof(uint));
cudaMalloc((void **)&Xs_d, SIZE * SIZE * sizeof(uint));
cudaMalloc((void **)&Ys_d, SIZE * SIZE * sizeof(uint));
dim3 grid_size(1, 1, 1);
dim3 block_size(4 * 4);
kernel<<<grid_size, block_size>>>(Xs_d, Ys_d, 4);
cudaMemcpy(Xs, Xs_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
cudaMemcpy(Ys, Ys_d, SIZE * SIZE * sizeof(uint), cudaMemcpyDeviceToHost);
cudaDeviceSynchronize();
for (int row = 0; row < SIZE; ++row) {
for (int col = 0; col < SIZE; ++col) {
std::cout << "[" << Xs[row * SIZE + col] << "|" << Ys[row * SIZE + col]
<< "] ";
}
std::cout << "\n";
}
cudaFree(Xs_d);
cudaFree(Ys_d);
free(Xs);
free(Ys);
}

View File

@@ -0,0 +1,14 @@
#pragma once
#include "kernels/10_kernel_warptiling.cuh"
#include "kernels/11_kernel_double_buffering.cuh"
#include "kernels/12_kernel_double_buffering.cuh"
#include "kernels/1_naive.cuh"
#include "kernels/2_kernel_global_mem_coalesce.cuh"
#include "kernels/3_kernel_shared_mem_blocking.cuh"
#include "kernels/4_kernel_1D_blocktiling.cuh"
#include "kernels/5_kernel_2D_blocktiling.cuh"
#include "kernels/6_kernel_vectorize.cuh"
#include "kernels/7_kernel_resolve_bank_conflicts.cuh"
#include "kernels/8_kernel_bank_extra_col.cuh"
#include "kernels/9_kernel_autotuned.cuh"

View File

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

View File

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

View File

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

View File

@@ -0,0 +1,29 @@
#pragma once
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
/*
Matrix sizes:
MxK * KxN = MxN
*/
__global__ void sgemm_naive(int M, int N, int K, float alpha, const float *A,
const float *B, float beta, float *C) {
const uint x = blockIdx.x * blockDim.x + threadIdx.x;
const uint y = blockIdx.y * blockDim.y + threadIdx.y;
// if statement is necessary to make things work under tile quantization
if (x < M && y < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) {
tmp += A[x * K + i] * B[i * N + y];
}
// C = α*(A@B)+β*C
C[x * N + y] = alpha * tmp + beta * C[x * N + y];
}
}

View File

@@ -0,0 +1,24 @@
#pragma once
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
template <const uint BLOCKSIZE>
__global__ void sgemm_global_mem_coalesce(int M, int N, int K, float alpha,
const float *A, const float *B,
float beta, float *C) {
const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE);
const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);
// if statement is necessary to make things work under tile quantization
if (cRow < M && cCol < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) {
tmp += A[cRow * K + i] * B[i * N + cCol];
}
C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol];
}
}

View File

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

View File

@@ -0,0 +1,80 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM>
__global__ void sgemm1DBlocktiling(int M, int N, int K, float alpha,
const float *A, const float *B, float beta,
float *C) {
// If we flip x and y here we get ~30% less performance for large matrices.
// The current, 30% faster configuration ensures that blocks with sequential
// blockIDs access columns of B sequentially, while sharing the same row of A.
// The slower configuration would share columns of A, but access into B would
// be non-sequential. So the faster configuration has better spatial locality
// and hence a greater L2 hit rate.
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// each warp will calculate 32*TM elements, with 32 being the columnar dim.
const int threadCol = threadIdx.x % BN;
const int threadRow = threadIdx.x / BN;
// allocate space for the current blocktile in SMEM
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// todo: adjust this to each thread to load multiple entries and
// better exploit the cache sizes
assert(BM * BK == blockDim.x);
assert(BN * BK == blockDim.x);
const uint innerColA = threadIdx.x % BK; // warp-level GMEM coalescing
const uint innerRowA = threadIdx.x / BK;
const uint innerColB = threadIdx.x % BN; // warp-level GMEM coalescing
const uint innerRowB = threadIdx.x / BN;
// allocate thread-local cache for results in registerfile
float threadResults[TM] = {0.0};
// outer loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
__syncthreads();
// advance blocktile
A += BK;
B += BK * N;
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// we make the dotproduct loop the outside loop, which facilitates
// reuse of the Bs entry, which we can cache in a tmp var.
float tmpB = Bs[dotIdx * BN + threadCol];
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
threadResults[resIdx] +=
As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
}
}
__syncthreads();
}
// write out the results
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
C[(threadRow * TM + resIdx) * N + threadCol] =
alpha * threadResults[resIdx] +
beta * C[(threadRow * TM + resIdx) * N + threadCol];
}
}

View File

@@ -0,0 +1,102 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void __launch_bounds__((BM * BN) / (TM * TN), 1)
sgemm2DBlocktiling(int M, int N, int K, float alpha, const float *A,
const float *B, float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
const uint totalResultsBlocktile = BM * BN;
// A thread is responsible for calculating TM*TN elements in the blocktile
const uint numThreadsBlocktile = totalResultsBlocktile / (TM * TN);
// ResultsPerBlock / ResultsPerThread == ThreadsPerBlock
assert(numThreadsBlocktile == blockDim.x);
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
const uint innerRowA = threadIdx.x / BK;
const uint innerColA = threadIdx.x % BK;
// calculates the number of rows of As that are being loaded in a single step
// by a single block
const uint strideA = numThreadsBlocktile / BK;
const uint innerRowB = threadIdx.x / BN;
const uint innerColB = threadIdx.x % BN;
// for both As and Bs we want each load to span the full column-width, for
// better GMEM coalescing (as opposed to spanning full row-width and iterating
// across columns)
const uint strideB = numThreadsBlocktile / BN;
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
// register caches for As and Bs
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
As[(innerRowA + loadOffset) * BK + innerColA] =
A[(innerRowA + loadOffset) * K + innerColA];
}
for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
Bs[(innerRowB + loadOffset) * BN + innerColB] =
B[(innerRowB + loadOffset) * N + innerColB];
}
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN] =
alpha * threadResults[resIdxM * TN + resIdxN] +
beta * C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN];
}
}
}

View File

@@ -0,0 +1,98 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmVectorize(int M, int N, int K, float alpha, float *A,
float *B, float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,103 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmResolveBankConflicts(int M, int N, int K, float alpha,
float *A, float *B, float beta,
float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
// "linearize" Bs while storing it
tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 0) * 16 + innerColB / 2] = tmp.x;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 1) * 16 + innerColB / 2] = tmp.y;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 2) * 16 + innerColB / 2] = tmp.z;
Bs[((innerColB % 2) * 4 + innerRowB * 8 + 3) * 16 + innerColB / 2] = tmp.w;
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[(dotIdx * 8 + i) * 16 + threadCol];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,103 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemmResolveBankExtraCol(int M, int N, int K, float alpha,
float *A, float *B, float beta,
float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// BN/TN are the number of threads to span a column
const int threadCol = threadIdx.x % (BN / TN);
const int threadRow = threadIdx.x / (BN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
const int extraCols = 5;
__shared__ float Bs[BK * (BN + extraCols)];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
// transpose A while loading it
float4 tmp =
reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;
tmp = reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 0] = tmp.x;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 1] = tmp.y;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 2] = tmp.z;
Bs[innerRowB * (BN + extraCols) + innerColB * 4 + 3] = tmp.w;
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * (BN + extraCols) + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[resIdxM * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
__syncthreads();
}
// write out the results
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0];
// perform GEMM update in reg
tmp.x = alpha * threadResults[resIdxM * TN + resIdxN] + beta * tmp.x;
tmp.y = alpha * threadResults[resIdxM * TN + resIdxN + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[resIdxM * TN + resIdxN + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[resIdxM * TN + resIdxN + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(
&C[(threadRow * TM + resIdxM) * N + threadCol * TN + resIdxN])[0] =
tmp;
}
}
}

View File

@@ -0,0 +1,127 @@
#pragma once
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
const int K9_NUM_THREADS = 256;
template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void __launch_bounds__(K9_NUM_THREADS)
sgemmAutotuned(int M, int N, int K, float alpha, float *A, float *B,
float beta, float *C) {
const uint cRow = blockIdx.y;
const uint cCol = blockIdx.x;
// size of warptile
constexpr int WM = TM * 16;
constexpr int WN = TN * 16;
// iterations of warptile
constexpr int WMITER = CEIL_DIV(BM, WM);
constexpr int WNITER = CEIL_DIV(BN, WN);
// Placement of the thread in the warptile
const int threadCol = threadIdx.x % (WN / TN);
const int threadRow = threadIdx.x / (WN / TN);
// allocate space for the current blocktile in smem
__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];
// Move blocktile to beginning of A's row and B's column
A += cRow * BM * K;
B += cCol * BN;
C += cRow * BM * N + cCol * BN;
// calculating the indices that this thread will load into SMEM
// we'll load 128bit / 32bit = 4 elements per thread at each step
const uint innerRowA = threadIdx.x / (BK / 4);
const uint innerColA = threadIdx.x % (BK / 4);
constexpr uint rowStrideA = (K9_NUM_THREADS * 4) / BK;
const uint innerRowB = threadIdx.x / (BN / 4);
const uint innerColB = threadIdx.x % (BN / 4);
constexpr uint rowStrideB = K9_NUM_THREADS / (BN / 4);
// allocate thread-local cache for results in registerfile
float threadResults[WMITER * WNITER * TM * TN] = {0.0};
float regM[TM] = {0.0};
float regN[TN] = {0.0};
// outer-most loop over block tiles
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
// populate the SMEM caches
for (uint offset = 0; offset + rowStrideA <= BM; offset += rowStrideA) {
float4 tmp = reinterpret_cast<float4 *>(
&A[(innerRowA + offset) * K + innerColA * 4])[0];
// transpose A while storing it
As[(innerColA * 4 + 0) * BM + innerRowA + offset] = tmp.x;
As[(innerColA * 4 + 1) * BM + innerRowA + offset] = tmp.y;
As[(innerColA * 4 + 2) * BM + innerRowA + offset] = tmp.z;
As[(innerColA * 4 + 3) * BM + innerRowA + offset] = tmp.w;
}
for (uint offset = 0; offset + rowStrideB <= BK; offset += rowStrideB) {
reinterpret_cast<float4 *>(
&Bs[(innerRowB + offset) * BN + innerColB * 4])[0] =
reinterpret_cast<float4 *>(
&B[(innerRowB + offset) * N + innerColB * 4])[0];
}
__syncthreads();
for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
// calculate per-thread results
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
// block into registers
for (uint i = 0; i < TM; ++i) {
regM[i] = As[dotIdx * BM + (wmIdx * WM) + threadRow * TM + i];
}
for (uint i = 0; i < TN; ++i) {
regN[i] = Bs[dotIdx * BN + (wnIdx * WN) + threadCol * TN + i];
}
for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
threadResults[(wmIdx * TM + resIdxM) * (WNITER * TN) +
wnIdx * TN + resIdxN] +=
regM[resIdxM] * regN[resIdxN];
}
}
}
}
}
__syncthreads();
// advance blocktile
A += BK; // move BK columns to right
B += BK * N; // move BK rows down
}
// write out the results
for (uint wmIdx = 0; wmIdx < WMITER; ++wmIdx) {
for (uint wnIdx = 0; wnIdx < WNITER; ++wnIdx) {
float *C_interim = C + (wmIdx * WM * N) + (wnIdx * WN);
for (uint resIdxM = 0; resIdxM < TM; resIdxM += 1) {
for (uint resIdxN = 0; resIdxN < TN; resIdxN += 4) {
// load C vector into registers
float4 tmp = reinterpret_cast<float4 *>(
&C_interim[(threadRow * TM + resIdxM) * N + threadCol * TN +
resIdxN])[0];
// perform GEMM update in reg
const int i =
(wmIdx * TM + resIdxM) * (WNITER * TN) + wnIdx * TN + resIdxN;
tmp.x = alpha * threadResults[i + 0] + beta * tmp.x;
tmp.y = alpha * threadResults[i + 1] + beta * tmp.y;
tmp.z = alpha * threadResults[i + 2] + beta * tmp.z;
tmp.w = alpha * threadResults[i + 3] + beta * tmp.w;
// write back
reinterpret_cast<float4 *>(&C_interim[(threadRow * TM + resIdxM) * N +
threadCol * TN + resIdxN])[0] =
tmp;
}
}
}
}
}

View File

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

View File

@@ -0,0 +1,26 @@
#pragma once
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
void cudaCheck(cudaError_t error, const char *file,
int line); // CUDA error check
void CudaDeviceInfo(); // print CUDA information
void range_init_matrix(float *mat, int N);
void randomize_matrix(float *mat, int N);
void zero_init_matrix(float *mat, int N);
void copy_matrix(const float *src, float *dest, int N);
void print_matrix(const float *A, int M, int N, std::ofstream &fs);
bool verify_matrix(float *mat1, float *mat2, int N);
float get_current_sec(); // Get the current moment
float cpu_elapsed_time(float &beg, float &end); // Calculate time difference
void run_kernel(int kernel_num, int m, int n, int k, float alpha, float *A,
float *B, float beta, float *C, cublasHandle_t handle);