基于矩阵乘的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.mzph.cn/news/593074.shtml

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

相关文章

C/C++汇编学习(一)——C++代码到汇编代码的转换

目录 一、 学习C编译过程 预处理器的作用 编译器如何将C代码转换为汇编 汇编代码示例&#xff08;假设&#xff09; 链接器的作用 二、 使用编译器探索C代码的汇编表示 使用GCC查看汇编输出 使用Visual Studio查看汇编输出 比较不同优化级别的汇编代码 三、 理解C控制…

判断两张图片是否完全一致

判断两张图片是否为完全相同的图片 批量判断尺寸 大小 图像展示内容体是否完全一致的图片 import os import hashlib from PIL import Imagedef check_img_repeat(directory):"""批量对图片进行重复性校验是检查一组图像中是否有相同或几乎相同的图像副本。一…

Linux C/C++ 获取CPUID

实现方式&#xff1a; INTEL CC 格式 AT^T CC 格式 GCC/C库 __cpuid 宏 大致讲义&#xff1a; AT^T 格式汇编很反人类&#xff0c;GCC可以改编译器选项为INTEL内嵌汇编&#xff0c;但一般在GCC还是按照默认的AT^T汇编来拽写把&#xff0c;不想用也可以让AI工具把INTEL内嵌…

神经网络模型常见的类型介绍

神经网络模型有很多种类型&#xff0c;下面是一些常见的神经网络模型&#xff1a; 前馈神经网络&#xff08;Feedforward Neural Network&#xff09;&#xff1a;这是最简单、最基础的神经网络形式。在这种网络中&#xff0c;信息只朝一个方向流动&#xff0c;从输入节点经过…

Redis命令---String篇 (超全)

目录 1.Redis Setnx 命令 - 只有在 key 不存在时设置 key 的值。简介语法可用版本: > 1.0.0返回值: 设置成功&#xff0c;返回 1 。 设置失败&#xff0c;返回 0 。 示例 2.Redis Getrange 命令 - 返回 key 中字符串值的子字符简介语法可用版本: > 2.4.0返回值: 截取得到…

2024年个人工作计划怎么写?新年待办计划这样写更方便

元旦的钟声还在耳边回响&#xff0c;2024年的新篇章已经开启。面对新的一年&#xff0c;我深知一个清晰、实用的个人工作计划是多么重要。它不仅是指引我前进的灯塔&#xff0c;更是我实现目标、提升效率的秘密武器。 但如何制定这样一个计划呢&#xff1f;在过去&#xff0c;…

嵌入式开发——ADC开发

学习目标 了解ADC开发流程掌握采样方式能够使用ADC进行芯片内部通道进行采样能够使用ADC对外部电路进行采样学习内容 GD32F4的ADC 特点: 16个外部模拟输入通道;1个内部温度传感通道(VSENSE);1个内部参考电压输入通道(VREFINT);1个外部监测电池VBAT供电引脚输入通道。ADC开…

《工具录》nslookup

工具录 1&#xff1a;nslookup2&#xff1a;选项介绍3&#xff1a;示例 - 命令行模式3.1&#xff1a;查询类型设置3.2&#xff1a;指定 DNS 服务器 4&#xff1a;示例 - 交互模式5&#xff1a;其他 本文以 kali-linux-2023.3-vmware-amd64 为例。 1&#xff1a;nslookup nsloo…

Python计算三角形的面积

Python 计算三角形的面积 以下实例为通过用户输入三角形三边长度&#xff0c;并计算三角形的面积&#xff1a; # 三角形第一边长 a 3 # 三角形第二边长 b 4 c float( input("输入三角形第三边长: ") ) # 计算半周长 s (a b c) / 2 # 计算…

算法29:不同路径问题(力扣62和63题)--针对算法28进行扩展

题目&#xff1a;力扣&#xff08;LeetCode&#xff09;官网 - 全球极客挚爱的技术成长平台 一个机器人位于一个 m x n 网格的左上角 &#xff08;起始点在下图中标记为 “Start” &#xff09;。 机器人每次只能向下或者向右移动一步。机器人试图达到网格的右下角&#xff0…

什么是安全信息和事件管理(SIEM),有什么用处

安全信息和事件管理&#xff08;SIEM&#xff09;对于企业主动识别、管理和消除安全威胁至关重要。SIEM 解决方案采用事件关联、AI 驱动的异常检测以及机器学习驱动的用户和实体行为分析 &#xff08;UEBA&#xff09; 等机制来检测、审查和应对网络安全威胁。这些功能使 SIEM …

AntDB设计之CheckPoint——引言与功能简述

1.引言 数据库服务能力提升是一项系统性的工程&#xff0c;在不同的应用场景下&#xff0c;用户对于数据库各项能力的关注点也不同&#xff0c;如&#xff1a;读写延迟、吞吐量、扩展性、可靠性、可用性等等。国内不少数据库系统通过系统架构优化、硬件设备升级等方式&#xf…

数据库课程设计报告——音乐管理系统

目录 省流版word文档需求分析系统目标业务需求及处理流程功能需求及数据需求分析业务规则分析 概念设计命名规范实体集及属性联系集及属性系统总ER图 逻辑设计关系的设计关系的优化数据库基本表设计 物理设计关系模式存取方式选择数据库的存储结构 数据库应用设计数据库脚本数据…

C语言链表、树、图的实现(结构体)

链表、树、图 链表树图邻接矩阵邻接表 链表 参看此线性表实现&#xff08;C语言——结构体&#xff09;博文 树 struct Tree{int val;struct Tree *left;struct Tree *right; };在上面的代码中&#xff0c;每一部分都是定义二叉树节点所必需的&#xff0c;所以没有多余的可以…

VS2022 Android NativeActivity 开发指南

几年前最初使用VS时&#xff0c;记得是有Android NativeActivity的&#xff0c;今天更新到了2022最新版&#xff0c;发现找不到这个创建选项。 然后确保安装了C 跨平台开发工具后&#xff0c;开始排查原因。 Visual Studio 2022 中没有“本机活动应用程序” - android - SO中…

LeetCode1534. Count Good Triplets

文章目录 一、题目二、题解 一、题目 Given an array of integers arr, and three integers a, b and c. You need to find the number of good triplets. A triplet (arr[i], arr[j], arr[k]) is good if the following conditions are true: 0 < i < j < k < …

Doris数仓开发规范

文章目录 一、字符集规范二、建表规范三、数据变更规范四、数据查询规范结尾 一、字符集规范 【强制】数据库字符集指定utf-8&#xff0c;并且只支持utf-8。 二、建表规范 【建议】库名统一使用小写方式&#xff0c;中间用下划线&#xff08;_&#xff09;分割&#xff0c;长…

Android 车联网——CarPackageManagerService介绍(十一)

CarPackageManagerService 主要用于车上使用场景扩充了一些包管理相关的接口。包括黑白名单的机制,这主要是出于安全的考虑,车上的应用有更严格的限制。结合用户体验限制对运行在Android Automotive OS 上的应用有一个更好的约束。 一、简介 CarPackageManagerService 是 An…

力扣42. 接雨水

双指针法 思路&#xff1a; 将数组前后设置为 left、right 指针&#xff0c;相互靠近&#xff1b;在逼近的过程中记录两端最大的值 leftMax、rightMax&#xff0c;作为容器的左右边界&#xff1b;更新指针规则&#xff1a; 如果数组左边的值比右边的小&#xff0c;则更新 left…

使用GO开发的IDE简介

一、IDE介绍 Goland Goland是由JetBrains公司开发的商业IDE&#xff0c;专门为Go语言开发设计。JetBrains是一家知名的软件开发公司&#xff0c;以其强大的IDE产品如IntelliJ IDEA而闻名。 优点&#xff1a; 基于IntelliJ平台&#xff0c;因此拥有与IntelliJ IDEA相似的强大功能…