基于矩阵乘的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.hqwc.cn/news/316115.html

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

相关文章

test dbtest-03-对比 Liquibase、flyway、dbDeploy、dbsetup

test dbtest-03-对比 Liquibase、flyway、dbDeploy、dbsetup

详细对比 Liquibase、flyway、dbDeploy、dbsetup&#xff0c;给出对比表格 下面是一个简要的对比表格&#xff0c;涵盖了 Liquibase、Flyway、dbDeploy 和 DbSetup 这四个数据库变更管理工具的一些主要特点。 特点/工具LiquibaseFlywaydbDeployDbSetup开发语言Java&#xff0…
阅读更多...
c++初阶-------类和对象

c++初阶-------类和对象

作者前言 &#x1f382; ✨✨✨✨✨✨&#x1f367;&#x1f367;&#x1f367;&#x1f367;&#x1f367;&#x1f367;&#x1f367;&#x1f382; ​&#x1f382; 作者介绍&#xff1a; &#x1f382;&#x1f382; &#x1f382; &#x1f389;&#x1f389;&#x1f389…
阅读更多...
【Spark精讲】记一个SparkSQL引擎层面的优化：SortMergeJoinExec

【Spark精讲】记一个SparkSQL引擎层面的优化：SortMergeJoinExec

SparkSQL的Join执行流程 如下图所示&#xff0c;在分析不同类型的Join具体执行之前&#xff0c;先介绍Join执行的基本框架&#xff0c;框架中的一些概念和定义是在不同的SQL场景中使用的。 在Spark SQL中Join的实现都基于一个基本的流程&#xff0c;根据角色的不同&#xff0…
阅读更多...
C# visual studio COM创建及调用

C# visual studio COM创建及调用

1.visual studuio操作 1.1visual Studio创建类库项目ClassLibrary1。 1.1.1ClassLibrary1项目Class1.cs内容如下&#xff1a; using System; using System.Collections.Generic; using System.Linq; using System.Runtime.InteropServices; using System.Text; using System.T…
阅读更多...
Spring技术内幕笔记之SpringMvc

Spring技术内幕笔记之SpringMvc

WebApplicationContext接口的类继承关系 org.springframework.web.context.ContextLoader#initWebApplicationContext 对IOC容器的初始化 SpringMvc如何设计 DispatcherServlet类继承关系 MVC处理流程图如下&#xff1a; DispatcherServlet的工作大致可以分为两个部分&#xf…
阅读更多...
使用GPTs+Actions自动获取第三方数据

使用GPTs+Actions自动获取第三方数据

目录 安装插件与GPT对话联网插件首先，创建GPTs。 Voxscript 官网：https://voxscript.awt.icu/index.htmlOpenAI Schema：https://voxscript.awt.icu/swagger/v1/swagger.yamlServer URL： servers: url: https://voxscript.awt.icu安装插件 要使用这个插件&
阅读更多...
RK3568 学习笔记 : 解决 linux_sdk 编译 python 版本报错问题

RK3568 学习笔记 : 解决 linux_sdk 编译 python 版本报错问题

前言 最近买了 【正点原子】 的 RK3568 开发板&#xff0c;下载了 开发板的资料&#xff0c;包括 Linux SDK&#xff0c;这个 Linux SDK 占用的空间比较大&#xff0c;扩展了一下 VM 虚拟机 ubuntu 20.04 的硬盘空间&#xff0c;编译才正常通过。 编译 RK3568 Linux SDK 时&am…
阅读更多...
文件监控软件丨文件权限管理工具

文件监控软件丨文件权限管理工具

文件已经成为企业最重要的资产之一。然而&#xff0c;文件的安全性和完整性经常受到威胁&#xff0c;如恶意软件感染、人为误操作、内部泄密等。 为了确保文件的安全&#xff0c;文件监控软件应运而生。本文将深入探讨文件监控软件的概念、功能、应用场景和未来发展等方面。 文…
阅读更多...
解决jenkins的Exec command命令不生效，或者执行停不下来的问题

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

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

webpack 5 loader

webpack 本身不能识别js&#xff0c;json外的资源&#xff0c;所以我们需要借助其他loader来处理对应的文件 CSS Loader&#xff0c;处理css 安装 npm i css-loader style-loader -D css-loader 负责讲css编译成webpack能识别的模块内容style-loader 动态创建<style&g…
阅读更多...
分类预测 | Matlab实现DBO-SVM蜣螂算法优化支持向量机多特征分类预测

分类预测 | Matlab实现DBO-SVM蜣螂算法优化支持向量机多特征分类预测

分类预测 | Matlab实现DBO-SVM蜣螂算法优化支持向量机多特征分类预测 目录 分类预测 | Matlab实现DBO-SVM蜣螂算法优化支持向量机多特征分类预测分类效果基本描述程序设计参考资料 分类效果 基本描述 1.Matlab实现DBO-SVM蜣螂算法优化支持向量机多特征分类预测&#xff08;完整…
阅读更多...
计算机网络【Cookie和session机制】

计算机网络【Cookie和session机制】

会话&#xff08;Session&#xff09;跟踪是Web程序中常用的技术&#xff0c;用来跟踪用户的整个会话。常用的会话跟踪技术是Cookie与Session。Cookie通过在客户端记录信息确定用户身份&#xff0c;Session通过在服务器端记录信息确定用户身份。 本章将系统地讲述Cookie与Sess…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023