目录
引言
在CUDA编程的初学过程中矩阵乘法是难点之一，本文就来一步一步进阶的讲解一下矩阵的乘法，全部代码放在了文章末尾，叠个小甲，算法成千上万，这里只展示我个人的理解，若有错误欢迎批评指正。
本文主要围绕以下三个任务展开，以下三个任务都分别使用全局内存以及共享内存实现。
1、实现两个4*4方阵的乘法
2、实现 [4,6]的矩阵以及[ 6,4 ]的矩阵相乘
3、实现两个任意维度的矩阵相乘
我们先介绍一下矩阵数据的读取方式，我们会在CPU端定义一个二维数组，然后使用cudaMalloc以及cudaMemcpy这两个函数会将其传入设备端，矩阵在设备中可以看作在存储器中按照行主序的方式进行线性存储的，所以在设备上索引A[i][j] 的值时则要用 A[i*4+j] (假设一行有4个数)
任务一：方阵乘法
全局内存
在这三个任务中，任务一算是矩阵乘法的一个特例，实现起来也较为简单。
先使用全局内存实现，先说一下计算资源分配的规则，因为两个4*4的矩阵AB相乘的结果C矩阵维度也是4*4，所以我们就分配1个block以及4*4的thread，一个线程就负责结果矩阵C对应位置的值的计算。如下图所示，线程一就计算A的第一行与B的第一列的内积，线程二就计算A的第一行与B的第二列的内积...
核函数代码如下。
// 分配了 4*4 的线程
__global__ void mission1_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y; // 第几行
int col = threadIdx.x; // 第几列
/*
若是分配的为2*2的block以及2*2的thread，则
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
其余不变
*/
for(int k=0; k<4; k++)
{
// d_a第row行所有元素与d_b第col列所有元素作点乘
d_result[row*4+col] += d_a[row*4+k] * d_b[k*4+col];
}
}共享内存
因为矩阵A，B的内容我们需要反复使用，为了加快运行速度，我们把他们放进共享内存中，我们同样分配4*4的thread以及一个block。注意！！！在使用共享内存时计算资源分配的规则与使用全局内存不同（虽然这里也是分配4*4的thread），这里先插个眼，在介绍任务二时再解释清楚。
下面说一下核函数的内容，我们定义两个共享的4*4二维数组（当然如果你觉得二维数组速度慢了你也可以定义长为16的一维数组,当然后续的程序可能得改变一点）目的是为了将A，B的数据全部放进去。因为这里定义的是二维数组，所以数据的读取方式为A[i][j]。
当我们赋值完共享内存后，计算方式与之前就大差不差了。核函数代码如下
// 分配4*4的线程
__global__ void mission1_shared(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
__shared__ int sh_a[4][4]; // 定义共享内存
__shared__ int sh_b[4][4];
for(int i=0; i<4; i++)
{
for(int j=0; j<4; j++)
{
sh_a[i][j] = d_a[i*4+j]; // 将数据放进共享内存中
sh_b[i][j] = d_b[i*4+j];
__syncthreads();
}
}
for(int i=0; i<4; i++)
{
// sh_a第row行所有元素与sh_b第col列所有元素作点乘
d_result[row*4+col] += sh_a[row][i] * sh_b[i][col];
}
}好的我们来增加一些难度，假设我们要运算的矩阵非常大，但是我的显卡最大支持一个block [1024,1024,64]的线程（你也可以自己查询最大的线程维度），如果一个block塞不下怎么办？这个时候就要涉及到矩阵的分块。
沿用之前的例子，我们将矩阵进行切分，这里切分为四块2*2大小的子矩阵，相应地我们分配的计算资源为2*2的block以及2*2的thread，在说明计算思路之前我们先要知道分块矩阵的计算规则，如下图所示。
那么我们核函数的计算思路也大差不差，因为我们分配了2*2的block，每一个block就代表了一个子矩阵。所以block1就负责计算A1*B1+A2*B3，block2就计算A1*B2+A2*B4，block3.....，block4......以此类推
我们通过两次for循环来计算加号（+）前和加号后的值。比如：第一次循环就计算block1的A1*B1，block2的A1*B2.....第二次循环就计算block1的A2*B3，block2的A2*B4......计算完之后累加到对应位置
而其中A1*B1同样是一个2*2的矩阵之间的计算，所以也同样需要进行两次for循环，那么2*2的矩阵之间的计算方式就可以参考我们之前计算4*4矩阵相乘的思路了。（超级简单，无非就是4变成了2）
好的，核函数代码如下：
// 分配 2*2->block 2*2->thread
__global__ void mission1_shared1(int* d_a, int* d_b, int* d_result)
{
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
__shared__ int sh_a[2][2], sh_b[2][2];
for(int i=0; i<2; i++)
{
// 先计算四个block中的加号前的内容（i=0时），然后计算加号后的内容(i=1时)
sh_a[threadIdx.y][threadIdx.x] = d_a[row*4 + i*2 + threadIdx.x];
sh_b[threadIdx.y][threadIdx.x] = d_b[(threadIdx.y + i*2)*4 + col];
__syncthreads();
for(int j=0; j<2; j++)
{
// 子矩阵是一个2*2矩阵之间的计算
d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}任务二：矩阵*其转置
在了解了方阵之间的乘法后我们再次增加一些难度，这次的任务为矩阵乘以它的转置。
以A*B ＝> [4,6] * [6,4]为例。
全局内存
老样子，我们先使用全局内存来实现它，有关全局内存下计算资源分配的规则上文已经说明，我们分配4*4的thread，然后循环六次来计算每个位置的值。话不多嗦代码如下：
// 分配线程4*4
__global__ void mission2_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
for(int j=0; j<6; j++)
{
// d_a第row行所有元素与d_b第col列所有元素作点乘
d_result[row*4 + col] += d_a[row*6 + j] * d_b[j*4 + col];
}
}共享内存
这里共享内存我们直接考虑矩阵分块的情况，这里我们将其进行四分块，分配2*2的block以及3*3的thread，对于矩阵A每小块2*3，对于矩阵B每小块3*2，下面我会解释在使用共享内存下计算资源如何分配的问题。
我们将矩阵分块后就可以将[4,6] * [6,4]问题看作是一个2*2矩阵之间的乘法，因为共享内存是对于一个块所共享的，所以C矩阵的block1就负责计算A1*B1+A2*B3, block2 负责A1*B2+A2*B4，与之前一致。
我们会将子矩阵放入共享内存中然后逐步计算加号前后的内容，因为A，B每一个小块的维度数不一样A为2*3，B为3*2，若我们分配2*3的thread那么B子矩阵最后一行没有线程负责计算它，若分配3*2的thread则A子矩阵的最后一列没有线程负责计算它,所以我们就干脆分配3*3的thread虽然会存在多余的线程，但是所有的数据都会得到计算。
那么在共享内存下计算资源分配的问题就迎刃而解了，目的就是让矩阵的每一块都有线程负责计算，也就是说线程的尺寸是一个方形，长度为分块后矩阵的最大边长（不是一定要分块）。举个例子，[16,24] * [24,16]，不分块时为 1个block，24*24的thread，分4块时为2*2的block，12*12的thread，分16块时为4*4的block和6*6的thread。
了解了共享内存下计算资源的分配问题后我们就可以开始编写核函数了。函数的逻辑和任务一的大差不差，代码如下
// 分配内存 block->2*2 thread->3*3
// 版本一
__global__ void mission2_shared(int* d_a, int* d_b, int* d_result)
{
__shared__ int sh_a[3][3], sh_b[3][3];
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
for(int i=0; i<2; i++)
{
// 先计算四个block中的加号前的内容（i=0时），然后计算加号后的内容(i=1时)
sh_a[threadIdx.y][threadIdx.x] = d_a[row*6 + i*3 + threadIdx.x];
sh_b[threadIdx.y][threadIdx.x] = d_b[col + (threadIdx.y + i*3)*4];
__syncthreads();
for(int j=0; j<3; j++)
{
// 向量点乘
d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}
// block->2*2 thread->3*3
// 版本二运行速度可能稍慢，但便于理解
__global__ void mission2_shared1(int* d_a, int* d_b, int* d_result)
{
__shared__ int sh_a[2][3], sh_b[3][2];
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
for(int i=0; i<2; i++)
{
// 线程多余的部分就不管他了
if(threadIdx.y<2) sh_a[threadIdx.y][threadIdx.x] = d_a[row*6 + i*3 + threadIdx.x];
if(threadIdx.x<2) sh_b[threadIdx.y][threadIdx.x] = d_b[col + (threadIdx.y + i*3)*4];
__syncthreads();
for(int j=0; j<3; j++)
{
if(threadIdx.x<2 && threadIdx.y<2) d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}任务三：任意维度矩阵乘法
终于到了任务三，这也是矩阵计算的一般情况，前两个任务都可以说是任务三的特例，搞懂了这个那么以后在进行任何矩阵运算的时候就能迎刃而解了。
我们将以A*B＝>[2,6]*[6*3]为例，C＝>[2,3]。
全局内存
首先使用全局内存来实现它，还记得计算资源如何分配吗，C矩阵维度是2*3所以我们分配一个block以及3*2的thread（注意不是2*3哦）。核函数思路就和之前一样，直接附上代码。
// 任务三 [2,6]*[6,3]->[2,3] block->1*1 thread->3*2
__global__ void mission3_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
for(int i=0; i<6; i++)
{
d_result[row*3 + col] += d_a[row*6 + i] * d_b[i*3 + col];
}
}共享内存
这次我们不考虑矩阵分块的情况，读者也可以自己构思一下如何编写分块下的核函数。
计算资源的分配问题也解释清楚了，这里就直接说了，我们分配6*6的thread和一个block。我们把A，B存进两个共享内存后，就直接使用共享内存进行计算了，代码如下。
// 任务三 [2,6]*[6,3]->[2,3] block->1*1 thread->6*6
__global__ void mission3_shared(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
__shared__ int sh_a[6][6], sh_b[6][6];
sh_a[row][col] = d_a[row*6 + col];
sh_b[row][col] = d_b[row*3 + col];
__syncthreads();
for(int i=0; i<6; i++)
{
// 只在规定范围内（result矩阵的最大尺寸下）计算
if(row<2 && col<3) d_result[row*3 + col] += sh_a[row][i] * sh_b[i][col];
}
}全部代码
// 矩阵乘积
/*
任务一：实现两个4*4的方阵相乘
任务二：实现4*6与6*4的矩阵相乘
任务三：实现任意维度矩阵相乘
子任务一：使用全局内存
子任务二：使用共享内存
*/
#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime.h>
// 分配了 4*4 的线程
__global__ void mission1_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y; // 第几行
int col = threadIdx.x; // 第几列
/*
若是分配的为2*2的block以及2*2的thread，则
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
其余不变
*/
for(int k=0; k<4; k++)
{
// d_a第row行所有元素与d_b第col列所有元素作点乘
d_result[row*4+col] += d_a[row*4+k] * d_b[k*4+col];
}
}
// 分配4*4的线程
__global__ void mission1_shared(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
__shared__ int sh_a[4][4]; // 定义共享内存
__shared__ int sh_b[4][4];
for(int i=0; i<4; i++)
{
for(int j=0; j<4; j++)
{
sh_a[i][j] = d_a[i*4+j]; // 将数据放进共享内存中
sh_b[i][j] = d_b[i*4+j];
__syncthreads();
}
}
for(int i=0; i<4; i++)
{
// sh_a第row行所有元素与sh_b第col列所有元素作点乘
d_result[row*4+col] += sh_a[row][i] * sh_b[i][col];
}
}
// 分配 2*2->block 2*2->thread
__global__ void mission1_shared1(int* d_a, int* d_b, int* d_result)
{
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
__shared__ int sh_a[2][2], sh_b[2][2];
for(int i=0; i<2; i++)
{
// 先计算四个block中的加号前的内容（i=0时），然后计算加号后的内容(i=1时)
sh_a[threadIdx.y][threadIdx.x] = d_a[row*4 + i*2 + threadIdx.x];
sh_b[threadIdx.y][threadIdx.x] = d_b[(threadIdx.y + i*2)*4 + col];
__syncthreads();
for(int j=0; j<2; j++)
{
// 子矩阵是一个2*2矩阵之间的计算
d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}
// 分配线程4*4
__global__ void mission2_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
for(int j=0; j<6; j++)
{
// d_a第row行所有元素与d_b第col列所有元素作点乘
d_result[row*4 + col] += d_a[row*6 + j] * d_b[j*4 + col];
}
}
// 分配内存 block->2*2 thread->3*3
// 版本一
__global__ void mission2_shared(int* d_a, int* d_b, int* d_result)
{
__shared__ int sh_a[3][3], sh_b[3][3];
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
for(int i=0; i<2; i++)
{
// 先计算四个block中的加号前的内容（i=0时），然后计算加号后的内容(i=1时)
sh_a[threadIdx.y][threadIdx.x] = d_a[row*6 + i*3 + threadIdx.x];
sh_b[threadIdx.y][threadIdx.x] = d_b[col + (threadIdx.y + i*3)*4];
__syncthreads();
for(int j=0; j<3; j++)
{
// 向量点乘
d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}
// block->2*2 thread->3*3
// 版本二运行速度可能稍慢，但便于理解
__global__ void mission2_shared1(int* d_a, int* d_b, int* d_result)
{
__shared__ int sh_a[2][3], sh_b[3][2];
int row = blockIdx.y*2 + threadIdx.y;
int col = blockIdx.x*2 + threadIdx.x;
for(int i=0; i<2; i++)
{
// 线程多余的部分就不管他了
if(threadIdx.y<2) sh_a[threadIdx.y][threadIdx.x] = d_a[row*6 + i*3 + threadIdx.x];
if(threadIdx.x<2) sh_b[threadIdx.y][threadIdx.x] = d_b[col + (threadIdx.y + i*3)*4];
__syncthreads();
for(int j=0; j<3; j++)
{
if(threadIdx.x<2 && threadIdx.y<2) d_result[row*4 + col] += sh_a[threadIdx.y][j] * sh_b[j][threadIdx.x];
}
}
}
// 任务三 [2,6]*[6,3]->[2,3] block->1*1 thread->3*2
__global__ void mission3_global(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
for(int i=0; i<6; i++)
{
d_result[row*3 + col] += d_a[row*6 + i] * d_b[i*3 + col];
}
}
// 任务三 [2,6]*[6,3]->[2,3] block->1*1 thread->6*6
__global__ void mission3_shared(int* d_a, int* d_b, int* d_result)
{
int row = threadIdx.y;
int col = threadIdx.x;
__shared__ int sh_a[6][6], sh_b[6][6];
sh_a[row][col] = d_a[row*6 + col];
sh_b[row][col] = d_b[row*3 + col];
__syncthreads();
for(int i=0; i<6; i++)
{
// 只在规定范围内（result矩阵的最大尺寸下）计算
if(row<2 && col<3) d_result[row*3 + col] += sh_a[row][i] * sh_b[i][col];
}
}
int main()
{
/*
任务一
*/
printf("----------------\n");
printf("任务一\n");
int *d_a, *d_b, *d_result;
int h_a[4][4], h_b[4][4], h_result[4][4];
//初始化两个矩阵
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
h_a[i][j] = i;
h_b[i][j] = j;
}
}
cudaMalloc((void**)&d_a, sizeof(int)*4*4);
cudaMalloc((void**)&d_b, sizeof(int)*4*4);
cudaMalloc((void**)&d_result, sizeof(int)*4*4);
cudaMemcpy(d_a, h_a, sizeof(int)*4*4, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, sizeof(int)*4*4, cudaMemcpyHostToDevice);
// 开始运算
// mission1_global<<<1,dim3(4,4,1)>>>(d_a, d_b, d_result);
// mission1_shared<<<1,dim3(4,4,1)>>>(d_a, d_b, d_result);
mission1_shared1<<<dim3(2,2,1),dim3(2,3,1)>>>(d_a, d_b, d_result);
// 返回计算结果
cudaMemcpy(h_result, d_result, sizeof(int)*4*4, cudaMemcpyDeviceToHost);
for(int i=0; i<4; i++)
{
for(int k=0; k<4; k++)
{
printf("h_result[%d][%d]:%d ", i, k, h_result[i][k]);
}
printf("\n");
}
/*
任务二
*/
printf("----------------\n");
printf("任务二 \n");
// 定义需要用到的数据
int h_a1[4][6], h_b1[6][4], h_result1[4][4];
int *d_a1, *d_b1, *d_result1;
for(int i=0; i<4; i++)
{
for(int j=0; j<6; j++)
{
h_a1[i][j] = i;
}
}
for(int i=0; i<6; i++)
{
for(int j=0; j<4; j++)
{
h_b1[i][j] = j;
}
}
cudaMalloc((void**)&d_a1, sizeof(int)*4*6);
cudaMalloc((void**)&d_b1, sizeof(int)*6*4);
cudaMalloc((void**)&d_result1, sizeof(int)*4*4);
cudaMemcpy(d_a1, h_a1, sizeof(int)*4*6, cudaMemcpyHostToDevice);
cudaMemcpy(d_b1, h_b1, sizeof(int)*6*4, cudaMemcpyHostToDevice);
// 开始计算
// mission2_global<<<1,dim3(4,4,1)>>>(d_a1, d_b1, d_result1);
mission2_shared<<<dim3(2,2,1),dim3(3,3,1)>>>(d_a1, d_b1, d_result1);
// mission2_shared1<<<dim3(2,2,1),dim3(3,3,1)>>>(d_a1, d_b1, d_result1);
// 返回计算结果
cudaMemcpy(h_result1, d_result1, sizeof(int)*4*4, cudaMemcpyDeviceToHost);
for(int i=0; i<4; i++)
{
for(int j=0; j<4; j++)
{
printf("h_result1[%d][%d]:%d ", i, j, h_result1[i][j]);
}
printf("\n");
}
/*
任务三
*/
printf("----------------\n");
printf("任务三 \n");
// 定义要用到的数据
int *d_a2, *d_b2, *d_result2;
int h_a2[2][6], h_b2[6][3], h_result2[2][3];
for(int i=0; i<2; i++)
{
for(int j=0; j<6; j++)
{
// printf("h_a2[%d][%d]:%d ", i, j, i);
h_a2[i][j] = i;
}
// printf("\n");
}
for(int i=0; i<6; i++)
{
for(int j=0; j<3; j++)
{
// printf("h_b2[%d][%d]:%d ", i, j, j);
h_b2[i][j] = j;
}
// printf("\n");
}
cudaMalloc((void**)&d_a2, sizeof(int)*2*6);
cudaMalloc((void**)&d_b2, sizeof(int)*6*3);
cudaMalloc((void**)&d_result2, sizeof(int)*2*3);
cudaMemcpy(d_a2, h_a2, sizeof(int)*2*6, cudaMemcpyHostToDevice);
cudaMemcpy(d_b2, h_b2, sizeof(int)*6*3, cudaMemcpyHostToDevice);
// 开始计算
// mission3_global<<<1,dim3(3,2,1)>>>(d_a2, d_b2, d_result2);
mission3_shared<<<1,dim3(6,6,1)>>>(d_a2, d_b2, d_result2);
// 返回计算结果
cudaMemcpy(h_result2, d_result2, sizeof(int)*2*3, cudaMemcpyDeviceToHost);
for(int i=0; i<2; i++)
{
for(int j=0; j<3; j++)
{
printf("h_result2[%d][%d]:%d ", i, j, h_result2[i][j]);
}
printf("\n");
}
// 程序结束 释放内存
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_result);
cudaFree(d_a1);
cudaFree(d_b1);
cudaFree(d_result1);
cudaFree(d_a2);
cudaFree(d_b2);
cudaFree(d_result2);
return 0;
}
结果展示
1024
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
