# cuda

**Repository Path**: chasoul/cuda

## Basic Information

- **Project Name**: cuda
- **Description**:  使用cuda来实现矩阵乘法及相关性能分析
- **Primary Language**: C++
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 1
- **Created**: 2023-10-20
- **Last Updated**: 2023-10-20

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# cuda学习笔记

* 使用cuda进行矩阵乘法，假设A矩阵的大小为M*K，B矩阵的大小为K*N，结果矩阵C的大小为M*N
* 一个entry即C中的结果，即C中包含M*N个entry
* 实验时假设M=4096+7,K=4096+8,N=4096+9，Time Cost单位为ms，Occupancy括号后面的是一个SM可以驻存的block数量，Reg/T表示一个线程所需寄存器数量
* 术语

| Term | Iterpretation |
|:-----|:--------------|
| GMEM | global memory |
| SMEM | shared memory |
------------------------

## 1. P100设备参数
 | Key                                           | Value |
 |:----------------------------------------------|:------|
  CUDA Driver Version / Runtime Version          | 11.4 / 11.0
  CUDA Capability Major/Minor version number:    | 6.0
  Total amount of global memory:                 | 16281 MBytes (17071734784 bytes)
  (56) Multiprocessors, ( 64) CUDA Cores/MP:     | 3584 CUDA Cores
  GPU Max Clock rate:                            | 1329 MHz (1.33 GHz)
  Memory Clock rate:                             | 715 Mhz
  Memory Bus Width:                              | 4096-bit
  L2 Cache Size:                                 | 4194304 bytes
  Maximum Texture Dimension Size (x,y,z)         | 1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)
  Maximum Layered 1D Texture Size, (num) layers  | 1D=(32768), 2048 layers
  Maximum Layered 2D Texture Size, (num) layers  | 2D=(32768, 32768), 2048 layers
  Total amount of constant memory:               | 65536 bytes
  Total amount of shared memory per block:       | 49152 bytes
  Total number of registers available per block: | 65536
  Warp size:                                     | 32
  Maximum number of threads per multiprocessor:  | 2048
  Maximum number of threads per block:           | 1024
  Max dimension size of a thread block (x,y,z):  | (1024, 1024, 64)
  Max dimension size of a grid size    (x,y,z):  | (2147483647, 65535, 65535)
  Maximum memory pitch:                          | 2147483647 bytes
  Texture alignment:                             | 512 bytes
  Concurrent copy and kernel execution:          | Yes with 2 copy engine(s)
  Run time limit on kernels:                     | No
  Integrated GPU sharing Host Memory:            | No
  Support host page-locked memory mapping:       | Yes
  Alignment requirement for Surfaces:            | Yes
  Device has ECC support:                        | Enabled
  Device supports Unified Addressing (UVA):      | Yes
  Device supports Compute Preemption:            | Yes
  Supports Cooperative Kernel Launch:            | Yes
  Supports MultiDevice Co-op Kernel Launch:      | Yes
  Device PCI Domain ID / Bus ID / location ID:   | 0 / 4 / 0
  -----------------------------------------------------------

## 2. 笔记
* CUDA的每个SM的寄存器相关：寄存器分配的粒度是256，由于最小调度单位是warp，因此每个线程寄存器分配的粒度是256/32=8个。每个寄存器是4字节。
* 内存模型：
    * [Local Memory](https://www.zhihu.com/tardis/bd/art/344739356?source_id=1001)不是寄存器，而是从Global Memory上为每一个线程分配的一块内存。其优势在于运行时自动分配；使用结束后自动释放；拥有更精细的缓存策略；强制合并访问等。需要注意的是，寄存器资源使用过多时，编译器可能将其溢出到Local Memory，从而影响性能。可以使用`-maxrregcount=32`来限制kernel使用寄存器的数量。
    * Constant Memory是用来存储常量数据的，其对应的数据需要从host端进行全局声明，每个GPU只可以声明不超过64KB的Constant Memory。Constant Memory必须使用`cudaMemcpyToSymbol`来进行初始化，如：
        ```C
        __constant__ float const_data[256];
        float data[256];
        cudaMemcpyToSymbol(const_data, data, sizeof(data));
        cudaMemcpyFromSymbol(data, const_data, sizeof(data));
        ```
    * [Texture Memory](https://zhuanlan.zhihu.com/p/569971933)和Constant Memory一样是一种只读内存。


## 3. SGEMM迭代过程

### 3.1 sgemm_shared_mem_block

* 配置和编译后

|      Variable    |       Value      |                     Note                    |
|:-----------------|:-----------------|:--------------------------------------------|
| BLOCK_SIZE       | 32               | the width and length of each block          |
| grid (x,y,z)     | (M/32, N/32, 1)  | grid size                                   |
| block (x,y,z)    | (32*32, 1, 1)    | block size                                  |
| register num     | 28               | the number of registers used                |
| shared memory    | 8192B            | the size of shared memory used per thread   |
| constant memory  | 368B             | the size of constant memory used per thread |
-------------------------------------------------------------------------------------

这里要求每一个grid中每一个block的线程数量是x的平方，如16，64，256等。每一个block中的每一个线程都将从Global Memory从读取两个数据，分别是
`A[blockRow * BLOCK_SIZE + threadRow, blockCol * BLOCK_SIZE + threadCol]`，这里的`BLOCK_SIZE`的平方就等于一个block中线程的数量，然后每个线程会将读取到的数据写入到`As[threadRow, threadCol]`。
* 由于CUDA在GMEM上的访存粒度是32B，因此只要BLOCK_SIZE是4的倍数，那么对GMEM的访存就是合并的，并且没有浪费。
* 由于线程号连续的32个线程会组成一个warp，而对As的访问也是连续的，且SMEM包含32个bank，因此不存在bank写冲突。

从线程的角度看，这里调整BLOCK_SIZE并不会影响寄存器使用的数量，由于一个SM最多包含65536个寄存器，因此一个SM最多可以容纳约min(2340,2048)个线程，即64个warp。由于一个SM最大可以容纳64个warp，因此occupancy可以达到100%。
从block的角度看，由于一个SM最多容纳2048个线程，因此一个SM可以容纳的BLOCK数量为 $\frac{2048}{\text{BLOCK\_SIZE}^2}$ 。一个BLOCK所需的SMEM为 $\text{BLOCK\_SIZE}^2 \times 2 \times sizeof(float)$ 。综上，一个SM的SMEM占用的计算公式可以描述如下：

$$SMEM=\frac{2048}{\text{BLOCK\_SIZE} ^2} \times \text{BLOCK\_SIZE}^2 \times 2 \times 4=16KB$$

在执行以下代码时：
```C
for (int idx = 0; idx < BLOCK_SIZE; ++idx) {
    sum += As[threadRow * BLOCK_SIZE + idx] * Bs[idx * BLOCK_SIZE + threadCol];
}
```
对于一个warp，每执行一条指令，可以看出对As和Bs的访问都不会存在bank读冲突。因为这里对As和Bs的访问不存对同一个bank不同数据的访问。

对于一个block，其计算访存比为：

$$\frac{K \times \text{BLOCK\_SIZE}^2}{(K \times \text{BLOCK\_SIZE} \times 2 + \text{BLOCK\_SIZE} ^ 2) \times sizeof(float)} = \frac{K \times \text{BLOCK\_SIZE}}{8 \times K + 4 \times \text{BLOCK\_SIZE}}$$

因此，BLOCK_SIZE越大，则计算访存比越大，考虑到一个block最大可以支持1024个线程，即BLOCK_SIZE最大值为32，计算访存比最大极限为4。

* 实验结果

|       Method        | Time Cost | Occupancy | FLOPs/Byte | Reg/T | SMEM |
|:--------------------|:----------|:----------|:-----------|:------|:-----|
| cublas sgemm        | 19.83     | --        | --         | --    | --   |
| sgemm_shared_mem_8  | 171.50    | 78% (25)  | ~1         | 40    | 512  |
| sgemm_shared_mem_16 | 66.03     | 100% (8)  | ~2         | 32    | 2048 |
| sgemm_shared_mem_32 | 62.71     | 100% (2)  | ~4         | 32    | 8192 |
---------------------------------------------------------------------------

- 当BLOCK_SIZE为16时，其性能比BLOCK_SIZE要高一倍多，但当BLOCK_SIZE增大到32时，其性能比BLOCK_SIZE为16要略好。猜测是因为BLOCK_SIZE为32时，一个SM只能驻留2个block，而当BLOCK_SIZE为16时，一个SM可以驻留8个block，从而导致同步(`__syncthreads()`)的耗时显著增大。
- 此方案大概达到cublas sgemm的30%。

### 3.2 sgemm_BM_BN_BK_TM

本方案将引入两个新的变量，分别是BK和TM，其中BK表示As的列数和Bs的行数，TM表示每个线程计算的结果的数量。

对于一个block，每个线程各需要从GMEM中读取 $\frac{K}{BK}$ 个A的值和B的值。每一次循环中，每个线程读取一个数据，从而构成BM * BK的As和BK * BN 的Bs。因此BM等于BN，BM * BK等于一个block的线程数。从而，有以下关于TM的公式成立：
$$TM = \frac{BM \times BN}{BM \times BK} = \frac{BM \times BN}{BK \times BN}=\frac{BM}{BK}=\frac{BN}{BK}$$

此时，每个block中的线程数量$NT$可以表示为：
$$NT=\frac{BM \times BN}{TM}=BM \times BK$$

对一个block，本方案的计算访存比为：

$$\frac{K \times BM \times BN}{(K \times (BM + BN) + BM \times BN) \times sizeof(float)}$$

可以看出，BM和BN越大，则计算访存比越大。由于BM等于BN，且一个block最大可以包含1024个线程，因此BM和BN的最大值$\frac{1024}{BK}$。

由于GMEM上的访存粒度是32B，即8个float，因此将BK最小设置为8可以达到最优，即BM和BN的最大值可以为128（$BM \times BK \le 1024$），从而得出计算访存比的最大极限值为64（暂不考虑寄存器和SMEM）。

从硬件上考虑，若将TM个结果都存储到SMEM，此时每个线程需要28个寄存器，则需要(BM*BN+BM*BK+BK*BN)*sizeof(float)字节的SMEM，由于BK等于8，因此BM和BN最大只能为64（每个SM的SMEM为48KB，BM和BN都是2的整数次幂的情况下）。

~~若将TM个结果存储在寄存器中，则每个线程需要56个寄存器，一个block需要占用$(BM \times BK + BN \times BK) \times sizeof(float)$的SMEM，因此一个SM可以最多可以支持约1170个线程。将上述BM=BN=64，BK=8带入，可得一个SM可以容纳2个block，即32个warp，得出occupancy为50%，而计算访存比的最大极限值为32。~~

* 实验结果

下表中`sgemm_BM_BN_BK_TM_`后跟的两个数字分别是BM和BK的取值。

|         Method          | Time Cost | Occupancy | FLOPs/Byte | Reg/T | SMEM | 
|:------------------------|:----------|:----------|:-----------|:------|:-----|
| cublas sgemm            | 19.83     | --        | --         | --    | --   |
| sgemm_BM_BN_BK_TM_16_8  | 77.10     | 75% (12)  | ~2         | 40    | 1024 |
| sgemm_BM_BN_BK_TM_32_8  | 52.63     | 50% (4)   | ~4         | 51    | 2048 |
| sgemm_BM_BN_BK_TM_64_8  | 37.21     | 50% (2)   | ~8         | 64    | 4096 |
| sgemm_BM_BN_BK_TM_64_16 | 46.50     | 50% (1)   | ~8         | 48    | 8192 |
-------------------------------------------------------------------------------

- 当BM等于64且BK等于8时效果最好，此时效果达到了cublas的50%。
- 当计算访存比为4时，此时BM等于32且BK等于8，此时耗时为52.63ms，相比3.1中BLOCK_SIZE等于32时，尽管Occupancy从100%降低到了50%，但性能提升了大约16%。差异在于3.1中一个SM只能放下2个block，而本方案中一个SM可以放下4个block，这侧面证明了3.1中我们的猜想。
- 当BM等于64时，将BK从8增大到16，此时性能反而有所降低，这也可以证明在block中如果存在线程同步，此时估计性能要同时考虑一个SM可以驻存的block数量和计算访存比。

### 3.3 sgemm_BM_BN_BK_TM_TN

在3.1和3.2中，我们发现增大计算访存比和一个SM可以驻存的block的数量是可以有效提高性能的。在3.1和3.2中，每个线程每次外部循环只读取A和B矩阵中的一个值，在不浪费带宽的情况下这会显著影响计算访存比。

因此，在本实验中，我们引入了TN，$TM \times TN$等于每个线程计算的结果的数量。每次外部循环时，分别从A和B中取出BM\*BK和BK\*BN个值，然后计算出BM*BN个结果，因此其计算访存比的公式是与3.2中时相同的，但不同的是，此时需要的线程数量$NT$是：
$$NT = \frac{BM \times BN}{TM \times TN}$$

可以看出，每个SM中可以驻存的block数量在不考虑寄存器和SMEM消耗的情况下是3.2中的TN倍。

同时，在每次外部循环中，每个线程需要从A和B中读取的数据量$NI$为：
$$NI=\frac{BM \times BK}{NT} = \frac{BK \times TM \times TN}{BM} = \frac{BK \times TM \times TN}{BN}$$

* 实验结果

下表中`sgemm_`后跟的四个数字分别是BM、BK、TM和TN的取值。`I/T`表示每个线程在外部循环中每次需要从A或B读取的数据量。

|     Method      | Time Cost | Occupancy | FLOPs/Byte | Reg/T | SMEM | NT  | NI |
|:----------------|:----------|:----------|:-----------|:------|:-----|:----|:---|
| cublas sgemm    | 18.86     | --        | --         | --    | --   | --  | -- |
| sgemm_64_8_4_4  | 21.65     | 50% (4)   | ~8         | 56    | 4096 | 256 | 2  |
| sgemm_64_8_8_4  | 19.01     | 31% (5)   | ~8         | 96    | 4096 | 128 | 4  |
| sgemm_64_8_8_8  | 21.76     | 25% (8)   | ~8         | 128   | 4096 | 64  | 8  |
| sgemm_128_8_8_8 | 21.26     | 25% (2)   | ~16        | 128   | 8192 | 256 | 4  |
----------------------------------------------------------------------------------

* 由于每个线程需要计算TM\*TN个结果，并且需要计算多轮（$\frac{K}{BK}$），因此只有将这些结果先暂存在SMEM或者寄存器中。若存在SMEM中，以BM为64为例，此时一个block需要16KB，因此一个SM只能容纳3个block，是不合理的，因此我们将这些结果都暂存在寄存器中。从sgemm_64_8_4_4和sgemm_64_8_8_8可以看出，增大TM和TN将使每个线程所需寄存器数量显著增大，从而使得occupancy降低，进而降低性能。
* 本实验中最佳性能已接近cublas的99%。

### 3.4 sgemm_BM_BN_BK_TM_TN_vertorize

## 4. TODO list

- [ ] 详细了解Texture Memory
- [ ] 可以对任意shape的A矩阵和B矩阵进行矩阵乘法
- [ ] 调研warp tile