矩阵规模

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

一、一维并行

1.一维线程并行(thread)

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

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

xid0=a1hang✖️b1lie	xid1=a1hang✖️b2lie	…	xid(m-1)=a1hang✖️bmlie
xid(m)=a2hang✖️b1lie	xid(m+1)=a2hang✖️b2lie	…	xid(2m-1)=a2hang✖️bmlie
.	…	…	…
xid((n-1)*m)=anhang✖️b1lie	xid((n-1)*m+1)=anhang✖️b2lie	…	xid((n-1)*(m-1))=anhang✖️bmlie

即可知一共需要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矩阵的列已每块线程数跳步.

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

	thread0	thread1	thread2	thread3	thread0+blockDim.x	hread1+blockDim.x	…	thread3+m/blockDim.x
block0	a1hang✖️b1lie	a1hang✖️b2lie	a1hang✖️b3lie	a1hang✖️b4lie	a1hang✖️b5lie	a1hang✖️b6lie	…	a1hang✖️bmlie
block1	a2hang✖️b1lie	a2hang✖️b2lie	a2hang✖️b3lie	a2hang✖️b4lie	a2hang✖️b5lie	a2hang✖️b6lie	…	a2hang✖️bmlie
.	…	…	…	…	…	…	…
block(n-1)	anhang✖️b1lie	anhang✖️b2lie	anhang✖️b3lie	anhang✖️b4lie	anhang✖️b5lie	anhang✖️b6lie	…	anhang✖️bmlie

//计算矩阵乘法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的每行给共享存储.

共享存储速的关键是数据重复利用,因为从全局读一次数据,然后重复多用几次最合适,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.共享存储二维分块

问题与改进
当A的列数>shared memory最大容量时->采用分块存在shared中;
一维时只存A矩阵->二维将B矩阵也分块存入;
解决办法:设置两个二维的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的元素,即每宽条每移动一块步长就为Bsizefor (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+iAs[ty][tx] = A[l * (BLOCK_SIZE * blockIdx.y + threadIdx.y) + (threadIdx.x+i)];//第几行:threadIdx.y + i.      第几列:BLOCK_SIZE * blockIdx.x + threadIdx.xBs[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.xC[(BLOCK_SIZE * blockIdx.y + threadIdx.y)*m + BLOCK_SIZE * blockIdx.x + threadIdx.x] = numsub;
}

cuda矩阵乘法(简单理解)相关推荐

MegEngine| CUDA 矩阵乘法终极优化
前言单精度矩阵乘法(SGEMM)几乎是每一位学习 CUDA 的同学绕不开的案例,这个经典的计算密集型案例可以很好地展示 GPU 编程中常用的优化技巧,而能否写出高效率的 SGEMM Kernel , ...
CUDA 矩阵乘法优化
这个完全是基础知识啊~~ 哪不对大佬们帮忙指出啊 CUDA 矩阵乘法优化手段详解 Naive 实现的分析:到底差在哪里? 笔者面试过不少具有 CUDA 编程经验的校招同学,当提问使用 CUDA 编 ...
CUDA: 矩阵乘法优化
矩阵乘法是有实用价值的程序,我们会使用浮点数. 虽然矩阵乘法有点老套,不过因为它相当简单,而且也可以用来介绍一些有关 CUDA 的有趣性质. 矩阵乘法为了单纯起见,我们这里以方形的矩阵为例子.基本上 ...
CUDA矩阵乘法优化
前言纸上的来终觉浅,绝知此事要躬行. naive写法一个矩阵的乘法简单如下:C=A*B, 一般用gemm(A,B,C,M,N,K)来表示,其中的m,n,k代表的位置如下,默认是k表示消失的纬度. ...
矩阵分块与矩阵乘法的理解
block matrix:分块矩阵: 线性组合(linear combination) ⇔ 矩阵乘法(matrix multiplication): 矩阵和向量相乘,自然容易立即: 矩阵和矩阵相乘,则 ...
雅克比矩阵的简单理解（手稿）
如有错误或者不足之处,欢迎大家留言指正!
【poj3734】矩阵乘法求解
[题意] 给N个方块排成一列.现在要用红.蓝.绿.黄四种颜色的油漆给这些方块染色.求染成红色方块和染成绿色方块的个数同时为偶数的染色方案的个数,输出对10007取余后的答案.(1<=n<= ...
CUDA加速计算矩阵乘法进阶玩法(共享内存)
CUDA加速计算矩阵乘法&进阶玩法~共享内存一.基础版矩阵乘法二.为什么可以利用共享内存加速矩阵乘法 1.CUDA内存读写速度比较 2.申请共享内存三.改进版矩阵乘法(利用共享内存) 一 ...
将矩阵转为一行_理解矩阵乘法
理解矩阵乘法考研需要考一门课<线性代数>,这门课其实是教矩阵. 刚学的时候,还蛮简单的,矩阵加法就是相同位置的数字加一下. 矩阵乘法也类似,矩阵乘以一个常数,就是所有位置都乘以这个数. ...

cuda矩阵乘法(简单理解)

CUDA矩阵乘法