高性能计算-SGEMM矩阵乘法(29)
1. 介绍
- 矩阵A（MK） B（KN）单精度浮点数进行矩阵乘法；
- 分别实现CPU串行，GPU多种并行计算算法，与 cublas 库 sgemm 函数效率对比。
2. gpu 并行算法简介
- 并行算法一：二维block,一个线程程计算一个C元素,缺点:访存次数过多
- 并行算法二:
优化一: 使用线程块 tile(将线程块计算数据划分tile) + 共享内存,线程块大小与C每个局部结果的大小相同
优化二: 在K方向划分数据块,降低共享内存的大小
优化三: 使用函数模板,减少寄存器数量
缺点: 块内线程仍有重复访存
- 并行算法三: block共享内存 + 线程tile,block共享内存再划分线程tile,一个线程计算一块线程tile数据
简介：
- A B矩阵分别按行列划分数据块,防止共享内存太大,在K方向划分数据块作为block要处理的数据,在K方向按步长 BK=8 划分数据块,A的block计算数据大小为 BM * BK=128 * 8;
- block 计算结果矩阵大小为 BMBN,每个线程计算数据大小为TYTX,blocksSize: (BM/TY) * (BN/TX)=16 * 16=256,gridSize:(M/BM,N/BN);
- 从全局内存取数据：每个block计算数据块用所有线程分别从A B矩阵取数据(使用float4,每个线程每次取 4个 float)，计算取数据到共享内存需要迭代的次数：BM * BK/(blockDim.y * blockDim.x);
- 大循环计算：在 K 方向上以 BK 为步长循环，将每次block计算的部分结果相加，最终得到最终 C 局部结果；
- 大循环中线程小循环
- 线程tile数据初始化：块内线程（threadIdx.x,threadIdx.y）初始化取 A全局内存中行索引为 blockIdx.y * blockDim.y+threadIdx.y * TY 到 offset + TY 的数据，列索引看在哪一次大循环中；初始化取 B全局内存中列索引为 blockIdx * blockIdm.x + threadIdx.x * TX 到 offset + TX 的数据，行索引看那一次大循环中。一个线程共需要初始化 A和B中各 TY * TX = 8 * 8的数据；
- 线程计算：一个线程tile计算的数据为 TY * TX = 8 * 8,共 TX 个小迭代。线程tile内小迭代每轮一个线程取 1 * TX=1 * 8 的数据。
优化一: 在K方向上共享内存可能太大,在K方向上划分 block tile 循环迭代；
优化二: 从全局内存取数据用寄存器中转,增加独立指令,提高并行性；
优化三: 矩乘使用外积,减少访存次数;每轮大迭代本线程计算线程tile中 TYTX 的数据,每轮大迭代结果相加得到最终本线程计算的 C中 TYTX矩阵的结果
优化四: 每个线程每轮小迭代处理 81的数据,进行转置处理为 18,增加访存缓存命中率,减少访存次数
优化五: 共享内存及线程计算数据寄存器double buffer的作用:读写分离,访存与计算 overlap
优化六: 预读第一轮大循环共享内存数据,大循环把warp1最后一轮小循环计算提出放最后,实现与 warp2 大循环刚开始的访问共享内存预读下轮计算数据overlap。
优化七: 大循环外预读本线程第一轮小循环寄存器数据,小循环每次迭代不同warp的计算和下一轮小迭代计算数据的访存实现 overlap
待优化:由于block每行线程数量为16,此处有 warpSize/(32/TX)= 8路 bank conflict 的问题,考虑分 bank 读取降低 bank conflict,但该方法不能完全解决问题。
NOTICE 注意: 寄存器的 double buffer 在一个warp中的读写可以分离(在一个warp中的不同线程的读写指令可以并行),但是一个线程中串行的读写指令无法并行。
3. 源码
- common.h
点击查看代码
#include <stdlib.h>
#include <sys/time.h>
#include <stdio.h>
// #define CUDA_ERROR_CHECK
#define PRINTF_LINE printf("%d final\n",__LINE__);\
fflush(stdout)
//宏定义检查API调用是否出错
#define CudaSafeCall(err) __cudaSafeCall(err,__FILE__,__LINE__)
inline void __cudaSafeCall(cudaError_t err,const char* file,const int line)
{
#ifdef CUDA_ERROR_CHECK
if(err!=cudaSuccess)
{
fprintf(stderr,"cudaSafeCall failed at %s:%d :(%d) %s\n",file,line,err,cudaGetErrorString(err));
exit(-1);
}
#endif
}
//宏定义检查获取流中的执行错误，主要是对核函数
#define CudaCheckError() _cudaCheckError(__FILE__,__LINE__)
inline void _cudaCheckError(const char * file,const int line)
{
#ifdef CUDA_ERROR_CHECK
cudaError_t err = cudaGetLastError();
if(err != cudaSuccess)
{
fprintf(stderr,"cudaCheckError failed at %s:%d :(%d) %s\n",file,line,err,cudaGetErrorString(err));
exit(-1);
}
#endif
}
//CPU时间，用于串行实现计时 s
double cpuSecond()
{
struct timeval tp;
gettimeofday(&tp, NULL);
return ((double)tp.tv_sec + (double)tp.tv_usec * 1e-6);
}
void CheckResult(float *arr,float *arrD,int len)
{
for(int i=0;i<len;i++)
{
if(fabs(arr[i] - arrD[i])>1e-6)
{
printf("第%d个数值不相等,arr:%.3f arrD:%.3f\n",i,arr[i],arrD[i]);
exit(0);
}
}
}
- gemm.cu
/*
对单精度方阵矩阵乘运算 A:M*K B:K*N
备注:编译参数 -lcubasl -arch=sm_XX(根据显卡实际参数填写)
执行参数: 0:kernel_SgemmBlock 1:kernel_SgemmBlockTileShared 3:kernel_SgemmBlockTileShared
*/
#include "common.h"
#include <cublas_v2.h>
//数据类型转换的优先级高于 []
#define GET_FLOAT4(pointer) (reinterpret_cast<float4*>(pointer)[0]) //转换指针用 float4 访问数据
#define OFFSET2D(row,col,width) ((row)*(width)+col) //根据二维行列索引计算一维索引
//并行算法一: 二维block,一个线程程计算一个C元素,缺点:访存次数过多
__global__ void kernel_SgemmBlock(float *A,float *B,float *C,int M,int K,int N)
{
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if(x < N && y <M)
{
float temp=0.0;
for(int k=0;k<K;k++)
temp += A[y*K+k] * B[k*N+x];
C[y*N+x] = temp;
}
}
//并行算法二:
//优化一: 使用线程块 tile(将线程块计算数据划分tile) + 共享内存,线程块大小与C每个局部结果的大小相同
//优化二: 在K方向划分数据块,降低共享内存的大小
//优化三: 使用函数模板,编译阶段可以优化代码
//缺点:块内线程仍有重复访存
template<
const int TileBM, //A 矩阵tile高度
const int TileBK, //A 矩阵tileK方向宽度
const int TileBN //B 矩阵tile宽度
>
__global__ void kernel_SgemmBlockTileShared(float *A,float *B,float *C,int M,int K,int N)
{
int idy = blockIdx.y * blockDim.y + threadIdx.y;
int idx = blockIdx.x * blockDim.x + threadIdx.x;
//共享内存初始化
__shared__ float SharedA[TileBM][TileBK];
__shared__ float SharedB[TileBK][TileBN];
float temp = 0;
for(int k=0;k<K/TileBK;k++)
{
//块内每个线程初始化一组 A B数据
SharedA[threadIdx.y][threadIdx.x] = A[idy*K+k*TileBK+threadIdx.x];
SharedB[threadIdx.y][threadIdx.x] = B[(k*TileBK+threadIdx.y)*N+idx];
__syncthreads();
for(int s=0;s<TileBK;s++)
temp += SharedA[threadIdx.y][s] * SharedB[s][threadIdx.x];
//tile 数据块计算要同步,避免共享内存数据错误
__syncthreads();
}
C[idy * N +idx] = temp;
}
/*并行算法三: block共享内存 + 线程tile,block共享内存再划分线程tile,一个线程计算一块线程tile数据
A B矩阵分别按行列划分数据块,防止共享内存太大,在K方向划分数据块作为block要处理的数据,在K方向按步长 BK=8 划分数据块,A的block计算数据大小为 BM * BK=128 * 8;
block 计算结果矩阵大小为 BM*BN,每个线程计算数据大小为TY*TX,blocksSize: (BM/TY) * (BN/TX)=16 * 16=256,gridSize:(M/BM,N/BN);
从全局内存取数据：每个block计算数据块用所有线程分别从A B矩阵取数据(使用float4,每个线程每次取 4个 float)，计算取数据到共享内存需要迭代的次数：BM * BK/(blockDim.y * blockDim.x);
大循环计算：在 K 方向上以 BK 为步长循环，将每次block计算的部分结果相加，最终得到最终 C 局部结果；
大循环中线程小循环
> 线程tile数据初始化：块内线程（threadIdx.x,threadIdx.y）初始化取 A全局内存中行索引为 blockIdx.y * blockDim.y+threadIdx.y * TY 到 offset + TY 的数据，列索引看在哪一次大循环中；\
初始化取 B全局内存中列索引为 blockIdx * blockIdm.x + threadIdx.x * TX 到 offset + TX 的数据，行索引看那一次大循环中。一个线程共需要初始化 A和B中各 TY * TX = 8 * 8的数据；
> 线程计算：一个线程tile计算的数据为 TY * TX = 8 * 8,共 TX 个小迭代。线程tile内小迭代每轮一个线程取 1 * TX=1 * 8 的数据。
优化一: 在K方向上共享内存可能太大,在K方向上划分 block tile 循环迭代
优化二: 从全局内存取数据用寄存器中转,增加独立指令,提高并行性
优化三: 矩乘使用外积,减少访存次数;每轮大迭代本线程计算线程tile中 TY*TX 的数据,每轮大迭代结果相加得到最终本线程计算的 C中 TY*TX矩阵的结果
优化四: 每个线程每轮小迭代处理 8*1的数据,进行转置处理为 1*8,增加访存缓存命中率,减少访存次数
优化五: 共享内存及线程计算数据寄存器double buffer的作用:读写分离,访存与计算 overlap
优化六: 预读第一轮大循环共享内存数据,大循环把warp1最后一轮小循环计算提出放最后,实现与 warp2 大循环刚开始的访问共享内存预读下轮计算数据overlap。
优化七: 大循环外预读本线程第一轮小循环寄存器数据,小循环每次迭代不同warp的计算和下一轮小迭代计算数据的访存实现 overlap
待优化一:由于block每行线程数量为16,此处有 warpSize/(32/TX)= 8路 bank conflict 的问题,考虑分 bank 读取降低 bank conflict,但该方法不能完全解决问题
NOTICE 注意: 寄存器的 double buffer 在一个warp中的读写可以分离(在一个warp中的不同线程的读写指令可以并行),但是一个线程中串行的读写指令无法并行
*/
template<
const int BM, //线程块计算A矩阵数据高度 128
const int BK, //线程块计算数据K方向宽度 8
const int BN, //线程块计算B矩阵数据宽度 128
const int TY, //单个线程计算A矩阵数据高度 8
const int TX //单个线程计算B矩阵数据宽度 8
>
__global__ void kernel_SgemmThreadTile(float *A,float *B,float *C,int M,int K,int N)
{
int idx = threadIdx.x;
int idy = threadIdx.y;
int id = idy * blockDim.x + idx; //块内id
int bidx = blockIdx.x;
int bidy = blockIdx.y;
//共享内存,存放转置BlockA 的计算数据
__shared__ float AS[2][BK][BM];
__shared__ float BS[2][BK][BN];
//计算blocksize 大小
const int BRows = BM/TY;
const int BCols = BN/TX;
/*
由于在K方向数据分块，这里先从全局内存取第一次block计算数据到寄存器再到共享内存
一个线程用float4指令一次取 4个float数据,每行 BK=8 个数据需要 BK/4=2 个线程取数据
防止BM过大所有线程一次取不完数据,需要循环读取,计算BM 方向需要循环次数 BM*BK/(BRows*BCols*4)=1,
*/
// 先取矩阵A
//所有线程每取一轮数据在BM方向的步长 =128
const int strideBM = (BRows*BCols*4)/BK;
//所有线程取一轮数据在B矩阵K 方向上的步长 = 8
const int strideBK = (BRows*BCols*4)/BN;
//定义从全局内存取数据中转的寄存器
float g2rA[BM/strideBM * 4]; // global memory to register
float g2rB[BK/strideBK * 4];
//每个循环本线程取数据起始位置: flot4 数据在tile的strideBM 行内的起始行列索引
int rowStartA = id/(BK/4);
int colStartA = id%(BK/4)*4;
#pragma unroll
for(int i=0;i<BM*BK/(BRows*BCols*4);i++)
{
GET_FLOAT4(g2rA+i*4) = GET_FLOAT4(A + ((bidy*BM +i*strideBM + rowStartA)*K + 0 + colStartA));
//转置放入共享内存,第一个float原索引为 [i*strideBM + rowStartA][colStartA]
AS[0][colStartA][i*strideBM + rowStartA] = g2rA[i*4+0];
AS[0][colStartA + 1][i*strideBM + rowStartA] = g2rA[i*4+1];
AS[0][colStartA + 2][i*strideBM + rowStartA] = g2rA[i*4+2];
AS[0][colStartA + 3][i*strideBM + rowStartA] = g2rA[i*4+3];
}
//再取矩阵B到共享内存
//flot4 数据在strideBK 行内的起始行列索引
int rowStartB = id/(BN/4);
int colStartB = id%(BN/4)*4;
#pragma unroll
for(int i=0;i<(BN*BK)/(BRows*BCols*4);i++)
{
//无转置,只初始化一轮,省略寄存器中转
GET_FLOAT4(&BS[0][rowStartB][colStartB]) = GET_FLOAT4(B + (0 + i*strideBK + rowStartB)*N + bidx*BN + colStartB);
}
//共享内存已初始化,块内数据同步
__syncthreads();
//当前线程小迭代计算第一次数据准备
//第一次大循环前每个线程从共享内存取出各自线程小迭代第一轮的计算数据放入临时寄存器,每个线程取 bk=8个数据;
//双buffer作用: 计算与访存overlap,第一次小循环时本轮计算与下轮小循环的数据预读在不同warp上重叠
float tempA[2][TY]; //2*8*1
float tempB[2][TX]; //2*1*8
#pragma unroll
for(int s=0;s<TY;s+=4)
GET_FLOAT4(&tempA[0][s]) = GET_FLOAT4(&AS[0][0][idy*TY+s]);
#pragma unroll
for(int s=0;s<TX;s+=4)
GET_FLOAT4(&tempB[0][s]) = GET_FLOAT4(&BS[0][0][idx*TX+s]);
//初始化共享内存double buffer 标志
int WFlag = 1;
// 大循环,K 方向遍历结果相加,最终得到块的最终结果
float tempC[TY][TX] = {0}; //存放本轮循环本线程计算结果
for(int k=0; k<K; k+=BK)
{
//1.对下一轮大循环,每个线程先从全局内存预读取block数据到寄存器;在每个block的每个 stride 内 colStart 相对位置不变
// colStartA += BK;
// rowStartB += BK;
//FIXME:ok. k+BK,而不是 (k+1)*BK
#pragma unroll
for(int i=0;i<BM*BK/(BRows*BCols*4);i++) //条件为 stride 的个数
GET_FLOAT4(g2rA+i*4) = GET_FLOAT4(A + ((bidy*BM +i*strideBM + rowStartA)*K + (k+BK) + colStartA));
//再取矩阵B
//这里从B中取数据为了提高并行性使用了独立指令,用寄存器做数据中转
#pragma unroll
for(int i=0;i<(BN*BK)/(BRows*BCols*4);i++)
GET_FLOAT4(g2rB+i*4) = GET_FLOAT4(B + ((k+BK) + i*strideBK + rowStartB)*N + bidx*BN + colStartB);
WFlag = 1-(k/BK)%2;
int RFlag = 1 - WFlag;
//2.小迭代,块内计算
// 进行 TY-1次计算
//由于每轮大迭代开始需要初始化下一轮的共享内存,延迟高,所以这里进行TY-1次小迭代计算,将最后一次小迭代与下一轮大循环的初始化共享内存并行
#pragma unroll
for(int r=0;r<TY-1;r++)
{
//2.1 从共享内存预读取下一轮小迭代计算数据到寄存器写缓冲区
//A矩阵 按float4读取,因为线程块大小是根据 BM/TY BN/TX计算的,所以每个线程每轮小迭代取 TY(2个float4)个计算数据 (B一样)
//将一个block计算数据划分为 BM/TY=128/8=16个 tile,
//每个线程根据线程的 idx idy 判断初始化共享内存中哪个 tileA tileB:
//确定每个线程要初始化的 tile数据:索引为 id 的线程在 AS中取tile索引为 idy 的数据,在 BS中取tile索引为 idx 的数据;
//第 i 次小迭代,取 tileA 第 i 列 TY*1=8*1;与 tileB 第 i 行 1*TX=1*8的数据做外积
//TODO 由于block每行线程数量为16,此处有 warpSize/(32/TX)= 8路 bank conflict 的问题,考虑分 bank 读取降低 bank conflict,但该方法似乎不能完全解决问题
//TODO ok,此处寄存器双buffer的用处??一次循环中读缓冲区读数据做计算，写缓冲区写数据，实现 warp2 访问共享内存与 warp1 计算的重叠。
#pragma unroll
for(int s=0;s<TY;s+=4)
GET_FLOAT4(&tempA[1-r%2][s]) = GET_FLOAT4(&AS[RFlag][r+1][idy*TY+s]); //这里对共享内存读,小迭代寄存器写入
#pragma unroll
for(int s=0;s<TX;s+=4)
GET_FLOAT4(&tempB[1-r%2][s]) = GET_FLOAT4(&BS[RFlag][r+1][idx*TX+s]);
//2.2 本轮小迭代计算,从寄存器读缓冲区取数据
#pragma unroll
for(int i=0;i<TY;i++)
{
#pragma unroll
for(int j=0;j<TX;j++)
tempC[i][j] += tempA[r%2][i] * tempB[r%2][j]; //计算 读寄存器
}
}
//3.第 1 步的中转寄存器数据 g2r 转移到共享内存(下一轮大循环初始化共享内存)
#pragma unroll
for(int i=0;i<BM*BK/(BCols*BRows*4);i++)
{
//这里对共享内存写入与小迭代对共享内存的读取,在不同线程束间读写分离,实现线程束级并行
AS[WFlag][colStartA][i*strideBM + rowStartA] = g2rA[i*4];
AS[WFlag][colStartA+1][i*strideBM + rowStartA] = g2rA[i*4+1];
AS[WFlag][colStartA+2][i*strideBM + rowStartA] = g2rA[i*4+2];
AS[WFlag][colStartA+3][i*strideBM + rowStartA] = g2rA[i*4+3];
}
#pragma unroll
for(int i=0;i<BN*BK/(BRows*BCols*4);i++)
GET_FLOAT4(&BS[WFlag][i*strideBK+rowStartB][colStartB]) = GET_FLOAT4(g2rB + i*4);
//4.共享内存,块内数据同步
__syncthreads();
// 5.取出下轮大迭代中小迭代第一次的计算数据到寄存器
//小迭代这里对寄存器写入,
#pragma unroll
for(int s=0;s<TY;s+=4)
GET_FLOAT4(&tempA[0][s]) = GET_FLOAT4(&AS[WFlag][0][idy*TY + s]);
#pragma unroll
for(int s=0;s<TX;s+=4)
GET_FLOAT4(&tempB[0][s]) = GET_FLOAT4(&BS[WFlag][0][idx*TX + s]);
// 6.小迭代最后一轮计算,与其他 warp 下一轮大迭代中第一步访问全局内存到寄存器重叠
//NOTICE 注意: 寄存器的 double buffer 在一个warp中的读写可以分离(在一个warp中的不同线程的读写指令可以并行),但是一个线程中串行的读写指令无法并行
#pragma unroll
for(int i=0;i<TY;i++)
{
#pragma unroll
for(int j=0;j<TX;j++)
tempC[i][j] += tempA[1][i] * tempB[1][j];
}
}
//把本线程block 计算结果(8*8)保存到 C
#pragma unroll
for(int i=0;i<TY;i++)
{
#pragma unroll
for(int j=0;j<TX;j+=4)
GET_FLOAT4(C + (bidy*BM+idy*TY+i)*N + bidx*BN+idx*TX+j) = GET_FLOAT4(&tempC[i][j]);
}
}
int main(int argc,char** argv)
{
int mode = 0;
int M,N,K;
M=N=K=4096;
if(argc>1)
mode = atoi(argv[1]);
float *A,*B,*C;
A = (float*)calloc(M*K,sizeof(float));
B = (float*)calloc(K*N,sizeof(float));
C = (float*)calloc(M*N,sizeof(float));
//矩阵初始化
for(int i=0;i<M*K;i++)
A[i] = i/10;
for(int i=0;i<K*N;i++)
B[i] = i%10;
double cpuTime = 0.0;
#if 0
double start = cpuSecond();
for(int m=0;m<M;m++)
{
for(int n=0;n<N;n++)
{
float temp = 0;
for(int k=0;k<K;k++)
temp += A[m*K+k] * B[k*N+n];
C[m*N+n] = temp;
}
}
cpuTime = (cpuSecond()-start)*1000;
#endif
float *AD,*BD,*CD,*CFD; //cfd 主存保存显存结果
cudaMalloc((void**)&AD,M*K*sizeof(float));
cudaMalloc((void**)&BD,K*N*sizeof(float));
cudaMalloc((void**)&CD,M*N*sizeof(float));
CFD = (float*)calloc(M*N,sizeof(float));
cudaMemcpy(AD,A,M*K*sizeof(float),cudaMemcpyHostToDevice);
cudaMemcpy(BD,B,K*N*sizeof(float),cudaMemcpyHostToDevice);
//cublas设置
cublasHandle_t blas;
cublasCreate(&blas);
float alpha = 1.0,beta=0;
float gpuTime =0;
cudaEvent_t startD,endD;
cudaEventCreate(&startD);
cudaEventCreate(&endD);
// cudaEventRecord(startD);//FIXME 计时应放在条件分支预测后面,否则时间增加大约 100ms
if(mode == 0) //二维block
{
cudaEventRecord(startD);
int width,height;
width = height = 16;
dim3 BlockSize(width,height);
dim3 GridSize((N-1)/BlockSize.x+1,(M-1)/BlockSize.y+1);
kernel_SgemmBlock<<<GridSize,BlockSize>>>(AD,BD,CD,M,K,N);
}
else if (mode == 1) //block tile+共享内存
{
cudaEventRecord(startD);
int width,height;
width = height = 16;
const int TileBM = 16;
const int TileBK = 16;
const int TileBN = 16;
dim3 BlockSize(width,height);
dim3 GridSize((N-1)/BlockSize.x+1,(M-1)/BlockSize.y+1);
kernel_SgemmBlockTileShared<TileBM, TileBK, TileBN><<<GridSize, BlockSize>>>(AD,BD,CD,M,K,N);
}
else if (mode == 2)//线程tile,双缓冲区
{
cudaEventRecord(startD);
//以下为参考的经验参数,性能较好
const int BM = 128;
const int BN = 128;
const int BK = 8;
const int TY = 8;
const int TX = 8;
dim3 BlockSize(BM/TY,BN/TX); //16*16
//FIXME ok.error: dim3 GridSize(M / BlockSize.y, N / BlockSize.x);
dim3 GridSize(M/BM,N/BN);
kernel_SgemmThreadTile<BM,BK,BN,TY,TX><<<GridSize,BlockSize>>>(AD,BD,CD,M,K,N);
}
else if(mode ==3)//cublas
{
cudaEventRecord(startD);
//数据按列存储的
//转置运算, C = A * B ,C 是按列存储的,后面需要对结果再转置
cublasSgemm(blas,CUBLAS_OP_T,CUBLAS_OP_T,M,N,K,&alpha,AD,K,BD,N,&beta,CD,N);
//不转置计算 C^T=(AB)^T=B^T * A^T C^T为按列存储的最终结果
// cublasSgemm(blas,CUBLAS_OP_N,CUBLAS_OP_N,N,M,K,&alpha,BD,N,AD,M,&beta,CD,K);
}
// CudaCheckError();
cudaEventRecord(endD);
cudaEventSynchronize(endD);
cudaEventElapsedTime(&gpuTime,startD,endD);
cudaEventDestroy(startD);
cudaEventDestroy(endD);
cublasDestroy(cublas);
cudaMemcpy(CFD,CD,M*N*sizeof(float),cudaMemcpyDeviceToHost);
#if 1
if(mode==3) //计算参数为 CUBLAS_OP_T需要对结果(按列存储)转置(按行存储)
{
float *CFDT = (float*)calloc(M*N,sizeof(float));
memcpy(CFDT,CFD,M*N*sizeof(float));
#pragma unroll
for(int i=0;i<M;i+=4) //行列分块增加缓存命中
{
#pragma unroll
for(int j=0;j<N;j+=4)
{
#pragma unroll
for(int m=0;m<4;m++)
{
#pragma unroll
for(int n=0;n<4;n++)
*(CFDT+(j+n)*M+i+m) = *(CFD+(i+m)*N+j+n);
}
}
}
memcpy(CFD,CFDT,M*N*sizeof(float));
free(CFDT);
}
#endif
printf("cpu串行耗时 %.3f ms,gpu:耗时 %.3f ms,加速比 %.3f\n",cpuTime,gpuTime,cpuTime/gpuTime);
//检查结果
CheckResult(C,CFD,M*N);
free(A);
free(B);
free(C);
free(CFD);
cudaFree(AD);
cudaFree(BD);
cudaFree(CD);
return 0;
}
- 备注: 编译参数 -lcubasl -arch=sm_XX(根据显卡实际参数填写，否则计算性能不能发挥最大)
4. 测试数据
- 4096*4096 矩阵测试数据如下
| 计算方法 | 耗时 ms | 与 cublasSgemm 浮点计算能力比较 |
|---|---|---|
| cpu 串行 | 596887 | 0 |
| kernel_SgemmBlock | 1537 | 0.47% |
| kernel_SgemmBlockTileShared | 641 | 11.34% |
| kernel_SgemmThreadTile | 117 | 62.39% |
| cublasSgemm | 73 | 100% |
5. 总结
- 本代码性能最好的 kernel_SgemmThreadTile 与练习参考代码仍有10%的效率差距，后期可以通过 Nsight compute 进行分析继续改进。
- 待优化:由于block每行线程数量为16,此处有 warpSize/(32/TX)= 8路 bank conflict 的问题,考虑分 bank 读取降低 bank conflict,但该方法不能完全解决问题。
- kernel_SgemmThreadTile 线程tile实现的过程中遇到了一些bug问题，现将调试经过记录一下。
问题1：选用 128*128 的矩阵进行调试，结果检查显示第128个数据计算错误?
- 使用 Nvidia Visual Stidio Edition 插件在Visual Stidio 中调试: 由于每个线程计算一块数据，最后赋值给结果数据C，调试中发现线程 0 最后计算结果正确，但是在 kernel 外打印数据错误，数据被污染，考虑到每个线程计算自己的数据，并且存放数据的索引再三确认无误，不应出现此问题。
- 次日想到，在调试过程中有出现焦点跑到别的 block的情况，才发现一点问题。128*128的数据量理分配一个block计算，考虑 grideSize 有错误，经检查确实发现问题。
- 原因分析：grideSize 错误的本质是对该并行算法并未深入理解，该部分代码编写较早，先期只是有个大概认识，后期在编写 kernel 代码时才慢慢了解了该算法，后期没有检查 grideSize 参数。
问题2：测试 cublasSgemm 代码GPU耗时时发现比参考代码中对 cublasSgemm 的测试多了100ms?
- 经排查发现，因为多个并行算法在一个 cu 文件中，为了测试流程控制使用了条件分支，gpu 计时放在了条件分支外，导致cpu 条件分支预测错误性能损耗。
浙公网安备 33010602011771号