CUDA 矩阵相乘

2024-06-17 16:18
文章标签 cuda 相乘 矩阵

本文主要是介绍CUDA 矩阵相乘，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

优化点：

1.block分成多个，可以不受图片大小的限制 （一个block内的线程数量有限）

2.每个block内使用shared momery 可以优化数据访问速度

const int TILE_WIDTH = 32;
__global__ void mulKernel(int *c,  uchar *a,  uchar *b, int Width)
{__shared__ uchar Mds[TILE_WIDTH][TILE_WIDTH];__shared__ uchar Nds[TILE_WIDTH][TILE_WIDTH];//blockDim.x = blockDim.y = TILE_WIDTHint bx = blockIdx.x;int by = blockIdx.y;int tx = threadIdx.x;int ty = threadIdx.y;int row = TILE_WIDTH * by + ty;int col = TILE_WIDTH * bx + tx;int pValue = 0;for (size_t i = 0; i < Width/TILE_WIDTH; ++i){Mds[ty][tx] = a[row*Width+(i*TILE_WIDTH + tx)];  //a global memory.Nds[ty][tx] = b[(i*TILE_WIDTH + ty)*Width + col]; //b global memory.__syncthreads();for (size_t j = 0; j < TILE_WIDTH; ++j){pValue += ((int)(Mds[ty][j])) * Nds[j][tx];__syncthreads();}}c[row*Width+col] = pValue;
}
cudaError_t mulWithCuda(int *&c, uchar *&a, uchar *&b, unsigned int size)
{uchar *dev_a = 0;uchar *dev_b = 0;int *dev_c = 0;cudaError_t cudaStatus;// Choose which GPU to run on, change this on a multi-GPU system.cudaStatus = cudaSetDevice(0);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");goto Error;}// Allocate GPU buffers for three vectors (two input, one output)    .cudaStatus = cudaMalloc((void**)&dev_c, size * sizeof(int));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!");goto Error;}cudaStatus = cudaMalloc((void**)&dev_a, size * sizeof(uchar));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!");goto Error;}cudaStatus = cudaMalloc((void**)&dev_b, size * sizeof(uchar));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!");goto Error;}// Copy input vectors from host memory to GPU buffers.cudaStatus = cudaMemcpy(dev_a, a, size * sizeof(uchar), cudaMemcpyHostToDevice);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!!");goto Error;}cudaStatus = cudaMemcpy(dev_b, b, size * sizeof(uchar), cudaMemcpyHostToDevice);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!!!");goto Error;}double t0 = (double)cvGetTickCount();for (size_t i = 0; i < 10; i++){//共有 16*16个blockdim3 grid(16, 16); // 512/tile_width//每个block 有 32*32=1024个线程dim3 block(TILE_WIDTH, TILE_WIDTH);// Launch a kernel on the GPU with one thread for each element.mulKernel << <grid, block >> >(dev_c, dev_a, dev_b, 512);}double t1 = (double)cvGetTickCount();cout << "GPU time is:" << (double)(t1 - t0) / ((double)cvGetTickFrequency() * 1000) <<  "ms"<< endl;// Check for any errors launching the kernelcudaStatus = cudaGetLastError();if (cudaStatus != cudaSuccess) {fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));goto Error;}// cudaDeviceSynchronize waits for the kernel to finish, and returns// any errors encountered during the launch./*cudaStatus = cudaDeviceSynchronize();if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching mulKernel!\n", cudaStatus);goto Error;}*/// Copy output vector from GPU buffer to host memory.cudaStatus = cudaMemcpy(c, dev_c, size * sizeof(int), cudaMemcpyDeviceToHost);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!!!!");goto Error;}Error:cudaFree(dev_c);cudaFree(dev_a);cudaFree(dev_b);return cudaStatus;
}int main()
{const int img_size = 512;uchar *data1 = nullptr;uchar *data2 = nullptr;int *data = nullptr;change_data(rimg1, rimg2, data1, data2, data, img_size);//其中 矩阵的宽和高分别为 512mulWithCuda(data, data1, data2, img_size *img_size);return 0;
}

 

1. CUDA学习步骤

  1. CPU实现 a*b = c 的矩阵乘法（矩阵尺寸是n*m的，n和m大于1000）
  2. 下载 https://developer.nvidia.com/cuda-downloads，安装好cuda
  3. 将cpu代码移植到cuda。将CPU值传入GPU，使用cuda计算，与cpu结果对比。
  4. 优化思路1：将矩阵分块进行计算
  5. 优化思路2：使用share memory进行优化
  6. 优化思路3：将数据绑定在texture上

2. CPU实现的矩阵乘法

废话不多说，直接上源码

/* CPUMatMultiply：CPU下矩阵乘法* a:第一个矩阵指针，表示a[M][N]* b:第二个矩阵指针，表示b[N][S]* result:结果矩阵，表示为result[M][S]*/
void CPUMatMultiply(const int * a,const int * b, int *result,const int M,const int N,const int S)
{for (int i = 0; i < M; i++){for (int j = 0; j < S; j++){int index = i * S + j;result[index] = 0;//循环计算每一个元素的结果for (int k = 0; k < N; k++){result[index] += a[i * N + k] * b[k * S + j];}}}
}

3. CUDA实现的矩阵乘法

在https://developer.nvidia.com/cuda-downloads中下载相应的CUDA安装程序，进行安装。

PS：有关CUDA的环境搭建，Hello_World工程的创建可以移步http://www.mamicode.com/info-detail-327339.html。

下面直接进入正题，矩阵乘法的移植。 
从CPU上直接移植矩阵乘法到GPU上是非常简单的，不需要for循环，直接通过CUDA线程的id号，即threadIdx.x和threadIdx.y即可操作相应的数据。

gpu矩阵乘法核函数-源代码：

/* gpuMatMultKernel：GPU下矩阵乘法核函数
*  a:第一个矩阵指针，表示a[M][N]
*  b:第二个矩阵指针，表示b[N][S]
*  result:结果矩阵，表示result[M][S]
*/
__global__ void gpuMatMultKernel(const int *a, const int *b, int *result, const int M, const int N, const int S)
{int threadId = (blockIdx.y * blockDim.y + threadIdx.y) * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;if (threadId < M * S){int row = threadId / S;int column = threadId % S;result[threadId] = 0;for (int i = 0; i < N; i++){result[threadId] += a[row * N + i] * b[i * S + column];}}
}

其中，blockIdx、blockDim、threadIdx、gridDim都是CUDA的内置变量。 
blockIdx、threadIdx：分别CUDA中线程块的ID、线程的ID。 
blockDim、gridDim：分别是CUDA中线程块的维度，线程网格的维度。 
if语句是为了判断是否是对矩阵中的数据进行操作，防止在当线程threadId超过了M*S时，使用result[threadId]出现越界行为。 
后续的操作则是矩阵乘法的简单实现。

4. 程序优化1：共享内存分块运算

CUDA C支持共享内存。关键字__shared__添加到声明中，这将使这个变量驻留在共享内存中。当声明变量为shared变量时，线程块中的每一个线程都共享这块内存，使得一个线程块中的多个线程能够在计算上进行通信和协作。而且，共享内存缓冲区驻留在物理GPU上，而不是驻留在GPU之外的系统内存中。因此，在访问共享内存时的延迟要远远低于访问普通缓冲区的延迟，提高了效率。

使用共享内存的核函数：

/* gpuMatMultWithSharedKernel：GPU下使用shared内存的矩阵乘法
*  a:第一个矩阵指针，表示a[height_A][width_A]
*  b:第二个矩阵指针，表示b[width_A][width_B]
*  result:结果矩阵，表示result[height_A][width_B]
*/
template<int BLOCK_SIZE>
__global__ void gpuMatMultWithSharedKernel(const int *a, const int *b, int *result, const int height_A, const int width_A, const int width_B)
{int block_x = blockIdx.x;int block_y = blockIdx.y;int thread_x = threadIdx.x;int thread_y = threadIdx.y;if ((thread_y + block_y * blockDim.y) * width_B + block_x * blockDim.x + thread_x >= height_A * width_B){return;}const int begin_a = block_y * blockDim.y * width_A;const int end_a = begin_a + width_A - 1;const int step_a = blockDim.x;const int begin_b = block_x * blockDim.x;const int step_b = blockDim.y * width_B;int result_temp = 0;for (int index_a = begin_a, int index_b = begin_b;index_a < end_a; index_a += step_a, index_b += step_b){__shared__ int SubMat_A[BLOCK_SIZE][BLOCK_SIZE];__shared__ int SubMat_B[BLOCK_SIZE][BLOCK_SIZE];SubMat_A[thread_y][thread_x] = a[index_a + thread_y * width_A + thread_x];SubMat_B[thread_y][thread_x] = b[index_b + thread_y * width_B + thread_x];__syncthreads();for (int i = 0; i < BLOCK_SIZE; i++){result_temp += SubMat_A[thread_y][i] * SubMat_B[i][thread_x];}__syncthreads();}int begin_result = block_y * blockDim.y * width_B + begin_b;result[begin_result + thread_y * width_B + thread_x] = result_temp;
}

矩阵乘法的并行运算，每次计算矩阵的一块数据。利用共享内存的共享功能，每次将一块数据保存到共享内存中使得一个线程块同时调用数据进行计算当前块相对应得矩阵乘法结果值。

  • 代码 __shared__ int SubMat_A中的__shared__声明变量SubMat_A为共享内存中保存的变量。然后将数组中的数据提取到变量SubMat_A中保存在共享内存。
  • __syncthreads()对线程块中的线程进行同步，确保对__shared__进行下面的操作时上面的操作已经完成。
  • 两个for循环完成了当前线程块对应矩阵子块的乘法结果计算。

注意：CUDA架构将确保，除非线程块中的每个线程都执行了__syncthreads()，否则没有任何线程能执行__syncthreads()之后的指令。例如：

if (cacheIndex < i){cache[cacheIndex] += cache[i];__syncthreads();
}

此时，如果有线程中的cacheIndex>i则会出现一直等待的情况。

5. 程序优化2：纹理内存的运用

纹理内存，CUDA C程序中的另一种只读内存。纹理内存是专门为那些在内存访问模式中存在大量空间局部性的图形应用程序而设计的。在某个计算应用程序中，这意味着一个线程读取的位置可能与邻近线程读取的位置“非常接近”，如下图所示。 
这里写图片描述 
从数学的角度看，这四个地址并非连续的，在一般的CPU缓存模式中，这些地址将不会缓存。但由于GPU纹理缓存是专门为了加速这种访问模式而设计的，因此如果在这种情况中使用纹理内存而不是全局内存，那么将会获得性能提升。

使用纹理内存的核函数-源代码：

/* gpuMatMultWithTextureKernel：GPU下使用texture内存的矩阵乘法
*  result：结果矩阵，表示为result[M][S];
*  M：表示为矩阵A与矩阵result的行数
*  N：表示矩阵A的列数，矩阵B的行数
*  S：表示矩阵B和矩阵result的列数
*/
__global__ void gpuMatMultWithTextureKernel(int * result, const int M, const int N, const int S)
{int x = threadIdx.x + blockIdx.x * blockDim.x;int y = threadIdx.y + blockIdx.y * blockDim.y;int offset = x + y * blockDim.x * gridDim.x;if (offset < M * S){int a = 0, b = 0;int temp_result = 0;for (int i = 0; i < N; i++){a = tex1Dfetch(texA, y * N + i);b = tex1Dfetch(texB, i * S + x);temp_result += a * b;}result[offset] = temp_result;}
}

纹理内存的运用与普通内存运用时的算法大致相同，只不过数据是在核函数中调用tex1Dfetch从纹理中提取。

在使用纹理内存时，主要注意的是纹理内存的使用。 
首先，需要将输入的数据声明为texture类型的引用。 
注意，输入的数据是什么类型，相应的纹理也应该与之一致。并且纹理引用必须声明为文件作用域内的全局变量。

//这些变量将位于GPU上
texture<int> texA;
//二维纹理引用，增加了代表维数的参数2
texture<float, 2> texB;

在为这两个缓冲区分配了GPU内存后，需要通过cudaBindTexture将这些变量绑定到内存缓冲区。这相当于告诉CUDA运行时两件事：

  • 我们希望将指定的缓冲区作为纹理来使用。
  • 我们希望将纹理引用作为纹理的“名字”。
cudaBindTexture(NULL, texA, dev_a, desc, M * N * sizeof(int));
cudaBindTexture(NULL, texB, dev_b, desc, N * S * sizeof(int));

在绑定纹理时，CUDA运行时要求提供一个cudaChannelFormatDesc。此时，需要调用cudaCreateChannelDesc<int>()。

最后，通过cudaUnbindTexture()函数来取消纹理的绑定。

6.总结

矩阵乘法计算中，运行时间上有：共享内存 < 纹理内存 ≈ 普通GPU移植 < CPU运算。 
具体源代码https://github.com/ZYMing/CUDA_Samples/blob/master/MatMulti/MatMulti.cu。

 

计算性能对比：

512*512的矩阵

CPU time is:174.119 ms

GPU time is:0.0288702ms  （cu核函数）

CUDABlas cost time is:0.667336 ms

该时间未包含 cuda内存申请和拷贝的时间。

这篇关于CUDA 矩阵相乘的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/1069932。

相关文章

推荐算法之矩阵分解实例

推荐算法之矩阵分解实例

矩阵分解的数据利用的上篇文章的数据，协同过滤 用到的知识 python的surprise k折交叉验证 SVD SVDpp NMF 算法与结果可视化 # 可以使用上面提到的各种推荐系统算法from surprise import SVD,SVDpp,NMFfrom surprise import Datasetfrom surprise import print_perf
阅读更多...
江协科技51单片机学习- p16 矩阵键盘

江协科技51单片机学习- p16 矩阵键盘

🚀write in front🚀   🔎大家好，我是黄桃罐头，希望你看完之后，能对你有所帮助，不足请指正！共同学习交流 🎁欢迎各位→点赞👍 + 收藏⭐️ + 留言📝​  💬本系列哔哩哔哩江科大51单片机的视频为主以及自己的总结梳理📚  前言： 本文是根据哔哩哔哩网站上“江协科技51单片机”视频的学习笔记，在这里会记录下江协科技51单片机开发板的配套视频教程所作的实验和学习
阅读更多...
使用matlab的大坑，复数向量转置！！！！！变量区“转置变量“功能（共轭转置）、矩阵转置（默认也是共轭转置）、点转置

使用matlab的大坑，复数向量转置！！！！！变量区“转置变量“功能（共轭转置）、矩阵转置（默认也是共轭转置）、点转置

近期用verilog去做FFT相关的项目，需要用到matlab进行仿真然后和verilog出来的结果来做对比，然后计算误差。近期使用matlab犯了一个错误，极大的拖慢了项目进展，给我人都整emo了，因为怎么做仿真结果都不对，还好整体的代码都是我来写的，慢慢往下找就找到了问题的来源，全网没有看到多少人把这个愚蠢的错误写出来，我来引入一下。 代码错误的表现：复数向量的虚部被取反，正数变成负数，负数
阅读更多...
「学转录组入门生信」第二周来获取表达量矩阵

「学转录组入门生信」第二周来获取表达量矩阵

我们第二周目标有四个: 整理数据RNA-seq格式了解数据质控数据比对read定量 首先，我们得要知道我们在转录组分析过程中会遇到很多格式，建议先通过搜索查找了解这些格式是什么 fasta/fas/fagtf/gffbedsam/bamcsv/tsv/txt 接着，我们会在分析过程中时刻检查我们的数据质量，所以你要尝试回答下面这几个问题 数据质控要在哪个阶段做不同阶段要看什么标准质控有哪
阅读更多...
「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析

「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析

这一篇文章是回答知识星球中一位星友的提问，她的电脑内存有限，无法直接使用所有数据，只能分析部分数据。 数据来源: https://content.cruk.cam.ac.uk/jmlab/atlas_data.tar.gz 解压缩之后，得到下面数据 数据清单 其中raw_counts.mtx是以稀疏矩阵格式存放的表达量数据，文件为6.5G, 用普通的文本编辑器无法打开，
阅读更多...
Strassen矩阵乘法简要解析（第4章：分治策略）

Strassen矩阵乘法简要解析（第4章：分治策略）

Strassen矩阵乘法简要解析 Strassen矩阵乘法具体描述如下： 两个n×n 阶的矩阵A与B的乘积是另一个n×n 阶矩阵C，C可表示为假如每一个C(i, j) 都用此公式计算，则计算C所需要的操作次数为n3 m+n2 (n- 1) a，其中m表示一次乘法，a 表示一次加法或减法。 为了使讨论简便，假设n 是2的幂（也就是说， n是1，2，4，8，1 6，...）。 首先，假设
阅读更多...
[leetcode] 43. Multiply Strings(大数相乘)

[leetcode] 43. Multiply Strings(大数相乘)

Multiply Strings 描述 Given two non-negative integers num1 and num2 represented as strings, return the product of num1 and num2. Note: 1. The length of both num1 and num2 is < 110. 2. Both num1 an
阅读更多...
leetcode 二分查找·系统掌握 搜索二维矩阵

leetcode 二分查找·系统掌握 搜索二维矩阵

题目： 题解： 一个可行的思路是使用~01~泛型对每一行的最后一个元素进行查找找到第一个大于等于target的那一行，判断查找结果如果“失败”返回false否则继续在改行进行常规二分查找target的值根据查找结果返回即可。 bool searchMatrix(vector<vector<int>>& matrix, int target) {int l=0,r=matrix.size()-
阅读更多...
【LeetCode热题 100】螺旋矩阵

【LeetCode热题 100】螺旋矩阵

leetcode原地址：https://leetcode.cn/problems/spiral-matrix/description 描述 给你一个 m 行 n 列的矩阵 matrix ，请按照 顺时针螺旋顺序 ，返回矩阵中的所有元素。 示例 1： 输入：matrix = [[1,2,3],[4,5,6],[7,8,9]] 输出：[1,2,3,6,9,8,7,4,5] 示例 2：
阅读更多...
LaTeX中添加矩阵分块虚线并设置虚线疏密

LaTeX中添加矩阵分块虚线并设置虚线疏密

对于大型矩阵，有时需要添加分块虚线。 方法为使用arydshln宏包，然后在array环境中设置虚线。需要注意的是，使用矩阵环境需要搭配amsmath宏包使用，且需放在amsmath宏包之后。即导言区设置为 \usepackage{amsmath}\usepackage{arydshln} % 导入arydshln包 给出示例 \[\begin{bmatrix}\begin{array
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2