1、CUDA矩阵乘
1.1 背景知识
1.1.1 相对矩阵，当前的二维索引值
**第一步，线程和块索引映射到矩阵坐标**
//计算当前线程在矩阵的行列索引值
int ix = blockIdx.x * blockDim.x + threadIdx.x;
int iy = blockIdx.y * blockDim.y + threadIdx.y;
1.1.2数组存储顺序（或叫二维矩阵用一维表示）
第二步，矩阵坐标映射到全局内存中的索引/存储单元上
根据上面的公式，
第ix行第i列的元素，行优先是
idx = ix×n+i
第i行第iy列的元素，行优先是
idx = i×n+iy
两个矩阵分别为（ix，i）、（i，iy）。因此（ix，i）*（i，iy）=（ix，iy），也就是输出矩阵中（ix，iy）索引的元素。i控制得当，控制矩阵A的列，矩阵B的行。
换句话说，要得到输出矩阵（ix，iy）索引的值，那么就要有个中间变量（ix，x），（x，iy）。x来当读取值，也就是一个线程来读取A的一行，B的一列
1.2 全局内存矩阵乘
for(int i=0; i<k; i++)
{
sum += A[k * ix + i] * B[ i * n + iy];
}
利用循环i来读取A的第ix行、B的第iy列中的全部元素，读取方式如下：
矩阵A中一行有k个元素，第ix行第i列坐标是A[k * ix + i]，也即是k、ix不变，i在变；
矩阵B中一行有n个元素，第i行第iy列的坐标是B[ i * n + iy]，也即是n、iy不变，i在变。
上下i一致，即可以对应元素。
__global__void matrixMul(int *A,int *B, int *C, int m, int n, int k)
{
//计算当前线程在矩阵的行列索引值
int ix = blockIdx.x * blockDim.x + threadIdx.x;
int iy = blockIdx.y * blockDim.y + threadIdx.y;
//初试化一个临时变量
int sum = 0;
// 横坐标索引＜m，纵坐标索引＜n
if(ix<m && iy<n)
{
for(int i=0; i<k; i++)
{ // 相对全局内存的一维索引
sum += A[k * ix + i] * B[ i * n + iy];
}
C[ix * n + iy] = sum;
}
}
当程序走到这个循环时，ix，iy就确定了。也就是相当于A的第ix行，B的iy列确定了。进一步通过i来控制读取A每行、B每列的数据。（一行一列相乘）
for(int i=0; i<k; i++)
{ // 相对全局内存的一维索引
sum += A[k * ix + i] * B[ i * n + iy];
}
一个循环结束，矩阵A的第ix行与矩阵B的第iy列对应元素相乘后的累加和，矩阵C的第ix行第iy列的元素值。
C[ix * n + iy] = sum;
完整代码
# include<stdio.h>
# define m 4
# define k 4
# define n 4
# define BLOCK_SIZE 2
// global void中间有个空格
__global__ void matrixMul(int *A,int *B, int *C)
{
//计算当前线程在矩阵的行列索引值
int ix = blockIdx.x * blockDim.x + threadIdx.x;
int iy = blockIdx.y * blockDim.y + threadIdx.y;
//初试化一个临时变量
int sum = 0;
// 确保在矩阵C的范围内
if(ix<m && iy<n)
{
for(int i=0; i<k; i++)
{
sum += A[k * ix + i] * B[ i * n + iy];
}
C[ix * n + iy] = sum;
}
}
int main()
{
//分配主机内存
// int *h_a, int *h_b不行
int *h_a, *h_b, *h_c;
h_a = (int *)malloc(m*k*sizeof(int));
h_b = (int *)malloc(k*n*sizeof(int));
h_c = (int *)malloc(m*n*sizeof(int));
//初始化矩阵a，b
for(int i=0;i<m;i++)
{
for(int j=0;j<k;j++)
{
h_a[i*k+j]=1;
}
}
for(int i=0;i<k;i++)
{
for(int j=0;j<n;j++)
{
h_b[i*n+j]=2;
}
}
//分配设备内存
int *d_a, *d_b, *d_c;
cudaMalloc((void **)&d_a,sizeof(int)*m*k);
cudaMalloc((void **)&d_b,sizeof(int)*k*n);
cudaMalloc((void **)&d_c,sizeof(int)*m*n);
//将主机的值传给设备
cudaMemcpy(d_a,h_a,sizeof(int)*m*k,cudaMemcpyHostToDevice);
cudaMemcpy(d_b,h_b,sizeof(int)*k*n,cudaMemcpyHostToDevice);
//定义执行配置：要求线程数等于或大于矩阵元素的个数
//BLOCK_SIZE每块线程块的线程数
//举个简单的例子，假设m=5，BLOCK_SIZE=4。如果直接计算5/4，结果是1，
//这意味着只有一个线程块来处理行，这显然不够。但是，使用(5+4-1)/4计算，
//结果是2，这确保了有两个线程块来处理所有的行。
unsigned int grid_rows = (m+BLOCK_SIZE-1)/BLOCK_SIZE;
unsigned int grid_cols = (n+BLOCK_SIZE-1)/BLOCK_SIZE;
// dim3 dimGrid(grid_rows,grid_cols); 错
//cuda是列优先的
dim3 dimGrid(grid_cols,grid_rows);
dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
//执行核函数
matrixMul<<<dimGrid,dimBlock>>>(d_a,d_b,d_c);
//结果返回给主机
cudaMemcpy(h_c,d_c,m*n*sizeof(int),cudaMemcpyDeviceToHost);
//打印
for(int i=0;i<m;i++)
{
for(int j=0;j<n;j++)
{
printf("%d ", h_c[i*n+j]);
}
printf("\n");
}
//释放主机和设备
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
free(h_a);
free(h_b);
free(h_c);
return 0;
}
1.3 共享内存矩阵乘
__global__ void gpu_matrix_mult_shared(int *d_a, int *d_b, int *d_result, int M, int N, int K)
{
/*
__shared__: 这是一个CUDA关键字，
表示声明的变量是共享内存中的变量。共享内存是一块对于线程块中的所有线程可见的内存。
int tile_a[BLOCK_SIZE][BLOCK_SIZE];这是一个二维整数数组，
表示在共享内存中分配的二维数组。tile_a 被用于存储一个矩阵的子矩阵。
怎么理解被用于存储一个矩阵的子矩阵？
矩阵乘法通常涉及到将大矩阵分解为子矩阵，
然后每个线程块负责计算一个或多个子矩阵的乘法。
在矩阵乘法的上下文中，tile_a 存储了一个矩阵的子矩阵。当每个线程块负责计算某个子矩阵的乘法时，
将这个子矩阵加载到 tile_a 中，以便线程块中的每个线程可以快速访问和共享这部分数据。
*/
// 意思就是一个block处理一个tile
__shared__ int tile_a[BLOCK_SIZE][BLOCK_SIZE];
__shared__ int tile_b[BLOCK_SIZE][BLOCK_SIZE];
// 全局内存的位置
//前面有多少行（也就是正常坐标的列了），与下面呼应！！！
int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;
//前面有多少列
int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;
int tmp = 0;
int idx;
//读取每个矩阵小块
// for (int sub = 0; sub < gridDim.x; ++sub)
// 将全局矩阵中的小块加载到共享内存tile_a
// tile_a[threadIdx.y][threadIdx.x]=d_a[idx];
for (int sub = 0; sub <= N/BLOCK_SIZE; ++sub)
{
// idx = row * n + sub * BLOCK_SIZE + threadIdx.x;
// tile_a[threadIdx.y][threadIdx.x] = row<n && (sub * BLOCK_SIZE + threadIdx.x)<n? d_a[idx]:0;
// idx = (sub * BLOCK_SIZE + threadIdx.y) * n + col;
// tile_b[threadIdx.y][threadIdx.x] = col<n && (sub * BLOCK_SIZE + threadIdx.y)<n? d_b[idx]:0;
//详细版
//横着算的
int r = row;
int c = sub * BLOCK_SIZE + threadIdx.x;
idx = r * N + c;//看背景知识的图1，你就恍然大悟了
if (r >= M || c >= N)
{
// 将矩阵a的元素赋给tile_a
tile_a[threadIdx.y][threadIdx.x]=0;
}
else
{
// 将矩阵a的元素赋给tile_a
tile_a[threadIdx.y][threadIdx.x]=d_a[idx];
}
//竖着算的
//
r = sub * BLOCK_SIZE + threadIdx.y;
c = col;
// 一列有k个元素
idx = r * K + c;
if (c >= K || c >= N)
{
tile_b[threadIdx.y][threadIdx.x]=0;
}
else
{
tile_b[threadIdx.y][threadIdx.x]=d_b[idx];
}
__syncthreads();//使用 __syncthreads() 函数确保所有线程完成当前子矩阵的加载后，再进行下一步操作。
// 到了共享内存，那么就开始运算呀
//for中的k遍历tile_a的列，tile_b中的行
//相乘累加到一个temp中
for (int k = 0; k < BLOCK_SIZE; ++k)
{
tmp += tile_a[threadIdx.y][k] * tile_b[k][threadIdx.x];
}
__syncthreads();
}
if(row < M && col < K)
{
d_result[row * K + col] = tmp;
}
}
x[threadIdx.y][threadIdx.x]
对照上图，我们把一个int类型（四字节）的1024个元素的数组放到共享内存A中，每个int的索引对应到蓝框中，假设我们的块大小是(32,32)那么我们第一个线程束就是 threadIdx.y=0,threadIdx.x=0……31。
在二维数组（或矩阵）中，元素按行优先的方式存储。也就是说，内存中连续的元素属于同一行。要访问二维数组中的某个元素，首先确定其所在的行，然后再确定其在该行中的列位置。这种访问方式与我们习惯的矩阵表示方式一致，即行和列是垂直和水平方向。
2、其它CUDA面试题目
2.0 参考链接
2.1 warp sum
template<const int kWarpSize = WARP_SIZE>
__device__ __forceinline__ float warp_reduce_sum(float val) {
#pragma unroll
for (int mask = kWarpSize >> 1; mask >= 1; mask >>= 1) //p79，折半归约，位操作比整数操作更高效，相当于/2
{
val += __shfl_xor_sync(0xffffffff, val, mask);//通过线程数洗牌函数 __shfl_xor_sync异或交换数据p111
}
return val;
}
1）函数模板声明
template<const int kWarpSize = WARP_SIZE>
假设有一个CUDA线程束，包含32个线程（kWarpSize = 32），每个线程都有一个局部浮点数val。
调用
val = warp_reduce_sum<WARP_SIZE>(val);
val = warp_reduce_sum<NUM_WARPS>(val);
2）设备函数
__device__ __forceinline__
__device__表示该函数是在设备端（即GPU）执行的，而非主机端（CPU）。__forceinline__则是一个编译器提示，请求编译器尽可能内联此函数，以减少函数调用开销，提高执行效率。
减少函数调用开销：在C++中，函数调用会产生一定的开销，包括保存和恢复调用者的寄存器状态、跳转到函数地址等。对于非常简单的小函数，这些开销可能会超过函数体本身的执行时间。将这样的函数声明为内联函数，编译器会在调用处展开函数体，从而省去了函数调用的开销。
提高执行速度：通过函数体展开，CPU可以直接执行函数内的语句，避免了函数调用时的栈帧切换，理论上可以提高执行速度。
3）函数签名
float warp_reduce_sum(float val)
函数接受一个浮点数val作为输入，并返回一个浮点数，即线程束中所有输入val对应值的累加和。
4）循环展开
for (int mask = kWarpSize >> 1; mask >= 1; mask >>= 1)
循环变量mask从kWarpSize除以2开始（右移一位>> 1），每次迭代都右移一位，直到mask等于1。
循环从kWarpSize >> 1开始，每次迭代将mask右移一位，直到mask变为0。这种循环方式确保了warp内的每个线程都会参与归约操作，直到只剩下一个线程持有最终的归约结果。
5）__shfl_xor_sync操作
针对线程束内的线程
val += __shfl_xor_sync(0xffffffff, val, mask);
0xffffffff: 这是__shfl_xor_sync的第一个参数，它表示一个同步掩码。在这个例子中，它被设置为0xffffffff，这通常意味着所有的线程都会参与同步。
val: 这是当前线程的一个变量。
mask: 这是一个整数，用于确定与哪个线程进行值交换。
这行代码的执行过程是：
当前线程的val值与根据XOR shuffle计算得到的另一个线程的val值进行交换。
交换后的值被加到当前线程的val上。
__shfl_xor_sync(0xffffffff, val, mask)执行一个XOR shuffle操作，它将当前线程的val值与warp中另一个线程的val值相加。这个“另一个线程”是通过当前线程索引与mask进行XOR运算得到的
warp shuffle操作仅在同一warp内的线程之间有效，因此这种归约操作只能在一个warp的范围内进行。如果需要跨多个warp进行归约，则需要使用其他技术，如使用原子操作或多次归约步骤。
当然，以下是一个使用表格形式的简单示例，演示了`warp_reduce_sum`函数在一个大小为4的线程束（Warp）内如何并行累加求和。假设每个线程持有以下初始`val`值：
初始值：
| 线程索引 | 初始 `val` 值 |
|----------|--------------|
| 0 | 1.0 |
| 1 | 2.0 |
| 2 | 3.0 |
| 3 | 4.0 |
接下来，我们将按照`warp_reduce_sum`函数的逻辑，逐步展示线程束内数据交换和累加的过程。
**第一次循环（`mask = 2`）：**
| 线程索引 | 初始 `val` 值 | 目标线程索引（`tid ^ mask`） | 获取的 `val` 值 | 更新后的 `val` 值 |
|----------|--------------|------------------------------|----------------|-----------------|
| 0 | 1.0 | 2 (0^2） | 3.0 | 4.0 |
| 1 | 2.0 | 3 （1^2） | 4.0 | 6.0 |
| 2 | 3.0 | 0 | 1.0 | 4.0 |
| 3 | 4.0 | 1 | 2.0 | 6.0 |
**第二次循环（`mask = 1`）：**
| 线程索引 | 更新后的 `val` 值 | 目标线程索引（`tid ^ mask`） | 获取的 `val` 值 | 最终 `val` 值 |
|----------|-----------------|------------------------------|----------------|--------------|
| 0 | 4.0 | 1 | 6.0 | 10.0 |
| 1 | 6.0 | 0 | 4.0 | 10.0 |
| 2 | 4.0 | 3 | 6.0 | 10.0 |
| 3 | 6.0 | 2 | 4.0 | 10.0 |
经过两次循环后，每个线程的`val`值都变成了线程束内所有`val`值的总和，即10.0。这就是`warp_reduce_sum`函数如何在一个CUDA线程束内并行地对每个线程的`val`值进行累加求和的过程。尽管这个示例使用了一个较小的线程束（4个线程），但原理同样适用于更大规模的线程束（如典型的32个线程）。
tid ^ mask，目标线程索引是通过将当前线程索引（tid）与mask值进行按位异或（XOR）运算得到的。
所以后面if(lane0)，取值出来，其实lane1也可以，因为所有tid都是一样的结果了，10.0
5） #pragma unroll复制循环体代码
`#pragma unroll` 是CUDA编程中的一个编译指令，主要用于指导编译器对循环进行展开（unroll）。在CUDA C/C++编程中，特别是在GPU并行计算场景下，循环展开是一种常见的优化手段，可以帮助提升设备端（GPU）的执行效率。
当你在CUDA内核函数中使用 `#pragma unroll` 时，编译器会对紧跟其后的循环进行完全展开，即将循环体的代码重复多份，每份对应一次循环迭代。这样做可以减少循环控制逻辑带来的开销，并可能使GPU能够更好地并行执行循环体内的任务。
例如：
```cpp
__global__ void myKernel(int n) {
int idx = threadIdx.x;
#pragma unroll
for (int i = 0; i < n; ++i) {
// 这个循环会被展开
sharedMemory[idx] += someComputation(idx, i);
}
}
在此例中，如果 n 是一个小的常数，编译器将会生成对应的 n 份循环体代码。然而，需要注意的是过度的循环展开可能会增加GPU的注册文件压力以及编译后代码的大小，因此在使用时应结合实际情况权衡利弊。
2.2 warp max
// Warp Reduce Max
template<const int kWarpSize = WARP_SIZE>
__device__ __forceinline__ float warp_reduce_max(float val) {
#pragma unroll
for (int mask = kWarpSize >> 1; mask >= 1; mask >>= 1) {
val = fmaxf(val, __shfl_xor_sync(0xffffffff, val, mask));
}
return val;
}
__shfl_xor_sync(0xffffffff, val, mask)执行一个XOR shuffle操作，将当前线程的val值与warp中通过XOR运算得到的另一个线程的val值进行比较。
fmaxf(val, …)函数用于比较当前线程的val值与shuffle操作得到的值，并返回两者中的较大值。这个较大值随后被赋回给val。
__shfl_xor_sync 函数在 CUDA 中用于在 GPU 的 warp 内的线程之间执行 XOR shuffle 操作。这个函数会将当前线程的某个值（在本例中是 val）与通过 XOR 运算得到的另一个线程的值进行交换。
下面是一个简单的例子来解释 __shfl_xor_sync 的工作原理：
假设我们有一个包含 32 个线程的 warp，线程的索引从 0 到 31。我们想要使用 __shfl_xor_sync 来交换每个线程的值。
考虑 mask 的不同值：
当 mask = 16（即 kWarpSize >> 1 对于一个 32 线程的 warp）：
线程 0 会与线程 16 交换值。
线程 1 会与线程 17 交换值。
线程 2 会与线程 18 交换值。
以此类推，直到线程 15 与线程 31 交换值。
注意，由于 XOR 运算的性质，线程 16 会与线程 0 交换值，线程 17 会与线程 1 交换值，依此类推。
当 mask 在后续的迭代中减半（例如 mask = 8, mask = 4, 等等）：
对于 mask = 8，线程 0 会与线程 8 交换值，线程 1 会与线程 9 交换值，依此类推。
对于 mask = 4，线程 0 会与线程 4 交换值，线程 1 会与线程 5 交换值，依此类推。
这个过程继续进行，直到 mask 减少到 1。在 mask = 1 时，每个线程都会与其索引 XOR 1 的线程交换值，即线程 0 与线程 1，线程 2 与线程 3，依此类推。
重要的是要理解，__shfl_xor_sync 并不真正“交换”值；它只是将一个线程的值复制到另一个线程。因此，在调用 __shfl_xor_sync 之后，每个线程都会有一个新值，这个新值可能来自它自己，也可能来自 warp 中的另一个线程。
XOR shuffle 是一种在 CUDA 编程中使用的数据交换技术，特别适用于 GPU 上的 warp 级别的并行操作。在 XOR shuffle 中，“XOR” 表示异或（exclusive OR）操作，而 “shuffle” 则指数据在 warp 中的线程之间被重新排列或交换。
具体来说，XOR shuffle 使用异或运算来确定哪些线程之间应该交换数据。给定一个 mask 值和一个线程的索引 tid，通过计算 tid ^ mask，可以得到一个新的索引，该索引对应于应该与当前线程交换数据的“伙伴”线程。由于 XOR 运算的性质，这个操作是对称的：如果线程 A 的索引 XOR mask 得到线程 B 的索引，那么线程 B 的索引 XOR mask 也将得到线程 A 的索引。这保证了数据交换的一致性。
2.3 Block sum
template<const int kWarpSize = WARP_SIZE>
__device__ __forceinline__ float warp_reduce_sum(float val) {
#pragma unroll
for (int mask = kWarpSize >> 1; mask >= 1; mask >>= 1) //p79，折半归约，位操作比整数操作更高效，相当于/2
{
val += __shfl_xor_sync(0xffffffff, val, mask);//通过线程数洗牌函数 __shfl_xor_sync异或交换数据p111
}
return val;
}
template<const int NUM_THREADS=128>
__device__ __forceinline__ float block_reduce_sum(float val) {
// always <= 32 warps per block (limited by 1024 threads per block)
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE; //整除，向上取整
int warp = threadIdx.x / WARP_SIZE; //threadIdx.x 是当前线程在所在线程块中的索引，除自然等于第几个warp
int lane = threadIdx.x % WARP_SIZE;
static __shared__ float shared[NUM_WARPS];//静态共享内存数组shared[NUM_WARPS]，用于存储每个warp的局部求和结果。static要不要都可、
//以？？？
val = warp_reduce_sum<WARP_SIZE>(val);
if (lane == 0) shared[warp] = val; //lane是上面当前warp的第几个thread，如果是第0个。假设warp=0时，val给shared[0]，总共//
//有num_warps,所以是可以覆盖完的，而且上面那个warp_reduce_sum例子中，最后每个thread的数字都是一样的，所以无所谓取那个
//把每个warp的值给share。
//第几个warp的局部结果都被0号lane存入了共享内存数组shared[warp]中
__syncthreads();//所有线程完成进行下一步
val = (lane < NUM_WARPS) ? shared[lane] : 0.0f;//共享内存声明周期是核函数内，要取出来，lan只是个变量，也就是用lan去除num_warps
//个val出来，定义一个其他的i也可以其实。
val = warp_reduce_sum<NUM_WARPS>(val);
return val;
}
1）计算线程块中的线程束（Warp）数量
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE;
C++11 标准中，定义常量或常量表达式时可以用constexpr修饰，从而使该变量获得在编译阶段即可计算出结果的能力。
通过将线程块中的线程总数（NUM_THREADS）除以线程束大小（WARP_SIZE），向上取整（通过加WARP_SIZE - 1再除以WARP_SIZE实现），得到线程块中包含的线程束数量（NUM_WARPS）（这个解释很好！！！）
2）计算当前线程在warp中的索引（warp）和warp内的索引（lane）
int warp = threadIdx.x / WARP_SIZE;
int lane = threadIdx.x % WARP_SIZE;
第几个warp，这个warp中的第几个lan
3）声明并初始化共享内存
static __shared__ float shared[NUM_WARPS];
声明一个大小为NUM_WARPS的静态共享内存数组shared，用于暂存每个线程束的累加和。
4）在当前线程束（warp）内进行累加求和
val = warp_reduce_sum<WARP_SIZE>(val);
每个线程首先使用warp_reduce_sum函数在其所在的warp内对val进行归约求和。这样，每个warp的第一个线程（lane索引为0的线程）将持有warp内所有线程val值的和。
5）写入共享内存
if (lane == 0) shared[warp] = val;
只有线程束内的第一个线程（lane == 0）将线程束累加和存入共享内存shared[warp]中。
warp是第几个warp的索引，num_warp是有几个warp，lane==1也是可以的，看上面warp_sum的解释。
6）同步线程块
__syncthreads();
使用__syncthreads()指令确保所有线程完成线程束内求和并将结果写入共享内存后，再继续执行后续操作。（线程内完成，然后进行线程间的运算，如下所示）
7）加载相邻线程束的累加和
val = (lane < NUM_WARPS) ? shared[lane] : 0.0f;
(lane < NUM_WARPS)，检查当前线程索引（lane）是否小于线程块中线程束的数量（NUM_WARPS）。如果成立（即true），说明当前线程索引对应的有效线程束（0到NUM_WARPS - 1）；如果不成立（即false），说明当前线程索引超出了有效线程束范围。
对于有效线程束中的线程，赋予它们所对应的线程束的累加和；对于无效线程束中的线程（通常是线程束首领之后的线程），赋予0.0f，确保它们在后续求和时不参与计算。
shared[0]给val，shared[1]给val，shared[2]给val,那么有好多给val，就是num_warps个。
8）在所有线程束之间进行累加求和
线程中
val = warp_reduce_sum<WARP_SIZE>(val);
vs
线程间
val = warp_reduce_sum<NUM_WARPS>(val);
这种两级归约（先线程束内，后线程束间）的方法利用了CUDA硬件的特性，提高了大规模并行计算的效率。这种分级的归约方法可以有效减少线程之间的通信开销，提高归约操作的性能。
问题：解释为何在val = warp_reduce_sum<WARP_SIZE>(val);之后还需要执行val = (lane < NUM_WARPS) ? shared[lane] : 0.0f;这一步。
在完成第一步（线程束内求和）后，每个线程束首领（lane == 0）的val变量包含了本线程束内所有线程的val之和。然而，其他线程（非首领线程）的val变量仍保留着线程束内求和之前的值，这些值在下一步线程束间求和中是不需要的。因此，需要进行如下操作：
线程束首领将本线程束的累加和写入共享内存。
所有线程从共享内存中读取一个适当的值，用于参与下一步的线程束间求和。对于有效线程束中的线程，读取它们所对应的线程束的累加和（都要读取）（即shared[lane]）；对于无效线程束中的线程，读取0.0f。
通过执行val = (lane < NUM_WARPS) ? shared[lane] : 0.0f;这一步，实现了上述目标。这样，每个线程都准备好了参与下一步的线程束间求和所需的数据，最终通过再次调用warp_reduce_sum函数模板完成整个线程块的累加求和。
两个效率比较
`warp_reduce_sum`和`block_reduce_sum`在效率上的比较取决于具体的硬件环境、数据规模以及代码实现细节。然而，通常情况下，`warp_reduce_sum`的效率要高于`block_reduce_sum`，原因如下：
1. **并行度**：`warp_reduce_sum`在单个线程束（Warp）内部进行操作，而线程束是CUDA硬件调度的基本单位，其内部的所有线程能够同时执行同一指令。这意味着在`warp_reduce_sum`中进行的并行累加求和操作能够充分利用硬件的并行能力，没有额外的同步开销。相比之下，`block_reduce_sum`涉及到多个线程束之间的协作，需要使用共享内存进行数据交换和同步，这会引入额外的同步成本和潜在的内存访问瓶颈。
2. **内存访问**：`warp_reduce_sum`仅依赖于寄存器和线程束内的隐式同步，不需要访问全局或共享内存。这种内存访问模式非常高效，因为寄存器访问速度远快于全局或共享内存。相反，`block_reduce_sum`在进行线程束间求和时必须使用共享内存来暂存线程束的累加和，这会增加内存访问的复杂性和潜在延迟。
3. **指令调度**：`warp_reduce_sum`通常使用循环展开和CUDA内建函数（如`__shfl_xor_sync`）来实现高效的数据交换和累加。这些操作通常能被GPU硬件很好地优化和并行执行。而`block_reduce_sum`中的线程束间求和需要额外的循环和条件判断，可能会影响指令流水线的效率。
综上所述，`warp_reduce_sum`由于其更高的并行度、更优的内存访问模式以及更利于硬件优化的指令调度，通常比`block_reduce_sum`具有更高的效率。然而，实际情况可能会受到多种因素的影响，如硬件架构、数据分布、内存访问模式优化等因素。在某些特定场景下，比如数据规模较小、线程块内的线程束数量有限时，两者之间的效率差异可能不会特别显著。在实际应用中，应根据具体任务需求和硬件特性选择合适的归约策略，并通过性能分析和优化来进一步提升效率。
2.4 SGEMV warp k32
单精度浮点数矩阵A与向量x的矩阵-向量乘法，且K为32的倍数的情况。
a: MxK, x: Kx1, y: Mx1, compute: y = a * x
总共32M个线程，一个warp（32个线程）生成一个y的元素，也即是负责a的一行（k控制的），b的一列（k控制的），y的一个元素/一行（k控制的）。
对于n=32的情况，我们将每个block设置为256个线程，4个warp，然后每个warp负责一行元素的计算。每个warp要对x进行访问，然后在warp内部进行一次reduce求和操作。
// SGEMV: Warp SGEMV K32
// 假设K为32的倍数，每个warp负责一行
// grid(M/4), block(32,4) blockDim.x=32=K, blockDim.y=4
// a: MxK, x: Kx1, y: Mx1, compute: y = a * x
__global__ void sgemv_k32(float* a, float* x, float* y, int M, int K) {
int tx = threadIdx.x; // 0~31
int ty = threadIdx.y; // 0~4
int bx = blockIdx.x; // 0~M/4
int lane = tx % WARP_SIZE; // 0~31
int m = bx * blockDim.y + ty; // (0~M/4) * 4 + (0~3)
if (m < M) { //这个条件判断用于确保线程不会访问超出矩阵a的行数。如果m小于M
float sum = 0.0f;
int NUM_WARPS = (K + WARP_SIZE - 1) / WARP_SIZE;
#pragma unroll
for (int w = 0; w < NUM_WARPS; ++w)//对每个warp进行遍历
{
// 若NUM_WARPS>=2，先将当前行的数据累加到第一个warp中
int k = w * WARP_SIZE + lane; //k代表当前warp中的全部索引
sum += a[m * K + k] * x[k]; //假设 a 是一个 4x8 的矩阵（即 4 行 8 列），m 是 2，k 是 3。那么 a[m * K + k] 就等价于 a[2 * 8 +
//3]，即访问矩阵 a 中第 2 行第 3 列的元素。
}
sum = warp_reduce_sum<WARP_SIZE>(sum);
if (lane == 0) y[m] = sum;
}
}
1）索引
int tx = threadIdx.x; // 0~31
int ty = threadIdx.y; // 0~4
int bx = blockIdx.x; // 0~M/4
int lane = tx % WARP_SIZE; // 0~31
int m = bx * blockDim.y + ty; // (0~M/4) * 4 + (0~3)
tx和ty分别表示当前线程在二维线程块内的横纵坐标（线程索引）。
线程块的大小是(32, 4)，因此tx的范围是0到31，ty的范围是0到3。
bx表示当前线程块在二维网格中的横坐标（blockIdx）。
由于网格的大小是(M/4, 1)（假设每个线程块处理4行），bx的范围是0到M/4 - 1。
lane表示当前线程在**所在线程束（Warp）**内的索引（Lane）。
m计算出当前线程对应的在整个矩阵A的行索引，用于访问矩阵A的相应行元素，或最后赋值y的行元素
2）NUM_WARPS
int NUM_WARPS = (K + WARP_SIZE - 1) / WARP_SIZE;
总数/每个warp数据量=warp数
3）warp循环
for (int w = 0; w < NUM_WARPS; ++w) {
int k = w * WARP_SIZE + lane;
sum += a[m * K + k] * x[k];
}
循环变量w从0遍历到NUM_WARPS - 1，对于每个warp，计算对应的k值，这个k值表示在当前warp中当前线程处理的列索引。
lane变量表示当前线程在其warp中的位置，范围通常是0到WARP_SIZE - 1。因此，w * WARP_SIZE + lane这个计算用于确定当前线程处理的矩阵列的具体索引。
在循环体内，代码执行矩阵a的行m与向量x的对应列k的乘积，并将结果累加到sum变量中。这样，每个warp中的线程都负责累加它们各自处理的列到sum中。
4）warp内规约
sum = warp_reduce_sum<WARP_SIZE>(sum);
这一步是将warp内所有线程的sum值相加，得到warp内所有线程计算的点积总和。
2.5 SGEMV K128 float4
对于n=128的情况，同样让warp负责一行元素的计算，但是因为每行的元素比较多，所以采用了float4进行向量化的访存。能够有更高的访存效率。
针对的是K为128的倍数的情况，并且使用了float4类型来加速内存访问和计算。float4是一个包含四个单精度浮点数的向量类型，它允许一次从全局内存中读取或写入四个连续的浮点数，从而提高了内存带宽利用率。
// SGEMV: Warp SGEMV K128 + Vec4
// 假设K为128的倍数 float4
// grid(M/4), block(32,4) blockDim.x=32=K, blockDim.y=4
// a: MxK, x: Kx1, y: Mx1, compute: y = a * x
__global__ void sgemv_k128(float* a, float* x, float* y, int M, int K) {
// 每个线程负责4个元素，一个warp覆盖128个元素
int tx = threadIdx.x; // 0~31
int ty = threadIdx.y; // 0~3
int bx = blockIdx.x; // 0~M/4
int lane = tx % WARP_SIZE; // 0~31
int m = blockDim.y * bx + ty; // (0~M/4) * 4 + (0~3)
if (m < M) {
float sum = 0.0f;
// process 4*WARP_SIZE elements per warp.
int NUM_WARPS = (((K + WARP_SIZE - 1) / WARP_SIZE) + 4 - 1) / 4;
#pragma unroll
for (int w = 0; w < NUM_WARPS; ++w) {
int k = (w * WARP_SIZE + lane) * 4; //当前线程索引k，负责四个元素，因此要×4，即四个连续浮点数单精度的数据存储在索引k的线程中。
float4 reg_x = FLOAT4(x[k]); //向量 float4来封装四个数，赋给reg_x，上面一个k代表4个数
float4 reg_a = FLOAT4(a[m * K + k]); //矩阵
sum += (reg_a.x * reg_x.x + reg_a.y * reg_x.y
+ reg_a.z * reg_x.z + reg_a.w * reg_x.w);
}
sum = warp_reduce_sum<WARP_SIZE>(sum);
if(lane == 0) y[m] = sum;
}
}
1）warp 数量
读取矩阵x需要多少warp
int NUM_WARPS = (K + WARP_SIZE - 1) / WARP_SIZE;
vs
int NUM_WARPS = (((K + WARP_SIZE - 1) / WARP_SIZE) + 4 - 1) / 4;
int NUM_WARPS = (K + (WARP_SIZE * 4 - 1)) / (WARP_SIZE * 4); //这样或许更好
所需线程束数量 = 总数据量 / 每个线程束处理的数据量，一般需要向上取整。
2）K
int k = (w * WARP_SIZE + lane) * 4;
乘以4的原因在于float4是一个四分量的矢量类型，它一次性存储并处理四个连续的单精度浮点数。因此，要正确访问向量x中对应的float4数据，索引必须对应这四个连续元素的起始位置。
我理解您可能是想问，在日常编程中，如果我们只需要保存一个单精度浮点数（而不是`float4`这类包含四个单精度浮点数的矢量类型），那么是否还需要乘以4来计算索引？
答案是**不需要**。在日常编程中，如果我们只需要保存和操作一个单精度浮点数，那么直接使用该浮点数的索引来访问和存储就足够了，无需乘以4。乘以4的操作是针对使用`float4`这类矢量类型时，为了正确访问和操作由四个连续单精度浮点数组成的矢量而进行的特殊处理。
举个例子，假设有一个数组`float singlePrecisionArray[]`，其中每个元素都是单精度浮点数。如果我们想访问第`i`个元素，直接使用索引`i`即可：
float value = singlePrecisionArray[i];
在这里，无需对索引`i`乘以4，因为数组中的每个元素本身就是单精度浮点数，不需要转换为`float4`的索引。
总结来说，乘以4的操作仅适用于需要处理`float4`这类包含多个单精度浮点数的矢量类型的情况，目的是为了正确访问和操作由多个连续单精度浮点数组成的矢量。而在日常编程中，如果我们只需要保存和操作单个单精度浮点数，直接使用对应的索引即可，无需乘以4
3）使用 float4 加载数据
float4 reg_x = FLOAT4(x[k]);
float4 reg_a = FLOAT4(a[m * K + k]);
这两行分别读取向量x和矩阵a的对应元素，并将它们存储在float4类型的变量中，以加速内存访问。
这行代码创建了一个float4类型的变量reg_x，并初始化它为向量x的第k个元素。这里假设x是一个包含多个浮点数的数组或向量，而k是当前的列索引（k里面包含了4个）。就有下面的
sum += (reg_a.x * reg_x.x + reg_a.y * reg_x.y
+ reg_a.z * reg_x.z + reg_a.w * reg_x.w);
4） #pragma unroll
这是一个CUDA编译指令，它告诉编译器在编译时展开循环，而不是在运行时动态执行循环。这通常可以提高性能，因为减少了循环控制相关的开销。
5）其它
在上述CUDA内核函数`sgemv_k128`中，对`k`的计算中乘以`4`的原因在于`float4`是一个四分量的矢量类型，它一次性存储并处理四个连续的单精度浮点数。因此，要正确访问向量`x`中对应的`float4`数据，索引必须对应这四个连续元素的起始位置。
当计算`k = (w * WARP_SIZE + lane)`时，这个表达式给出的是当前线程在矩阵`A`中对应列的索引。然而，`A`中的每个元素是单精度浮点数，而我们实际上需要加载和操作的是`float4`类型的向量元素。所以，为了得到向量`x`中对应`float4`的起始地址，需要将这个单精度浮点数的索引转换为`float4`元素的索引。
由于每个`float4`包含4个单精度浮点数，要定位到`float4`的起始位置，就需要将计算出的单精度浮点数索引乘以`4`。这样做确保了：
- 当使用`FLOAT4(x[k])`来加载向量`x`时，会正确地提取出由四个连续单精度浮点数组成的`float4`数据。
- 在随后的乘加运算中，如`reg_a.x * reg_x.x + reg_a.y * reg_x.y + reg_a.z * reg_x.z + reg_a.w * reg_x.w`，每个`reg_a`的分量与对应的`reg_x`分量相乘，并累加，这些操作都是基于`float4`类型的向量进行的。
总结来说，对`k`乘以`4`是为了将原本计算出的单精度浮点数列索引转换为适用于`float4`数据类型的操作索引，以便正确地访问和处理向量`x`中的`float4`元素，实现高效的SIMD计算。
6） 并行，不需要m++
if (m < M)
只要比M小，一直都并行。
2.6 dot product
计算两个同型向量a和b的点积
// Dot Product
// grid(N/128), block(128)
// a: Nx1, b: Nx1, y=sum(elementwise_mul(a,b))
template<const int NUM_THREADS = 128>
__global__ void dot(float* a, float* b, float* y, int N) {
int tid = threadIdx.x;
int idx = blockIdx.x * NUM_THREADS + tid;
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE;
__shared__ float reduce_smem[NUM_WARPS];
// keep the data in register is enougth for warp operaion.
float prod = (idx < N) ? a[idx] * b[idx] : 0.0f;
int warp = tid / WARP_SIZE;
int lane = tid % WARP_SIZE;
// perform warp sync reduce.
prod = warp_reduce_sum<WARP_SIZE>(prod);
// warp leaders store the data to shared memory.
if (lane == 0) reduce_smem[warp] = prod;
__syncthreads(); // make sure the data is in shared memory.
// the first warp compute the final sum.
prod = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
if (warp == 0) prod = warp_reduce_sum<NUM_WARPS>(prod);
if (tid == 0) atomicAdd(y, prod);
}
1）
__shared__ float reduce_smem[NUM_WARPS];
声明了一个名为 reduce_smem 的数组，其类型为 float，数组长度为 NUM_WARPS。这意味着数组 reduce_smem 包含 NUM_WARPS 个连续的 float 类型元素。
2）点积计算与Warp级归约
// keep the data in register is enougth for warp operaion.
float prod = (idx < N) ? a[idx] * b[idx] : 0.0f;
int warp = tid / WARP_SIZE;
int lane = tid % WARP_SIZE;
// perform warp sync reduce.
prod = warp_reduce_sum<WARP_SIZE>(prod);
当前线程idx计算其负责的向量元素对a[idx]和b[idx]的乘积，并存储在prod中。
根据tid计算当前线程所在的Warp索引（warp）和在该Warp内的线程索引（lane）。
使用warp_reduce_sum模板函数对prod进行Warp级别的局部求和。
3）共享内存写入与同步
if (lane == 0) reduce_smem[warp] = prod;
__syncthreads(); // make sure the data is in shared memory.
线程索引 lane 为0时，表示当前线程是其所在 Warp 的第一个线程，也被称为 Warp 的“首领线程”。写入共享内存的相应位置。
只有每个Warp的第一个线程（lane == 0）将其局部求和结果存入共享内存的相应位置。然后使用__syncthreads()确保所有线程都完成了写入操作。
4）最终全局归约
prod = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
if (warp == 0) prod = warp_reduce_sum<NUM_WARPS>(prod);
接下来，每个线程读取与其所在Warp对应的共享内存中的中间结果（对于非第一个Warp的线程，这步无意义）。然后，仅由第一个Warp（warp == 0）再次使用warp_reduce_sum对所有Warp的中间结果进行全局归约，得到最终的点积结果。
5）原子性地更新输出
if (tid == 0) atomicAdd(y, prod);
atomicAdd 被用于将最终计算得到的点积结果（存储在 prod 变量中）原子性地累加到全局内存中的目标变量 y。
最后，仅由第一个线程块的第一个线程（tid == 0）使用atomicAdd函数将最终点积结果原子性地累加到全局输出变量y上。
2.7 softmax
这段代码定义了一个名为softmax的CUDA内核函数，用于计算给定向量x的softmax值，并将结果存储在向量y中。向量x和y的长度均为N。
// Softmax x: N, y: N
// grid(N/128), block(K=128)
template<const int NUM_THREADS = 128>
__global__ void softmax(float* x, float* y, float* total, int N) {
const int tid = threadIdx.x;
const int idx = blockIdx.x * blockDim.x + tid;
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE;
__shared__ float reduce_smem[NUM_WARPS];
float sum = (idx < N) ? expf(x[idx]) : 0.0f;
int warp = tid / WARP_SIZE;
int lane = tid % WARP_SIZE;
sum = warp_reduce_sum<WARP_SIZE>(sum);
if (lane == 0) reduce_smem[warp] = sum;
__syncthreads();
// compute the final sum in each warp
sum = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
sum = warp_reduce_sum<NUM_WARPS>(sum); // sum(e^x_0,...,e^x_n-1)
// get the total sum of all blocks.
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
// e^x_i/sum(e^x_0,...,e^x_n-1)
if (idx < N) y[idx] = expf(x[idx]) / (*total);
}
1）Softmax计算与Warp级归约
sum为每个warp的指数和
float sum = (idx < N) ? expf(x[idx]) : 0.0f;
int warp = tid / WARP_SIZE;
int lane = tid % WARP_SIZE;
sum = warp_reduce_sum<WARP_SIZE>(sum);
2）模板
if (lane == 0) reduce_smem[warp] = sum;
__syncthreads();
// compute the final sum in each warp
sum = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
只有每个Warp的第一个线程（lane == 0）将其局部求和结果存入共享内存的相应位置。然后使用__syncthreads()确保所有线程都完成了写入操作。
存到共享存储，然后从共享存储取值。
如果lane为3，该线程读取reduce_smem[3]（即Warp 3的中间结果sum_3）并赋值给sum。
为了更好地解释这段代码的工作原理，我们假设一个具体的例子：
- 线程块大小（`blockDim.x`）为128（即`NUM_THREADS`），Warp大小（`WARP_SIZE`）为32。
- 假设当前线程块包含4个Warp（因为`NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE = (128 + 32 - 1) / 32 = 4`）。
**局部归约结果写入共享内存**：
- 假设当前线程是Warp 0的首领线程（即`lane`为0），其局部归约结果为`sum_0`。此时，该线程将`sum_0`写入`reduce_smem`的索引为0的位置，即`reduce_smem[0] = sum_0`。
- 同理，Warp 1、Warp 2和Warp 3的首领线程分别将它们的局部归约结果`sum_1`、`sum_2`、`sum_3`写入`reduce_smem`的索引为1、2、3的位置。
**同步线程**：
- 执行`__syncthreads()`，确保所有线程都完成了上述写入操作，共享内存中的数据已经准备好供后续读取。
**读取Warp中间结果**：
- 假设当前线程属于Warp 0（即`lane`值小于`NUM_WARPS = 4`），其`lane`值为`lane_0`（0到31之间的一个整数）。
- 如果`lane_0`为0，该线程读取`reduce_smem[0]`（即Warp 0的中间结果`sum_0`）并赋值给`sum`。
- 如果`lane_0`为1，该线程读取`reduce_smem[1]`（即Warp 1的中间结果`sum_1`）并赋值给`sum`。
- ...
- 如果`lane_0`为3，该线程读取`reduce_smem[3]`（即Warp 3的中间结果`sum_3`）并赋值给`sum`。
- 同理，Warp 0的其他线程（`lane`值从1到31）也会依次读取所有Warp的中间结果，并将它们存储在各自的`sum`变量中。
总结：在这个例子中，通过上述三步操作，Warp 0的所有线程（包括首领线程和其他线程）都将各自负责的Warp的中间结果存储在`sum`变量中，为后续的跨Warp归约或进一步计算做好准备。其他Warp（如Warp 1、Warp 2、Warp 3）的线程也会执行类似的操作。这样，通过共享内存和线程同步，实现了在每个线程块内高效地收集和传递各个Warp的中间结果。
3）全局归约
sum = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
sum = warp_reduce_sum<NUM_WARPS>(sum); // sum(e^x_0,...,e^x_n-1)
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
接着，仅由第一个线程（tid == 0）使用atomicAdd函数将当前线程块的求和结果原子性地累加到全局变量total上。
2.8 更简单的softmax
// Softmax x: N, y: N
// grid(N/128), block(K=128)
template<const int NUM_THREADS = 128>
__global__ void softmax_v2(float* x, float* y, float* total, int N) {
const int tid = threadIdx.x;
const int idx = blockIdx.x * blockDim.x + tid;
float exp_val = (idx < N) ? expf(x[idx]) : 0.0f;
float sum = block_reduce_sum<NUM_THREADS>(exp_val);
// get the total sum of all blocks.
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
// e^x_i/sum(e^x_0,...,e^x_n-1)
if (idx < N) y[idx] = exp_val / (*total);
}
2.8.1 解释版
template<const int kWarpSize = WARP_SIZE>
__device__ __forceinline__ float warp_reduce_sum(float val)
{
#pragma unroll
for(int mask= kWarpSize>>1; mask>=1;mask>>1)
{
//tid与mask按位与，得到目标线程，然后与现在的线程相加，得到新的val
val += __shfl_xor_sync(0xffffffff,val,mask);
}
return val;
}
template<const int NUM_THREADS=128>
__device__ __forceinline__ float block_reduce_sum(float val)
{
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1)/WARP_SIZE;
int warp = threadIdx.x/WARP_SIZE; //num_warp是有几个，warp是索引
int lane = threadIdx.x%WARP_SIZE; //Warp内的索引，一个warp有32个线程
__shared__ float shared[NUM_WARPS];
//warp内
val = warp_reduce_sum<WARP_SIZE>(val);
if(lane==0) shared[warp] = val; //如csdn里面一样，最后所有tid都是一个结果，所以随便哪个都可以
__syncthreads();//保证块内所有线程完成以上操作
//warp间
val = (lane<NUM_WARPS) ? shared[lane]:0.0f;
val = warp_reduce_sum<NUM_WARPS>(val);
return val;
}
//softmax
template<const int NUM_THREADS = 128>
__global__ void softmax(float* x, float* y,float* total,int N)
{
const int tid = threadIdx.x;
const int idx = blockDim.x*blockIdx.x + tid;
// 分子
float exp_val = (idx<N)? expf(x[idx]):0.0f;
// 分母,block级的sum
float sum = block_reduce_sum<NUM_THREADS>(exp_val);
if(tid==0) atomicAdd(total,sum);
__threadfence(); //块间所有线程,否则拿不到全局的exp sum作为分母项。
if(idx<N) y[idx] = exp_val / (*total);
}
2.8.1 softmax代码小结
// Warp Reduce Sum
template<const int kWarpSize = WARP_SIZE>
__device__ __forceinline__ float warp_reduce_sum(float val) {
#pragma unroll
for (int mask = kWarpSize >> 1; mask >= 1; mask >>= 1) {
val += __shfl_xor_sync(0xffffffff, val, mask);
}
return val;
}
// Softmax x: N, y: N
// grid(N/128), block(K=128)
template<const int NUM_THREADS = 128>
__global__ void softmax(float* x, float* y, float* total, int N) {
const int tid = threadIdx.x;
const int idx = blockIdx.x * blockDim.x + tid;
constexpr int NUM_WARPS = (NUM_THREADS + WARP_SIZE - 1) / WARP_SIZE;
__shared__ float reduce_smem[NUM_WARPS];
float sum = (idx < N) ? expf(x[idx]) : 0.0f;
int warp = tid / WARP_SIZE;
int lane = tid % WARP_SIZE;
sum = warp_reduce_sum<WARP_SIZE>(sum);
if (lane == 0) reduce_smem[warp] = sum;
__syncthreads();
// compute the final sum in each warp
sum = (lane < NUM_WARPS) ? reduce_smem[lane] : 0.0f;
sum = warp_reduce_sum<NUM_WARPS>(sum); // sum(e^x_0,...,e^x_n-1)
// get the total sum of all blocks.
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
// e^x_i/sum(e^x_0,...,e^x_n-1)
if (idx < N) y[idx] = block_smem[tid] / (*total);
}
// Softmax x: N, y: N
// grid(N/128), block(K=128)
template<const int NUM_THREADS = 128>
__global__ void softmax_v2(float* x, float* y, float* total, int N) {
const int tid = threadIdx.x;
const int idx = blockIdx.x * blockDim.x + tid;
float exp_val = (idx < N) ? expf(x[idx]) : 0.0f;
float sum = block_reduce_sum<NUM_THREADS>(exp_val);
// get the total sum of all blocks.
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
// e^x_i/sum(e^x_0,...,e^x_n-1)
if (idx < N) y[idx] = exp_val / (*total);
}
// Softmax Vec4 x: N, y: N
// grid(N/128), block(128/4)
template<const int NUM_THREADS = 128/4>
__global__ void softmax_v2_vec4(float* x, float* y, float* total, int N) {
const int tid = threadIdx.x;
const int idx = (blockIdx.x * blockDim.x + tid) * 4;
float4 reg_x = FLOAT4(x[idx]);
float4 reg_exp;
reg_exp.x = (idx < N) ? expf(reg_x.x) : 0.0f;
reg_exp.y = (idx < N) ? expf(reg_x.y) : 0.0f;
reg_exp.z = (idx < N) ? expf(reg_x.z) : 0.0f;
reg_exp.w = (idx < N) ? expf(reg_x.w) : 0.0f;
float exp_val = (reg_exp.x + reg_exp.y + reg_exp.z + reg_exp.w);
float sum = block_reduce_sum<NUM_THREADS>(exp_val);
// get the total sum of all blocks.
if (tid == 0) atomicAdd(total, sum);
__threadfence(); // grid level memory fence 注意这里需要网格级别的内存同步
// e^x_i/sum(e^x_0,...,e^x_n-1)
if (idx < N) {
float4 reg_y;
reg_y.x = reg_exp.x / (*total);
reg_y.y = reg_exp.y / (*total);
reg_y.z = reg_exp.z / (*total);
reg_y.w = reg_exp.w / (*total);
FLOAT4(y[idx]) = reg_y;
}
}
2.8.2 矩阵乘代码小结
不同维度的thread、block的求法
简洁明了点
参考图更清晰
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <vector>
#include <algorithm>
#include <cuda_runtime.h>
#define WARP_SIZE 32
#define INT4(value) (reinterpret_cast<int4*>(&(value))[0])
#define FLOAT4(value) (reinterpret_cast<float4*>(&(value))[0])
// SGEMM: Block Tile + K Tile, with smem
// Block Tile (BM, BN) + K Tile (BK=32)
// grid((N + BN - 1) / BN, (M + BM - 1) / BM), block(BN, BM)
// a: MxK, b: KxN, c: MxN, compute: c = a * b, all row major
__global__ void sgemm(float* a, float* b, float* c, int M, int N, int K) {
// [1] Block Tile: 32x32的block处理c上一块32x32的元素计算
// [2] K Tile: 使用共享内存，并将K分块为BK大小的块
constexpr int BM = 32;
constexpr int BN = 32;
constexpr int BK = 32;
__shared__ float s_a[BM][BK], s_b[BK][BN];
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int tid = threadIdx.y * blockDim.x + tx; // tid within the block
// load values to shared memory, 32x32 threads working together
// to fetch data along the row direction of a and b both for s_a
// and s_b 32x32x4x2=8KB, we use 32x32 threads within block to
// load 32x32 elements from global memory to shared memory, namely,
// each thread will load 1 element.
int load_smem_a_m = tid / 32; // 0~31, tid / 32, tid / BM, threadIdx.y
int load_smem_a_k = tid % 32; // 0~31, tid % 32, tid % BK, threadIdx.x
int load_smem_b_k = tid / 32; // 0~31, tid / 32, tid / BK, threadIdx.y
int load_smem_b_n = tid % 32; // 0~31, tid % 32, tid % BN, threadIdx.x
int load_gmem_a_m = by * BM + load_smem_a_m; // global row of a and c
int load_gmem_b_n = bx * BN + load_smem_b_n; // global col of b and c
// if (load_gmem_a_m >= M || load_gmem_b_n >= N) return;
float sum = 0.f;
for (int bk = 0; bk < (K + BK - 1) / BK; ++bk) {
int load_gmem_a_k = bk * BK + load_smem_a_k;
int load_gmem_a_addr = load_gmem_a_m * K + load_gmem_a_k;
s_a[load_smem_a_m][load_smem_a_k] = a[load_gmem_a_addr];
int load_gmem_b_k = bk * BK + load_smem_b_k;
int load_gmem_b_addr = load_gmem_b_k * N + load_gmem_b_n;
s_b[load_smem_b_k][load_smem_b_n] = b[load_gmem_b_addr];
__syncthreads();
#pragma unroll
for (int k = 0; k < BK; ++k) {
int comp_smem_a_m = load_smem_a_m;
int comp_smem_b_n = load_smem_b_n;
sum += s_a[comp_smem_a_m][k] * s_b[k][comp_smem_b_n];
}
__syncthreads();
}
int store_gmem_c_m = load_gmem_a_m;
int store_gmem_c_n = load_gmem_b_n;
int store_gmem_c_addr = store_gmem_c_m * N + store_gmem_c_n;
c[store_gmem_c_addr] = sum;
}
// SGEMM: Block Tile + Thread Tile + K Tile + Vec4, with smem
// BK:TILE_K=8 BM=BN=128
// TM=TN=8 增加计算密度 BM/TM=16 BN/TN=16
// dim3 blockDim(BN/TN, BM/TM);
// dim3 gridDim((N + BN - 1) / BN, (M + BM - 1) / BM)
__global__ void sgemm_thread_tile_vec4(
float* a, float* b, float* c, int M, int N, int K) {
// [1] Block Tile: 一个16x16的block处理C上大小为128X128的一个目标块
// [2] Thread Tile: 每个thread负责计算TM*TN(8*8)个元素，增加计算密度
// [3] K Tile: 将K分块，每块BK大小，迭代(K+BK-1/BK)次，
// 每次计算TM*TN个元素各自的部分乘累加
// [4] Vectorize: 减少load和store指令，使用float4
constexpr int BM = 128;
constexpr int BN = 128;
constexpr int BK = 8;
constexpr int TM = 8;
constexpr int TN = 8;
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int tid = threadIdx.y * blockDim.x + tx; // tid within the block
__shared__ float s_a[BM][BK], s_b[BK][BN]; // 2*128*8*4=8KB
// 0. 先计算shared memory中的索引
// tid和需要加载的smem s_a[BM][BK] 之间的索引关系 BM=128 BK=8 按行读取 A行主序
// 对于s_a每行8个数据，每个线程读取4个，需要2个线程；总共128行，需要128x2刚好256线程
int load_smem_a_m = tid / 2; // tid/2 (128/8)*(128/8)=256 threads per block, tid/2->[0,128), BM=128 0~127
int load_smem_a_k = (tid % 2 == 0) ? 0 : 4; // (tid%2 == 0) ? 0 : 4, col of s_a 0,4
// tid和需要加载的smem s_b[BK][BN] 之间的索引关系 BK=8 BN=128 按行读取 B行主序
// 对于s_b每行128个数据，每个线程读4个数据，需要32个线程；总共8行，需要32x8=256个线程
int load_smem_b_k = tid / 32; // tid/32, row of s_b 256/32=8 行 0~7
int load_smem_b_n = (tid % 32) * 4; // (tid % 32) * 4, col of s_b 0,4,...,124
// 1. 再计算全局内存中的索引
// 要加载到s_a中的元素对应到A全局内存中的行数 每个block负责出C中大小为BM*BN的块
int load_gmem_a_m = by * BM + load_smem_a_m; // global row of a and c
int load_gmem_b_n = bx * BN + load_smem_b_n; // global col of b and c
float r_c[TM][TN] = {0.0}; // 8x8
// 2. 先对K进行分块，每块BK大小
for (int bk = 0; bk < (K + BK - 1) / BK; ++bk) {
// 加载数据到共享内存smem s_a BM*BK 128*8 vectorize float4
int load_gmem_a_k = bk * BK + load_smem_a_k; // global col of a
int load_gmem_a_addr = load_gmem_a_m * K + load_gmem_a_k;
FLOAT4(s_a[load_smem_a_m][load_smem_a_k]) = FLOAT4(a[load_gmem_a_addr]);
// 加载数据到共享内存smem s_b BK*BN 8*128 vectorize float4
int load_gmem_b_k = bk * BK + load_smem_b_k; // global row of b
int load_gmem_b_addr = load_gmem_b_k * N + load_gmem_b_n;
FLOAT4(s_b[load_smem_b_k][load_smem_b_n]) = FLOAT4(b[load_gmem_b_addr]);
__syncthreads();
#pragma unroll
for (int k = 0; k < BK; k++) {
// 3. 每个线程负责计算BM*BN(12x128)中的TM*TN(8x8)个元素
#pragma unroll
for (int m = 0; m < TM; m++) {
#pragma unroll
for (int n = 0; n < TN; n++) {
// k from 0~7，0 ~ BK, ty and tx range from 0 to 15, 16x8=128
int comp_smem_a_m = ty * TM + m; // 128*8 128/TM(8)=16 M方向 16线程
int comp_smem_b_n = tx * TN + n; // 8*128 128/TN(8)=16 N方向 16线程
r_c[m][n] += s_a[comp_smem_a_m][k] * s_b[k][comp_smem_b_n];
}
}
}
__syncthreads();
}
#pragma unroll
for (int m = 0; m < TM; ++m) {
int store_gmem_c_m = by * BM + ty * TM + m;
#pragma unroll
for (int n = 0; n < TN; n += 4) {
int store_gmem_c_n = bx * BN + tx * TN + n;
int store_gmem_c_addr = store_gmem_c_m * N + store_gmem_c_n;
FLOAT4(c[store_gmem_c_addr]) = FLOAT4(r_c[m][n]);
}
}
}
2.8.3 矩阵乘解释版
__global__ void sgemm(float* a,float* b,float* c,int M,int N,int K)
{//只运用了share memory，block级别的tile，没有用到thread的再次划分
constexpr BM=32;
constexpr BN=32;
constexpr BK=32;
__shared__ float s_a[BM][BK],s_b[BK][BN];
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
/*
假设我们有一个块，其维度为blockDim.x = 16和blockDim.y = 8。这意味着这个块中有16列（x方向）
和8行（y方向）的线程。如果我们想要唯一地标识这个块中的每个线程，我们需要一个从0到127
（总共128个线程）的整数。
对于第一行的线程（threadIdx.y = 0），它们的ID将是0到15（因为blockDim.x = 16）。
对于第二行的线程（threadIdx.y = 1），它们的ID将是16到31（因为它们的y索引是1，
所以它们的ID是从1*blockDim.x开始的）。
以此类推，我们可以发现每个线程的ID是其y索引乘以块的x维度，再加上其x索引。
这正是int tid = threadIdx.y*blockDim.x + tx;这行代码所做的事情。
*/
int tid = threadIdx.y*blockDim.x + tx;
// 共享内存地址，thread
// 一个warp读取一行,32*32个线程搬运32*32个数
int load_smem_a_m = tid / 32;
int load_smem_a_k = tid % 32;
int load_smem_b_k = tid / 32;
int load_smem_b_n = tid % 32;
//全局内存地址，block
/*
by * BM给出了当前线程块处理的C矩阵的行在A矩阵中的起始位置，
而load_smem_a_m给出了当前线程在块中需要加载的A矩阵的行的相对偏移量。
两者相加，就得到了应该从全局内存a中加载的确切行索引。
*/
int load_gmem_a_m = by * BM + load_gmem_a_m;
int load_gmem_b_n = bx * BN + load_smem_b_n;
float sum = 0.f;
for(int bk=0; bk<(K+BK-1)/BK; bk++)
{
//大循环，A1的纵坐标
int load_gmem_a_k = bk * BK + load_gmem_a_k;
int load_gmem_a_addr = load_gmem_a_m * K + load_gmem_a_k;
//全局给共享
s_a[load_smem_a_m][load_gmem_a_k] = a[load_gmem_a_addr];
//b1
int load_gmem_b_k = bk * BK + load_smem_b_k;
int load_gmem_b_addr = load_gmem_b_k * N + load_gmem_b_n;
s_b[load_smem_b_k][load_smem_b_n] = b[load_gmem_b_addr];
__syncthreads();
#pragma unroll
// 共享内存计算,所以是上步获取的行comp_seme_a_m,列comp_seme_b_n，
// 用k去遍历它
for(int k=0;k<BK;k++)//BK/1
{
int comp_seme_a_m = load_smem_a_m;
int comp_seme_b_n = load_smem_b_n;
sum += s_a[comp_seme_a_m][k] * s_b[k]*[comp_seme_b_n];
}
__syncthreads();
}
int store_gmem_c_m = load_gmem_a_m;
int store_gmem_c_n = load_gmem_b_n;
int store_gmem_c_addr = store_gmem_c_m * N + store_gmem_c_n;
c[store_gmem_c_addr] = sum;
}
__global__ void sgemm_thread_tile_vec4(
float* a,float* b,float* c,int M,int N,int K)
{
//定义块大小
constexpr int BM = 128;
constexpr int BN = 128;
constexpr int BK = 8;
constexpr int TM = 8;
constexpr int TN = 8;
// 索引
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int tid = blockDim.x*threadIdx.y + tx;
__shared__ float s_a[BM][BK],s_b[BK][BN];
//先计算shared memory中的索引,share可以算完
//1%2==1,then 4,col=4,也就是一行两个，要么是0的位置，要么4的位置
int load_smem_a_m = tid / 2;
int load_smem_a_k = (tid % 2 == 0) ? 0:4;
//一行32个，隔4个一个线程
int load_smem_b_k = tid / 32;
int load_smem_b_n = (tid % 32) * 4;
//全局只能算a的行，b的列
int load_gmem_a_m = by * BM + load_gmem_a_m;
int load_gmem_b_n = bx * BN + load_smem_b_n;
//写for，将全局数据给share
//大循环，对应第一次拆分，从全局内存搬数据到共享内存
for(int bk=0; bk<(K+BK-1)/BK; bk++)
{
//A1和B1的列由bk控制，然后算出地址，然后给share
int load_gmem_a_k = bk * BK+ load_gmem_a_k;
int load_gmem_a_addr = load_gmem_a_m * K + load_gmem_a_k;
//Float4读取数据，减少访问
FLOAT4(s_a[load_gmem_a_m][load_gmem_a_k]) = FLOAT4(a[load_gmem_a_addr]);
//B1来
int load_gmem_b_k = bk * BK + load_smem_b_k; // global row of b
int load_gmem_b_addr = load_gmem_b_k * N + load_gmem_b_n;
FLOAT4(s_b[load_smem_b_k][load_smem_b_n]) = FLOAT4(b[load_gmem_b_addr]);
__syncthreads();
#pragma unroll
//开始计算，并保存到reg,再次细分就要多这些。
//第二次拆分，涉及到共享内存中的计算
for(int k=0; k<BK; k++)
{
//先写出TM，TN的for遍历
for(int m=0;m<TM;m++)
{
for(int n =0;n<TN;n++)
{
//由C11的行列倒推A11，B11的行列索引
//C11的行comp_smem_a_m，确实是这样，在共享内存上，当前线程=ty*TM+m
int comp_smem_a_m = ty * TM + m;
int comp_smem_b_n = tx * TN + n;
r_c[m][n] += s_a[comp_smem_a_m][k] * s_b[k][comp_smem_b_n];
}
}
}
__syncthreads();
}
//由于存在更小的分块，则行和列均由3部分构成：全局行号store_c_gmem_m等于
//大分块的起始行号by * BM+小分块的起始行号ty * TM+小分块内部的相对行号i。列同理。
#pragma unroll
for (int m = 0; m < TM; ++m) {
//
int store_gmem_c_m = by * BM + ty * TM + m;
#pragma unroll
for (int n = 0; n < TN; n += 4) {
int store_gmem_c_n = bx * BN + tx * TN + n;
int store_gmem_c_addr = store_gmem_c_m * N + store_gmem_c_n;
FLOAT4(c[store_gmem_c_addr]) = FLOAT4(r_c[m][n]);
}
}
}
2.9 rms_norm
// RMS Norm: x: NxK(K=128<1024), y': NxK, y'=x/rms(x) each row
// 1/rms(x) = rsqrtf( sum(x^2)/K ) each row
// grid(N*K/K), block(K<1024) N=batch_size*seq_len, K=hidden_size
// y=y'*g (g: scale)
template<const int NUM_THREADS=128>
__global__ void rms_norm(float* x, float* y, float g, int N, int K) {
int tid = threadIdx.x; // 0..K-1
int bid = blockIdx.x; // 0..N-1
int idx = bid * blockDim.x + threadIdx.x;
const float epsilon = 1e-5f;
__shared__ float s_variance; // shared within block
float value = (idx < N * K) ? x[idx] : 0.0f; // load once only //rms的左式
float variance = value * value;
variance = block_reduce_sum<NUM_THREADS>(variance); //rms的根号里面的平方求和
if (tid == 0) s_variance = rsqrtf(variance / (float) K + epsilon);
// wait for s_variance in shared memory to be ready for all threads
__syncthreads();
if (idx < N * K) y[idx] = (value * s_variance) * g;
}
1）
const float epsilon = 1e-5f;
__shared__ float s_variance; // shared within block
定义一个常量epsilon用于数值稳定性。接着，声明一个共享内存变量s_variance，用于存储当前线程块内所有数据元素平方和的归一化系数。
2）x: NxK(K=128<1024)
float value = (idx < N * K) ? x[idx] : 0.0f; // load once only
这部分检查当前线程的全局索引idx是否小于矩阵x的总元素数N * K
3）线程块内线程间的全局同步归约
float variance = value * value;
variance = block_reduce_sum<NUM_THREADS>(variance);
4）
variance = block_reduce_sum<NUM_THREADS>(variance); //rms的根号里面的平方求和
if (tid == 0) s_variance = rsqrtf(variance / (float) K + epsilon);
variance 全局和，它相当于公式中的 Σ(x[i]^2)。
(variance / (float) K) 是对平方和进行平均，得到当前线程块内所有元素平方的平均值（均方值），对应公式中的 (1/K) * Σ(x[i]^2)。
+ epsilon 是为了增加数值稳定性，防止除以接近零的数（尤其是在均方值较小的情况下），对应公式中的 ε。
rsqrtf(...) 函数计算给定数值的倒数平方根，相当于公式中的 rsqrt(mean_squared)。
矩阵乘:
矩阵转置: 访存密集型算子的处理
一维reduce-sum：重点是如何处理bank confict
二维reduce-sum
卷积
将单stream改成多stream
6.10 矩阵转置
__global__ void transpose2(const real *A, real *B, const int N)
{
const int nx = blockIdx.x * blockDim.x + threadIdx.x;
const int ny = blockIdx.y * blockDim.y + threadIdx.y;
if (nx < N && ny < N)
{
B[ny * N + nx] = A[nx * N + ny];
}
}
2.10 向量求和
2.10.1普通向量求和
#include <cuda.h>
#include <cuda_runtime.h>
#include <iostream>
__global__ void sumElements(int* d_array, int n, int* result) {
unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < n) {
atomicAdd(result, d_array[tid]); // 每个线程将其负责的元素加到结果上
}
}
int main() {
const int n = 10; // 数组元素数量
int host_array[n];
for (int i = 0; i < n; ++i) host_array[i] = i + 1; // 初始化数组
int host_sum = 0; // 初始化主机端的和为0
int* device_array;
int* device_sum;
cudaMalloc((void**)&device_array, n * sizeof(int));
cudaMalloc((void**)&device_sum, sizeof(int));
cudaMemcpy(device_array, host_array, n * sizeof(int), cudaMemcpyHostToDevice);
cudaMemset(device_sum, 0, sizeof(int)); // 清零设备端的和
dim3 block(1); // 假设每个block只有一个线程（简化示例）
dim3 grid((n + block.x - 1) / block.x); // 确保有足够多的线程覆盖所有数组元素
sumElements<<<grid, block>>>(device_array, n, device_sum); // 启动核函数
cudaMemcpy(&host_sum, device_sum, sizeof(int), cudaMemcpyDeviceToHost); // 获取结果
cudaFree(device_array);
cudaFree(device_sum);
std::cout << "The sum is: " << host_sum << std::endl;
return 0;
}
2.10.2归约的向量求和
void __global__ reduce_global(real* d_x, real* d_y)
{
const int tid = threadIdx.x;
//定义指针X，右边表示 d_x 数组第 blockDim.x * blockIdx.x个元素的地址
//该情况下x 在不同线程块，指向全局内存不同的地址---》使用不同的线程块对dx数组不同部分分别进行处理
real* x = d_x + blockDim.x * blockIdx.x;
//blockDim.x >> 1 等价于 blockDim.x /2，核函数中，位操作比 对应的整数操作高效
for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1)
{
if (tid < offset)
{
x[tid] += x[tid + offset];
}
//同步语句，作用：同一个线程块内的线程按照代码先后执行指令（块内同步，块外不用同步）
__syncthreads();
}
if (tid == 0)
{
d_y[blockIdx.x] = x[0];
}
}
解释：
当然可以。为了简化例子，我们假设blockDim.x为4（即每个线程块中有4个线程），并且我们有一个包含16个元素的d_x数组。这样，我们需要4个线程块来处理整个数组。
假设d_x的初始值为：
d_x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
线程块和线程的分配如下：
线程块0（线程0-3）：处理d_x[0]到d_x[3]
线程块1（线程0-3）：处理d_x[4]到d_x[7]
线程块2（线程0-3）：处理d_x[8]到d_x[11]
线程块3（线程0-3）：处理d_x[12]到d_x[15]
现在，我们以线程块0为例，解释代码的执行过程：
初始化：
tid（线程ID）对于线程块0中的线程来说，分别是0、1、2、3。
x指针对于线程块0将指向d_x的起始地址。
归约操作：
初始offset为blockDim.x >> 1，即4 >> 1 = 2。
在第一轮归约中，tid为0和1的线程会执行累加操作：
线程0（tid == 0）将d_x[0]（值为1）和d_x[2]（值为3）相加，结果4存储在d_x[0]。
线程1（tid == 1）将d_x[1]（值为2）和d_x[3]（值为4）相加，结果6存储在d_x[1]。
执行__syncthreads()确保所有线程完成操作。
第二轮归约，offset减半为1。
线程0将d_x[0]（现在的值为4）和d_x[1]（现在的值为6）相加，结果10存储在d_x[0]。
再次执行__syncthreads()。
存储结果：
最后，线程0将d_x[0]的值（现在的值为10）存储到d_y[blockIdx.x]，即d_y[0]。
同样的过程会在其他线程块中并行进行。最终，d_y数组将包含每个线程块归约后的结果：
d_y = [10, 26, 42, 58]
这些值是d_x数组中每个线程块所处理部分的元素和：
线程块0：1 + 2 + 3 + 4 = 10
线程块1：5 + 6 + 7 + 8 = 26
线程块2：9 + 10 + 11 + 12 = 42
线程块3：13 + 14 + 15 + 16 = 58
3、CPU矩阵乘向量 +openmp
void matmul(float* xout, float* x, float* w, int n, int d) {
// W (d,n) @ x (n,) -> xout (d,)
// by far the most amount of time is spent inside this little function
int i;
#pragma omp parallel for private(i)
for (i = 0; i < d; i++) {
float val = 0.0f;
for (int j = 0; j < n; j++) {
val += w[i * n + j] * x[j]; //w行不变，列j变，乘对应的x中第j个元素
}
xout[i] = val;
}
}
1） openmp
OpenMP是作为共享存储标准而问世的。它是为在多处理机上编写并行程序而设计的一个应用编程接口。它包括一套编译指导语句和一个用来支持它的函数库。
当今双核、四核的 CPU 当道，而六核的CPU也已经面世多时，所以在多处理机上编写、运行并行程序会变得相当普遍。
对于一般单线程（single thread）的程序，多核心的处理器并没有办法提升它的处理效能；不过对于多线程（multi thread）的程序，就可以通过不同的核心同时计算，来达到加速的目的了！简单的例子，以单线程的程序来说，一件事做一次要十秒的话，要做十次，都丢给同一颗核心做的话，自然就是 10 秒 * 10 次，也就是 100 秒了；但是以多线程的程序来说，它可以把这一件事，分给两颗核心各自做，每颗核心各做 5 次，所以所需要的时间就只需要 50 秒！
4、CPU矩阵乘
输入矩阵 A、B 分别由指针 h_a 和 h_b 指向的一维数组表示，输出矩阵 C 由指针 h_result 指向的一维数组表示。三个矩阵的维度分别为：A（m 行 × n 列）、B（n 行 × k 列）和 C（m 行 × k 列）。
void matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) {
for (int i = 0; i < m; ++i) //遍历输出矩阵 C 的每一行，用变量 i 作为行索引，范围从 0 到 m-1
{
for (int j = 0; j < k; ++j) //遍历输出矩阵 C 的每一列，用变量 j 作为列索引，范围从 0 到 k-1。
{
int tmp = 0.0;
//算当前 C[i][j] 元素的乘积累加项。用变量 h 作为临时索引，遍历矩阵 A 的当前行（对应于 i 行）和矩阵 B 的当前列（对应于 j
// 列）相交的所有元素。范围从 0 到 n-1。
//计算 A[i][h] 与 B[h][j] 的乘积，并累加到临时变量 tmp 上。
// C的第[i][j]元素就是A的第A[i][h] 与 B[h][j] 的乘积，i是外循环，j是中间循环，h是内循环，作为A的第i行和B的第j列的临时索引。
for (int h = 0; h < n; ++h)
{
tmp += h_a[i * n + h] * h_b[h * k + j];
}
//将累加结果 tmp 赋给输出矩阵 C 的对应元素 C[i][j]，即 h_result[i * k + j]
h_result[i * k + j] = tmp;
}
}
}
// C的第[i][j]元素就是A的第A[i][h] 与 B[h][j] 的乘积，i是外循环，j是中间循环，h是内循环，作为A的第i行和B的第j列的临时索引。
在CPU上，矩阵乘法通常使用三层嵌套循环来进行逐元素计算，如您提供的matrix_mult函数所示。然而，在GPU编程中，尤其是CUDA编程，我们尝试最大化并行性以利用GPU的众多核心同时进行计算。__global__ void matrixMul 函数正是采用了这种方式，它利用CUDA线程模型，将矩阵乘法的问题划分给许多并行运行的线程，每个线程负责计算输出矩阵的一个元素。
虽然在概念上，GPU版本也包含了类似“for”循环的思想（即，每个线程内在进行一次隐含的循环来计算矩阵乘积），但GPU编程模型下的并行化方法显著减少了实际的串行循环层数，并让每个线程独立地完成其任务，从而提高了计算速度。
在GPU版的matrixMul函数中，ix和iy是线程在矩阵C中的全局索引，而内部并没有显式的第三层循环去遍历矩阵A和B的列。这是因为每个GPU线程负责计算矩阵C中的单个元素，通过索引ix和iy直接定位到相应的乘法操作，然后在一个隐含的循环中（即 for(int i=0; i<k; i++)）计算这个元素的值。由于所有的线程并行执行这一操作，所以整个矩阵乘法的计算时间得以显著缩短。
简而言之，GPU编程之所以能够减少对传统for循环的依赖，是因为它可以利用大量的并行线程同时处理多个计算任务，而不是像CPU那样依次执行循环体内容。在这种情况下，尽管看起来只有一层明显的循环，但实际上，成千上万的线程各自在其内部执行着类似的计算流程，从而整体上加速了矩阵乘法的执行速度。
5、 自注意力机制（单头注意力）
https://zhuanlan.zhihu.com/p/659332084
https://zhuanlan.zhihu.com/p/366592542
为什么要缩放点积：对于较大的值，点积的大小变大，反向传播会将softmax函数推到具有极小梯度的区域。为了抵消这个影响，将点积进行缩放。
当输入向量的某些元素特别大或者特别小时，Softmax函数的导数可能会变得非常小，这就是所谓的“梯度消失”问题。因为是平的在0处和1处
https://www.bilibili.com/video/BV18B4y1k7d3/?spm_id_from=333.337.search-card.all.click&vd_source=8ab58b351436e118646830a9dac56e72
import torch
import math
import torch.nn as nn
class ScaledDotProductAttention(nn.Module):
def __init__(self):
super(ScaledDotProductAttention,self).__init__()
self.softmax = nn.Softmax(dim=-1)
def forward(self, q, k, v):
scores = torch.matmul(q, k.transpose(-2, -1)) / math.sqrt(k.shape[-1])
attn_weights = self.softmax(scores)
output = torch.matmul(attn_weights, v)
return output
6、 多头注意力
https://www.bilibili.com/video/BV1yY4y1i7v3?t=383.1&p=14
x1wq=q1, x2wq=q2
将x和w叠加起来，放到一个矩阵中去。得到一个qkv需要一个wq，wk，wv，那么得到多个，需要多个wq，wk，wv，如下面多头
由于w0q，w1q等维度是一样，因此在代码里面可以将它们叠加起来
八个头的情况
1代表1个句子，seq_len是5个字，4代表每个字用4个元素的向量表示。
import torch
import torch.nn as nn
import torch.nn.functional as F
import math
class MultiHeadAttention(nn.Module):
def __init__(self, embed_dim, num_heads,max_seq_len):
super(MultiHeadAttention,self).__init__()
self.num_heads = num_heads
self.embed_dim = embed_dim
self.head_dim = embed_dim // num_heads
self.wq = nn.Linear(embed_dim, self.head_dim * self.num_heads)
self.wk = nn.Linear(embed_dim, self.head_dim * self.num_heads)
self.wv = nn.Linear(embed_dim, self.head_dim * self.num_heads)
self.wo = nn.Linear(self.head_dim * self.num_heads, embed_dim)
self.register_buffer(
# `torch.tril`是方法，作用在已经有的tensor上
"mask", torch.tril(torch.ones(max_seq_len, max_seq_len))
)
def forward(self,x):
bsz, seqlen, _ = x.shape
# [bsz,seqlen,self.head_dim * self.num_heads]
xq, xk, xv = self.wq(x), self.wk(x), self.wv(x)
# Separate different heads，将self.head_dim * self.num_heads拆开
# [bsz, seqlen, self.num_heads, self.head_dim]
# [bsz, self.num_heads,m]
# transpose（1,2）都是分头进行计算attention，每个头单独计算，所以要交换
xq = xq.view(bsz, seqlen, self.num_heads, self.head_dim).transpose(1,2)
xk = xk.view(bsz, seqlen, self.num_heads, self.head_dim).transpose(1,2)
xv = xv.view(bsz, seqlen, self.num_heads, self.head_dim).transpose(1,2)
# qk相乘是最后两位相乘
scores = torch.matmul(xq, xk.transpose(2, 3))/math.sqrt(self.head_dim)
# mask
mask = self.mask[:seq_len, :seq_len]
scores = scores.masked_fill(mask == 0, float('-inf'))
scores = F.softmax(scores, dim=-1)
output = torch.matmul(scores, xv)
# 拼接起来，还原回去
output = output.transpose(1,2).contiguous()
output = output.view(bsz, seq, self.embed_dim)
output = self.wo(output)
return output
参考：https://zhuanlan.zhihu.com/p/675271011
https://zhuanlan.zhihu.com/p/669330242
https://zhuanlan.zhihu.com/p/365386753
7、nms
import numpy as np
bboxes = np.array([[100, 100, 210, 210, 0.72],
[250, 250, 420, 420, 0.8],
[220, 220, 320, 330, 0.92],
[100, 100, 210, 210, 0.72],
[230, 240, 325, 330, 0.81],
[220, 230, 315, 340, 0.9]])
def nms(iou_thresh=0.5, conf_threash = 0.5):
x1, y1, x2, y2, confidence = bboxes[:, 0], bboxes[:, 1], bboxes[:, 2], bboxes[:, 3], bboxes[:, 4]
# 计算面积
area = (x2 - x1) * (y2 - y1)
indices = confidence.argsort()[::-1]
keep = []
while indices.size > 0:
idx_self, idx_other = indices[0], indices[1:]
# 如果置信度小于阈值的话，那么后面的bbox就都不符合要求了，直接退出就行了
if confidence[idx_self] < conf_threash:
break
keep.append(idx_self)
# 计算交集
xx1, yy1 = np.maximum(x1[idx_self], x1[idx_other]), np.maximum(y1[idx_self], y1[idx_other])
xx2, yy2 = np.minimum(x2[idx_self], x2[idx_other]), np.minimum(y2[idx_self], y2[idx_other])
w, h = np.maximum(0, xx2 - xx1), np.maximum(0, yy2 - yy1)
intersection = w * h
# 计算并集(两个面积和-交集)
union = area[idx_self] + area[idx_other] - intersection
iou = intersection / union
# 只保留iou小于等于阈值的元素
print(np.where(iou <= iou_thresh))
keep_idx = np.where(iou <= iou_thresh)[0]
indices = indices[keep_idx + 1]
return np.array(keep)
if __name__ == '__main__':
print(nms())
8、冒泡排序
void Swap(int* pa, int* pb)
{
int tmp = *pa;
*pa = *pb;
*pb = tmp;
}
//冒泡排序
void BubbleSort(int* a, int n)
{
for (int i = 0; i < n; i++)
{
//单趟
for (int j = 1; j < n - i; j++)
{
if (a[j - 1] > a[j])
{
Swap(&a[j], &a[j - 1]);
}
}
}
}
📝：冒泡排序是最简单排序之一，简单意味着容易理解同时也意味着效率低。
🍒：时间复杂度：O(N^2)
🍒：空间复杂度：O(1)
🍒：稳定性：稳定
9、快速排序
//hoare法
int PartSort(int* a, int left, int right)
{
int keyi = left; //选最左边为key
while (left < right)
{
//右边先走，找到比key小的停下
while (left<right && a[right]>a[keyi]) //防止数组越界，用left<right控制
{
right--;
}
//左边走，找到比key大的停下
while (left < right && a[left] < a[keyi])
{
left++;
}
//交换数据
Swap(&a[left], &a[right]);
}
//left和right相遇，交换值
Swap(&a[left], &a[keyi]);
return left;
}
//快速排序
void QuickSort(int* a, int begin, int end)
{
//子区间只有一个值或者不存在就是最小子问题
if (begin >= end)
{
return;
}
int keyi = PartSort(a, begin, end);
//分成[begin,keyi-1]和[keyi+1,end]两个左右区间进行递归
QuickSort(a, begin, keyi - 1);
QuickSort(a, keyi + 1, end);
}
📝：快速排序整体的综合性能和使用场景都是比较好的，所以才敢叫快速排序。
🍒：时间复杂度：O(N*logN)
🍒：空间复杂度：O(logN)
🍒：稳定性：不稳定
⭐时间复杂度分析：
最好情况：每次选的key都是中位数，此时key左边序列和后边序列相同。时间复杂度为O(N*logN)
最坏情况：每次选的key都是最大或者最小的数。时间复杂度为O(N^2)
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
