基于矩阵乘的CUDA编程优化过程

背景:网上很多关于矩阵乘的编程优化思路,本着看理论分析万遍,不如实际代码写一遍的想法,大概过一下优化思路。

矩阵乘的定义如下,约定矩阵的形状及存储方式为: A[M, K], B[K, N], C[M, N]。

C_{i,j}=\sum_{k=0}^{n}A_{ik}\times B_{kj}

CPU篇

朴素实现方法

        按照常规的思路,实现矩阵乘时如下的3层for循环。

#define OFFSET(row, col, ld) ((row) * (ld) + (col))
void cpuSgemm(float *a, float *b, float *c, const int M, const int N, const int K) 
{for (int m = 0; m < M; m++) {for (int n = 0; n < N; n++) {float psum = 0.0;for (int k = 0; k < K; k++) {psum += a[OFFSET(m, k, K)] * b[OFFSET(k, n, N)];}c[OFFSET(m, n, N)] = psum;}}
}

数据访存连续的优化

        矩阵B的存储默认为N方向连续,所以可以将上面的第2,3层循环互换顺序,这样B的取数就不会跨行了,而是连续取数,达到访问连续的效果。

void cpuSgemm_1(float *a, float *b, float *c, const int M, const int N, const int K) 
{for (int m = 0; m < M; m++) {for (int k = 0; k < K; k++) {for (int n = 0; n < N; n++){c[OFFSET(m, n, N)] += a[OFFSET(m, k, K)] * b[OFFSET(k, n, N)];}           }}
}

数据重排/数据复用的优化

        上面将M,N,K的for循环调整为M,K,N的循环顺序,导致我们K方向累加不能缓存了,增加了多次访问C矩阵的开销,所以我们不放先直接将B矩阵转置处理,然后再按照原始的M,N,K的for循环来处理。

void cpuSgemm_2(float *a, float *b, float *c, const int M, const int N, const int K) 
{float* b1=(float*) malloc(sizeof(float)*K*N);for(int i=0; i<K; i++){for (int j=0; j<N; j++){b1[OFFSET(j,i,K)]= b[OFFSET(i,j,N)];}}for (int m = 0; m < M; m++) {for (int n = 0; n < N; n++) {float psum = 0.0;for (int k = 0; k < K; k++) {psum += a[OFFSET(m, k, K)] * b1[OFFSET(n, k, K)];}c[OFFSET(m, n, N)] = psum;}}
}

性能表现

        如下是测试CPU环境下这几种方法的时间情况,其中M=N=512, K =256。可以发现经过优化后的代码在时间上是逐步减少的。

        CPU的优化思路还有其他的,比如循环展开,intrinsic函数,基于cache的矩阵切分等,注意本文并没有都实现出来。

cpuSgemm, Time measured: 416889 microseconds.
cpuSgemm_1, Time measured: 405259 microseconds.
cpuSgemm_2, Time measured: 238786 microseconds.

GPU篇

grid线程循环矩阵乘法

        输出矩阵C有M*N个点,每个点是K个数的乘积和,所以可以定义每个线程计算K个点的乘积和,即grid线程循环矩阵乘法。

__global__ void matrix_multiply_gpu_0(float*a, float*b, float*c, int M, int N, int K)
{int tidx =threadIdx.x;int bidx = blockIdx.x;int idx = bidx * blockDim.x +tidx;int row = idx/N;int col = idx%N;if(row<M && col < N){float tmp =0.0;for(int k=0; k<K; k++){tmp+=a[row*K+k] * b[k*N+col];}c[row*N+col] = tmp;}
}

block线程循环矩阵乘法

        grid内线程循环的矩阵乘法有如下缺憾:一个block内线程可能需要计算C矩阵不同行的矩阵元素,block内thread对相应的A矩阵访存不一致,导致无法广播和额外的访存开销,导致执行时间增加。

        针对这个问题,可以做如下改进:每个block计算C矩阵的一行,block内的thread以固定跳步步长blockDim.x的方法循环计算C矩阵的一行,每一行启动一个block,共计M个block。

__global__ void matrix_multiply_gpu_1(float*a, float*b, float*c, int M, int N, int K)
{int tidx =threadIdx.x;int bidx = blockIdx.x;float tmp;for(;bidx<M; bidx += gridDim.x){for(;tidx<N; tidx+=blockDim.x ){tmp=0.0;for(int k=0; k<K; k++){tmp+=a[bidx*K +k] * b[k*N+tidx];}c[bidx*N+tidx] = tmp;}              }
}

行共享存储矩阵乘法

        共享存储与L1 Cache同级,其访存延迟较全局存储小一个量级。用共享存储代替全局存储是GPU最重要的优化手段之一。采用共享存储优化的关键是数据复用,数据复用次数越多,共享存储优化可获得的收益也越高。

        在block循环乘法中,1个block内所有thread都会用到A矩阵的一行,此时与B矩阵每一列相乘,A矩阵中该行复用了N次。故可以考虑将A矩阵的一行读入shared memory,运算时候从shared memory读取相应的数据。

        注意代码中TILE_WIDTH>=K。

#define TILE_WIDTH 256
__global__ void matrix_multiply_gpu_2(float*a, float*b, float*c, int M, int N, const int K)
{__shared__ float data[TILE_WIDTH];int tid = threadIdx.x;int row = blockIdx.x;int i,j;for(i=tid; i<K; i+=blockDim.x){data[i]=a[row*K +i];}__syncthreads();float tmp;for(j=tid; j<N; j+=blockDim.x){tmp=0.0;for(int k=0; k<K; k++){tmp += data[k]*b[k*N+j];}c[row*N+j] = tmp;}
}

分块共享存储矩阵乘法

        根据上面共享存储的理解,我们很自然的想到把B矩阵也考虑数据复用,所以可以同时把A,B矩阵都分成棋盘似的小尺寸的数据块,从全局内存读取到共享内存,这样可以有效降低数据访问时间,充分复用矩阵乘的局部数据。

#define TILE_SIZE 32
__global__ void matrix_multiply_gpu_3(float*a, float*b, float*c, int M, int N, const int K)
{__shared__ float matA[TILE_SIZE][TILE_SIZE];__shared__ float matB[TILE_SIZE][TILE_SIZE];int bx = blockIdx.x;int by = blockIdx.y;int tx = threadIdx.x;int ty = threadIdx.y;int Col = bx * TILE_SIZE + tx;int Row = by * TILE_SIZE + ty;float Pervalue = 0.0;for(int i = 0;i < K / TILE_SIZE;i++)  {matA[ty][tx] = a[Row * K + (i * TILE_SIZE + tx)];matB[ty][tx] = b[Col + (i * TILE_SIZE + ty) * N];__syncthreads();for(int k = 0;k < TILE_SIZE;k++) Pervalue += matA[ty][k] * matB[k][tx];__syncthreads();}c[Row * N + Col] = Pervalue;}

性能表现

利用nvprof工具,统计各个核函数的执行时间如下,可以发现每一步优化思路都能直观的带来的性能提升。

完整代码:

GitHub - Briwisdom/study_CUDA_examples: some demos for study CUDA program.

#include <iostream>
#include <chrono>using namespace std;#define OFFSET(row, col, ld) ((row) * (ld) + (col))void initDate(float *arr,int Len, bool randFlag=true)
{if (randFlag){for (int i = 0; i < Len; i++) {arr[i] = rand()/1000000;}}else{float value =0.0;for (int i = 0; i < Len; i++) {arr[i] = value;}}  
}void compare_result(float *x, float *y, int n, char *name)
{int cnt=0;for (int i=0; i<n; i++){if (x[i]!=y[i]){cnt++;printf("x= %f, y= %f\n", x[i],y[i]);}}printf("%s, ", name);if(cnt ==0)printf("result matched.\n");elseprintf("something error! result not match number = %d int total number: %d .\n", cnt, n);}void cpuSgemm(float *a, float *b, float *c, const int M, const int N, const int K) 
{for (int m = 0; m < M; m++) {for (int n = 0; n < N; n++) {float psum = 0.0;for (int k = 0; k < K; k++) {psum += a[OFFSET(m, k, K)] * b[OFFSET(k, n, N)];}c[OFFSET(m, n, N)] = psum;}}
}void cpuSgemm_1(float *a, float *b, float *c, const int M, const int N, const int K) 
{for (int m = 0; m < M; m++) {for (int k = 0; k < K; k++) {for (int n = 0; n < N; n++){c[OFFSET(m, n, N)] += a[OFFSET(m, k, K)] * b[OFFSET(k, n, N)];}           }}
}void cpuSgemm_2(float *a, float *b, float *c, const int M, const int N, const int K) 
{float* b1=(float*) malloc(sizeof(float)*K*N);for(int i=0; i<K; i++){for (int j=0; j<N; j++){b1[OFFSET(j,i,K)]= b[OFFSET(i,j,N)];}}for (int m = 0; m < M; m++) {for (int n = 0; n < N; n++) {float psum = 0.0;for (int k = 0; k < K; k++) {psum += a[OFFSET(m, k, K)] * b1[OFFSET(n, k, K)];}c[OFFSET(m, n, N)] = psum;}}
}void operation(void (*func)(float*,float*, float*, int, int, int), float *a, float *b, float *c, const int M, const int N, const int K, int repeat, char* name)
{auto begin0 = std::chrono::high_resolution_clock::now();for(int i=0; i<repeat; i++){(*func)(a,b,c, M, N, K);}auto end0 = std::chrono::high_resolution_clock::now();auto elapsed0 = std::chrono::duration_cast<std::chrono::microseconds>(end0 - begin0);printf("%s, Time measured: %d microseconds.\n", name, int(elapsed0.count()/repeat));
}__global__ void matrix_multiply_gpu_0(float*a, float*b, float*c, int M, int N, int K)
{int tidx =threadIdx.x;int bidx = blockIdx.x;int idx = bidx * blockDim.x +tidx;int row = idx/N;int col = idx%N;if(row<M && col < N){float tmp =0.0;for(int k=0; k<K; k++){tmp+=a[row*K+k] * b[k*N+col];}c[row*N+col] = tmp;}
}__global__ void matrix_multiply_gpu_1(float*a, float*b, float*c, int M, int N, int K)
{int tidx =threadIdx.x;int bidx = blockIdx.x;float tmp;for(;bidx<M; bidx += gridDim.x){for(;tidx<N; tidx+=blockDim.x ){tmp=0.0;for(int k=0; k<K; k++){tmp+=a[bidx*K +k] * b[k*N+tidx];}c[bidx*N+tidx] = tmp;}              }
}#define TILE_WIDTH 256
__global__ void matrix_multiply_gpu_2(float*a, float*b, float*c, int M, int N, const int K)
{__shared__ float data[TILE_WIDTH];int tid = threadIdx.x;int row = blockIdx.x;int i,j;for(i=tid; i<K; i+=blockDim.x){data[i]=a[row*K +i];}__syncthreads();float tmp;for(j=tid; j<N; j+=blockDim.x){tmp=0.0;for(int k=0; k<K; k++){tmp += data[k]*b[k*N+j];}c[row*N+j] = tmp;}
}#define TILE_SIZE 32
__global__ void matrix_multiply_gpu_3(float*a, float*b, float*c, int M, int N, const int K)
{__shared__ float matA[TILE_SIZE][TILE_SIZE];__shared__ float matB[TILE_SIZE][TILE_SIZE];int bx = blockIdx.x;int by = blockIdx.y;int tx = threadIdx.x;int ty = threadIdx.y;int Col = bx * TILE_SIZE + tx;int Row = by * TILE_SIZE + ty;float Pervalue = 0.0;for(int i = 0;i < K / TILE_SIZE;i++)  {matA[ty][tx] = a[Row * K + (i * TILE_SIZE + tx)];matB[ty][tx] = b[Col + (i * TILE_SIZE + ty) * N];__syncthreads();for(int k = 0;k < TILE_SIZE;k++) Pervalue += matA[ty][k] * matB[k][tx];__syncthreads();}c[Row * N + Col] = Pervalue;}int main()
{int M=512;int N=512;int K=256;float *a = (float*) malloc(M*K * sizeof(float));float *b = (float*) malloc(N*K * sizeof(float));float *c = (float*) malloc(M*N * sizeof(float));float *c1 = (float*) malloc(M*N * sizeof(float));float *c2 = (float*) malloc(M*N * sizeof(float));float *c_gpu_0 = (float*) malloc(M*N * sizeof(float));float *c_gpu_1 = (float*) malloc(M*N * sizeof(float));float *c_gpu_2 = (float*) malloc(M*N * sizeof(float));float *c_gpu_3 = (float*) malloc(M*N * sizeof(float));initDate(a,M*K);initDate(b,N*K);initDate(c, M*N, false);initDate(c1, M*N, false);initDate(c2, M*N, false);initDate(c_gpu_0, M*N, false);initDate(c_gpu_1, M*N, false);initDate(c_gpu_2, M*N, false);initDate(c_gpu_3, M*N, false);//ensure result is right.cpuSgemm(a,b,c,M,N,K);cpuSgemm_1(a,b,c1,M,N,K);cpuSgemm_2(a,b,c2,M,N,K); compare_result(c, c1, M*N,"sgemm1");compare_result(c, c2,  M*N,"sgemm2");//test the prerformance.int repeat =10;operation(cpuSgemm,a,b,c,M,N,K,repeat,"cpuSgemm");operation(cpuSgemm_1,a,b,c1,M,N,K,repeat,"cpuSgemm_1");operation(cpuSgemm_2,a,b,c2,M,N,K,repeat,"cpuSgemm_2");float* d_a, *d_b, *d_c0, *d_c1, *d_c2, *d_c3;cudaMalloc((void**) &d_a, sizeof(float)*(M*K));cudaMalloc((void**) &d_b, sizeof(float)*(N*K));cudaMalloc((void**) &d_c0, sizeof(float)*(M*N));cudaMalloc((void**) &d_c1, sizeof(float)*(M*N));cudaMalloc((void**) &d_c2, sizeof(float)*(M*N));cudaMalloc((void**) &d_c3, sizeof(float)*(M*N));cudaMemcpy(d_a, a, sizeof(float)*M*K, cudaMemcpyHostToDevice);cudaMemcpy(d_b, b, sizeof(float)*N*K, cudaMemcpyHostToDevice);int threadnum=64;int blocks =(M*N+threadnum-1)/threadnum;cudaMemcpy(d_c0, c_gpu_0, sizeof(float)*M*N, cudaMemcpyHostToDevice);matrix_multiply_gpu_0<<<blocks, threadnum>>>(d_a, d_b, d_c0, M, N, K);cudaMemcpy(c_gpu_0, d_c0, sizeof(float)*M*N, cudaMemcpyDeviceToHost);compare_result(c, c_gpu_0,  M*N,"gpu_0");cudaFree(d_c0);cudaMemcpy(d_c1, c_gpu_1, sizeof(float)*M*N, cudaMemcpyHostToDevice);matrix_multiply_gpu_1<<<M, threadnum>>>(d_a, d_b, d_c1, M, N, K);cudaMemcpy(c_gpu_1, d_c1, sizeof(float)*M*N, cudaMemcpyDeviceToHost);compare_result(c, c_gpu_1,  M*N,"gpu_1");cudaFree(d_c1);cudaMemcpy(d_c2, c_gpu_2, sizeof(float)*M*N, cudaMemcpyHostToDevice);matrix_multiply_gpu_2<<<M, threadnum>>>(d_a, d_b, d_c2, M, N, K);cudaMemcpy(c_gpu_2, d_c2, sizeof(float)*M*N, cudaMemcpyDeviceToHost);compare_result(c, c_gpu_2,  M*N,"gpu_2");cudaFree(d_c2);threadnum=32;dim3 gridSize(M / threadnum,N / threadnum);dim3 blockSize(threadnum,threadnum);cudaMemcpy(d_c3, c_gpu_3, sizeof(float)*M*N, cudaMemcpyHostToDevice);matrix_multiply_gpu_3<<<gridSize, blockSize>>>(d_a, d_b, d_c3, M, N, K);cudaMemcpy(c_gpu_3, d_c3, sizeof(float)*M*N, cudaMemcpyDeviceToHost);compare_result(c, c_gpu_3,  M*N,"gpu_3");cudaFree(d_c3);free(a);free(b);free(c);free(c1);free(c2);free(c_gpu_0);free(c_gpu_1);free(c_gpu_2);free(c_gpu_3);cudaFree(d_a);cudaFree(d_b);}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/229615.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

刷算法-- leetcode 96. 不同的二叉搜索树

思路 观察树的组成&#xff0c;可以发现n3时的二叉搜索树可以由&#xff0c;头节点分别为1、2、3时的所有结果组成&#xff01;定义dp[i]为由i个节点组成的二叉搜索树的个数。确定递推公式&#xff0c;dp[i] 由1为头节点组成的二叉搜索树个数由2为头组成的个数…由i为头节点组…

使用 go-elasticsearch v8 基本请求

使用 go-elasticsearch 请求示例 你可以通过参考Go 官方文档找到简单的示例&#xff0c;所以我认为先看看这个是个好主意。 连接客户端有两种方式&#xff0c;如下图。 至于两者的特点&#xff0c;TypedClient有类型&#xff0c;更容易编写&#xff0c;但文档较少。另外&…

Linux文件的扩展属性 attr cap

文件属性 Linux文件属性分为常规属性与扩展属性&#xff0c;其中扩展属性有两种&#xff1a;attr与xattr. 一般常规的文件属性由stat API 读取&#xff0c;一般是三种权限&#xff0c;ower, group&#xff0c;时间等。 扩展属性attr 用户态API ioctl(fd, FS_IOC32_SETFLAGS…

QQ邮件发送(PHP的Laravel)

1. 开启 QQ 邮箱的 SMTP 支持 2.里面会一个类似于密码之类&#xff08;复制一下&#xff09; 3.然后再 .env文件里面配置一下 MAIL_DRIVERsmtp —— 使用支持 ESMTP 的 SMTP 服务器发送邮件&#xff1b; MAIL_HOSTsmtp.qq.com —— QQ 邮箱的 SMTP 服务器地址&#xff0c;必…

OCP NVME SSD规范解读-3.NVMe管理命令-part2

NVMe-AD-8&#xff1a;在某些情况下&#xff08;如Sanitize命令、Format NVM命令或TCG Revert方法后数据被清除&#xff09;&#xff0c;设备应允许读取已清除的LBAs而不产生错误&#xff0c;并在最后一次清除完成后&#xff0c;对未写入LBAs的读取返回所有零值给主机 NVMe-AD…

什么牌子的护眼灯好用?2024好用护眼台灯分享

不良的光线、长时间的用眼都会给眼睛带来压力&#xff0c;影响视力健康&#xff01; 本人就是一个因为工作原因需要长时间坐电脑前码字和P图的打工人&#xff0c;对于出现眼睛酸痛、疲劳以及眼球出现红血丝的情况有多难受我是深有体会&#xff0c;在此之前我搜索了好多缓解眼睛…

golang并发编程-channel

在golang 并发编程里&#xff0c;经常会听到一句话&#xff1a;不要通过共享内存进行通信&#xff0c;通过通信来共享内存。下面我们会介绍下channel, 通过源码的方式去了解channel是怎么工作的。 基本结构 流程图 代码解读 type hchan struct {qcount uint // …

(NeRF学习)NeRFStudio安装win11

参考&#xff1a; 【深度学习】【三维重建】windows11环境配置tiny-cuda-nn详细教程nerfstudio介绍及在windows上的配置、使用NeRFStudio官网githubRuntimeError: PytorchStreamReader failed reading zip archive: failed finding central directory原因及解决 目录 requireme…

不同角度深入探讨Maya和Blender这两款软件的差异

当我们面对三维建模软件的选择时&#xff0c;许多初学者可能会感到迷茫。今天&#xff0c;我们将从不同角度深入探讨Maya和Blender这两款软件的差异&#xff0c;特别是对于游戏建模领域的用户来说&#xff0c;这将有助于您更好地理解两者之间的区别。 软件授权与开发背景&#…

Python爬虫中的协程

协程 基本概念 协程&#xff1a;当程序执行的某一个任务遇到了IO操作时&#xff08;处于阻塞状态&#xff09;&#xff0c;不让CPU切换走&#xff08;就是不让CPU去执行其他程序&#xff09;&#xff0c;而是选择性的切换到其他任务上&#xff0c;让CPU执行新的任务&#xff…

引导过程的解析以及教程za

bios加电自检------mbr--------grub-------加载内核文件------启动第一个进程 bios的主要作用&#xff1a;检测硬件是否正常&#xff0c;然后根据bios中的启动项设置&#xff0c;去找内核文件 boot开机启动项顺序&#xff0c;你可以把内核文件放在何处&#xff1f; 1.硬盘 …

MySQL将多条数据合并成一条的完整示例

数据库中存的是多条数据&#xff0c;展示的时候需要合并成一条 数据表存储形式如下图 以type分组&#xff0c;type相同的算一条&#xff0c;且保留image和link的所有数据&#xff0c;用groupBy只保留一条数据 解决方案&#xff1a;用GROUP_CONCAT 完整语法如下 group_concat…

基于YOLOv8深度学习的人脸面部表情识别系统【python源码+Pyqt5界面+数据集+训练代码】深度学习实战

《博主简介》 小伙伴们好&#xff0c;我是阿旭。专注于人工智能、AIGC、python、计算机视觉相关分享研究。 ✌更多学习资源&#xff0c;可关注公-仲-hao:【阿旭算法与机器学习】&#xff0c;共同学习交流~ &#x1f44d;感谢小伙伴们点赞、关注&#xff01; 《------往期经典推…

10 个值得收藏的顶级手机数据恢复软件【2024年最新】

手机数据恢复&#xff0c;不要担心&#xff0c;今天就给大家分享10款数据恢复软件&#xff01; 现代人的手机中存储了许多重要数据&#xff0c;如照片、视频、消息、联系人等文件&#xff0c;如果手机损坏或数据丢失&#xff0c;这是一件非常烦恼的事情。此时&#xff0c;一款好…

解决jenkins的Exec command命令不生效,或者执行停不下来的问题

Jenkins构建完后将war包通过 Publish Over SSH 的插件发布到服务器上&#xff0c;在服务器上执行脚本时&#xff0c;脚本中的 nohup 命令无法执行&#xff0c;并不生效&#xff0c;我配置的Exec command命令是后台启动一个war包&#xff0c;并输出日志文件。 nohup java -jar /…

nginx源码分析-4

这一章内容讲述nginx的模块化。 ngx_module_t&#xff1a;一个结构体&#xff0c;用于描述nginx中的各个模块&#xff0c;其中包括核心模块、HTTP模块、事件模块等。这个结构体包含了一些模块的关键信息和回调函数&#xff0c;以便nginx在运行时能够正确地加载和管理这些模块。…

《动手学深度学习》学习笔记 第5章 深度学习计算

本系列为《动手学深度学习》学习笔记 书籍链接&#xff1a;动手学深度学习 笔记是从第四章开始&#xff0c;前面三章为基础知道&#xff0c;有需要的可以自己去看看 关于本系列笔记&#xff1a; 书里为了让读者更好的理解&#xff0c;有大篇幅的描述性的文字&#xff0c;内容很…

算法学习系列(十四):并查集

目录 引言一、并查集概念二、并查集模板三、例题1.合并集合2.连通块中点的数量 引言 这个并查集以代码短小并且精悍的特点&#xff0c;在算法竞赛和面试中特别容易出&#xff0c;对于面试而言&#xff0c;肯定不会让你去写一两百行的代码&#xff0c;一般出的都是那种比较短的…

[GKCTF 2020]ez三剑客-eztypecho

[GKCTF 2020]ez三剑客-eztypecho 考点&#xff1a;Typecho反序列化漏洞 打开题目&#xff0c;发现是typecho的CMS 尝试跟着创建数据库发现不行&#xff0c;那么就搜搜此版本的相关信息发现存在反序列化漏洞 参考文章 跟着该文章分析来&#xff0c;首先找到install.php&#xf…

Unable to connect to Redis server

报错内容&#xff1a; Exception in thread "main" org.redisson.client.RedisConnectionException: java.util.concurrent.ExecutionException: org.redisson.client.RedisConnectionException: Unable to connect to Redis server: 175.24.186.230/175.24.186.230…