cuda矩阵乘法(简单理解)

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

矩阵规模

矩阵A: n行,l列
矩阵B:l行,m列
结果矩阵C=A✖️B: n行,m列

一、一维并行

1.一维线程并行(thread)

一个线程并行一个A1行✖️B1列

1. 需要线程个数:n✖️m个,一个线程并行一个A1行✖️B1列, 即一个线程对应结果矩阵的C的一个点,e.g. c(1,1)=xid0的计算结果,所以一共需要:

即可知一共需要n✖️m个线程并行,每个线程处理一个向量✖️向量
2. 总共的第几个线程:xid=block.x*blockDim.x+threadIdx.x;
3. 因为线程总数为n✖️m个,所以当前线程xid对应A的第几行=xid/n,
当前线程xid对应B的第几列=xid%n;

 const int idx = blockIdx.x*blockDim.x+threadIdx.x;
    const int Arow = idx / n;
    const int Bcolumn = idx % n;

    //计算矩阵乘法
    if (Arow < n && Bcolumn < m)
    {
        float t = 0;

        for (i = 0; i < l; i++)
        {
            t += A[Arow * l + i] * B[i * m + Bcolumn];
        }
        C[Arow * m + Bcolumn] = t;
    }

2.一维块线程并行(block)

一个block代表矩阵A1行,块中一个线程代表B1列,B矩阵的列已每块线程数跳步.

1. 一共需要n个block,如果可以给的最大块数小于A的行数则已gridDim.x为block的跳步.
   线程数为j个,如果j小于B的列数则以blockDim.x为跳步

    
    //计算矩阵乘法
    for (int Arow=blockIdx.x;Arow<n;Arow+=gridDim.x)
    {
    for (int Bcolumn=threadIdx.x; Bcolumn < m;Bcolumn+=blockDim.x)
    {
        float num = 0.0;

        for (int i = 0; i < l; i++)
        {
            num += A[Arow * l + i] * B[i * m + Bcolumn];
        }
        C[Arow * m + Bcolumn] = num;
    }
    }

3.一维共享A一行程并行(shared)

上个方法中A的每行一个block,把A的每行给共享存储.

1. 共享存储速的关键是数据重复利用,因为从全局读一次数据,然后重复多用几次最合适,e.g.矩阵乘法中,A的每一行就成了B的m个列,即A1的每行都用重复用了m次.可将A的每行放进共享存储.

__shared__ shA[l];
int Arow=blockIdx.x;
for(int Acolumn=threadIdx.x; Acolumn < l;Bcolumn+=blockDim.x)
{
    shA[Acolumn]=A[Arow*l+Acolumn];//把A的每行给shared
}
__syncthreads();

for(int Bcolumn=threadIdx.x; Bcolumn<m; Bcolumn+=blockDim.x)//A的每行乘B的每列跳步为没块线程数
{
    float num = 0.0;
    for (int i = 0; i < l; i++)//A1行(已在shared中)乘B1列
    {
        num += shA[i] * B[i * m + Bcolumn];
    }
     C[Arow * m + Bcolumn] = num;
}    
   

二、二维并行

1.共享存储二维分块

1. 问题与改进
   当A的列数>shared memory最大容量时->采用分块存在shared中;
   一维时只存A矩阵->二维将B矩阵也分块存入;
2. 解决办法:设置两个二维的share memory : shA[][Tl_size]放A的shared子矩阵,shB[Tl_size][]放B的shared子矩阵,A横着更新shA子矩阵,B竖着更新shB子矩阵

template <int BLOCK_SIZE> __global__ void
matrixMulCUDA(float *C, float *A, float *B, int m, int n, int l)
{
 // 累加，得到行 * 列的值
 float numsub = 0.0;

 // 循环次数等于widthA / 16，把长向量点积运算转化为两个短向量点积后的和
 //因为B的宽条是竖着的,每宽条中宽是往下移动的,往下移动即代表动一行要增长m个元素,动一块即移动块宽Bsize✖️m个元素,所以步长为每次BLOCK_SIZE * m; 因为A每宽条是横着的,所以宽条中的块是横着移动,当前行每往右移动一块,增长块的块度Bsize的元素,即每宽条每移动一块步长就为Bsize
 for (int i = 0; i <l; i += BLOCK_SIZE)
 {
 // 定义A的共享子矩阵变量,因为__shared__声明，所以同一个block中的所有threads都可见，
 //每个thread填充一个元素，并计算一个行列乘积，减小带宽使用
 __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];

 // 定义A的共享子矩阵变量
 __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

 // 每个block包含16 * 16 个线程，所以每个线程负责一个矩阵元素的拷贝（注意同步）
 //A的shared块更新方式:行:1.块横着一块一块更新一宽条:l * BLOCK_SIZE * blockIdx.y为前面有的宽条:blockIdx.y表示第几个宽条,宽条的宽度为Bsize,则BLOCK_SIZE * blockIdx.y为前面一共有几行,每行有l列(BLOCK_SIZE * blockIdx.y)*l则为一共前面有多少元素. 2.hreadIdx.y 表示当前宽条的第几行,l * hreadIdx.y表示本块中本行在本块前的前面行一共的元素.所以一共第  (BLOCK_SIZE * blockIdx.y)+hreadIdx.y 行. 列:3.threadIdx.x表示本块第几个元素,本行一共第几个元素为threadIdx.x+i,所以l * BLOCK_SIZE * blockIdx.y + l * hreadIdx.y + threadIdx.x表示前面宽行中元素个数➕本宽条中前几行元素个数➕本行第几个元素
 //第几行:BLOCK_SIZE * blockIdx.y + threadIdx.y。 第几列:threadIdx.x+i
 As[ty][tx] = A[l * (BLOCK_SIZE * blockIdx.y + threadIdx.y) + (threadIdx.x+i)];
 //第几行:threadIdx.y + i.      第几列:BLOCK_SIZE * blockIdx.x + threadIdx.x
 Bs[ty][tx] = B[ m * (threadIdx.y + i) + (BLOCK_SIZE * blockIdx.x + threadIdx.x)];

 // Synchronize to make sure the matrices are loaded
 __syncthreads();

 // 每个线程计算 子矩阵的行列乘积，大循环外边还有累加，累加的是不同子矩阵点积和
 for (int k = 0; k < BLOCK_SIZE; ++k)
 {
 numsub += As[ty][k] * Bs[k][tx];
 }

 // 再次同步
 __syncthreads();
 }

if(BLOCK_SIZE * blockIdx.y + threadIdx.y<n && BLOCK_SIZE * blockIdx.x + threadIdx.x<m)
//C矩阵为n行m列的矩阵:C的行id=A的行id;C的列id=B的列id
//行:之前本行所在块上面块的一共行: blockIdx.y*Bsize,当前块第几行threadIdx.y,一共第几行:blockIdx.y*Bsize+threadIdx.y
//列:本列所在块之前块所有列:blockIdx.x*Bsize,当前块第几列:threadIdx.x;一共第几列:blockIdx.x*Bsize+threadIdx.x
//C中一共第几个元素:C行id✖️C列数+C列id=(blockIdx.y*Bsize+threadIdx.y)*m+blockIdx.x*Bsize+threadIdx.x

 C[(BLOCK_SIZE * blockIdx.y + threadIdx.y)*m + BLOCK_SIZE * blockIdx.x + threadIdx.x] = numsub;
}