parallel-processing - 矩阵乘法:在 CUDA 中合并全局内存访问后性能下降
问题描述
我最近开始使用 CUDA 处理 GPU。作为一个入门程序,我尝试有效地实现一个简单的矩阵乘法
C = AB
,从朴素的矩阵乘法开始(每个线程为 C 中的一个元素加载 A 和 B 的所有元素),平铺实现(线程在共享内存中的一个平铺中协同加载来自 A 和 B 的元素平铺以减少全局内存交通)提供良好的加速。然而,在平铺实现中,对全局内存的访问也不是按合并顺序进行的。因此,为了提高性能,最好转置矩阵 B 然后相乘。下面是我的代码,
#include<stdio.h>
#include<stdlib.h>
#include<cuda_runtime.h>
#include <time.h>
#include <sys/time.h>
void querydeviceprop();
void allocate_matrix(float *h_a, float *h_b, int matDim);
void verify(float *h_c, float *h_c_check, int matDim);
void print_matrix(float *ha, int m,int n);
void transpose_matrix(float *ha, int matDim);
void mat_mul();
#define TILE_WIDTH 16 //should be equal to numThread for tiling implementation
__global__ void MatrixMult_tiling(float *d_a,float *d_b,float *d_c, int dim){
__shared__ float ta[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
__shared__ float tb[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
int bx,by,tx,ty,i,j;
float res;
int row, col;
bx=blockIdx.x; by=blockIdx.y;
tx=threadIdx.x; ty=threadIdx.y;
row=by*TILE_WIDTH+ty;
col=bx*TILE_WIDTH+tx;
res=0;
for(i=0;i<dim/TILE_WIDTH;i++){
//collaboratively load the elements. Each thread loads a single element.
ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
tb[ty][tx]=d_b[(ty+i*TILE_WIDTH)*dim+col];
__syncthreads();
for(j=0;j<TILE_WIDTH;j++){
res=res+ta[ty][j]*tb[j][tx];
}
__syncthreads();
}
d_c[row*dim+col]=res;
}
__global__ void MatrixMult_tiling_coalesced(float *d_a,float *d_b,float *d_c, int dim){
__shared__ float ta[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
__shared__ float tb[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
int bx,by,tx,ty,i,j;
float res;
int row, col;
bx=blockIdx.x; by=blockIdx.y;
tx=threadIdx.x; ty=threadIdx.y;
row=by*TILE_WIDTH+ty;
col=bx*TILE_WIDTH+tx;
res=0;
for(i=0;i<dim/TILE_WIDTH;i++){
//collaboratively load the elements. Each thread loads a single element.
ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
__syncthreads();
for(j=0;j<TILE_WIDTH;j++){
res=res+ta[ty][j]*tb[tx][j];
}
__syncthreads();
}
d_c[row*dim+col]=res;
}
__global__ void MatrixMult_naive(float *d_a,float *d_b,float *d_c, int dim){
int row,col,i;
col=blockIdx.y*blockDim.y+threadIdx.y;
row=blockIdx.x*blockDim.x+threadIdx.x;
float res;
if(row<dim && col<dim){
res=0;
for(i=0;i<dim;i++){
res=res+(d_a[row*dim+i]*d_b[i*dim+col]);
}
d_c[row*dim+col]=res;
}
}
int main(){
mat_mul();
return 0;
}
void mat_mul(){
cudaSetDevice(0);
time_t t;
cudaError_t err = cudaSuccess;
srand((unsigned) time(&t));
cudaEvent_t start, stop;
float milliseconds=0;
cudaEventCreate(&start);
cudaEventCreate(&stop);
int matDim = 8192;
float *h_a, *h_b, *h_c, *h_c_check;
/*declare the host memories*/
h_a=(float *)malloc(matDim*matDim*sizeof(float));
h_b=(float *)malloc(matDim*matDim*sizeof(float));
h_c=(float *)malloc(matDim*matDim*sizeof(float));
h_c_check=(float *)malloc(matDim*matDim*sizeof(float));
// Verify that allocations succeeded
if (h_a == NULL || h_b == NULL || h_c == NULL || h_c_check ==NULL)
{
fprintf(stderr, "Failed to allocate host vectors!\n");
exit(EXIT_FAILURE);
}
allocate_matrix(h_a,h_b,matDim); // allocate memory to hold the matrix
//allocate cuda memory
float *d_a=NULL;
float *d_b=NULL;
float *d_c=NULL;
err=cudaMalloc((void **)&d_a, matDim*matDim*sizeof(float));
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err=cudaMalloc((void **)&d_b, matDim*matDim*sizeof(float));
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err=cudaMalloc((void **)&d_c, matDim*matDim*sizeof(float));
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
printf("Matrix dimension is : %d\n",matDim);
// Copy the host input matrix A and B in host memory to the device matrix in device memory
//printf("Copy input data from the host memory to the CUDA device\n");
cudaEventRecord(start);
err = cudaMemcpy(d_a, h_a, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
//printf("GPU memcpy matrix A %10.10f ms\n",milliseconds);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
cudaEventRecord(start);
err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
//printf("GPU memcpy matrix B %10.10f ms\n",milliseconds);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy vector B from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
/*constants for kernel launch*/
int numThread=16; //number of threads per Block axis
int numBlocks=matDim/numThread;
if(matDim%numThread)
numBlocks++;
dim3 dimGrid(numBlocks,numBlocks);
dim3 dimBlock(numThread,numThread);
//-------------------------------------------------------------
//-------transpose and copy to GPU-------
transpose_matrix(h_b, matDim);//transpose first the b matrix
err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
cudaEventSynchronize(stop);
if (err != cudaSuccess){
fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
//--------transpose and copy ends-------------
cudaEventRecord(start);
MatrixMult_tiling_coalesced<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
cudaEventRecord(stop);
err = cudaGetLastError();
if (err != cudaSuccess){
fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU time tiled & coalesced %10.10f ms\n",milliseconds);
//printf("Copy output data from the CUDA device to the host memory\n");
cudaEventRecord(start);
err = cudaMemcpy(h_c_check, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
//printf("GPU memcpy time tiled & coalesced %10.10f ms\n",milliseconds);
//------------transpose back the original B matrix----------------
transpose_matrix(h_b, matDim);//transpose first the b matrix
err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
cudaEventSynchronize(stop);
if (err != cudaSuccess){
fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
//------------transpose back the original matrix ends-------------
//-------------------------------------------------------------
cudaEventRecord(start);
MatrixMult_tiling<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
cudaEventRecord(stop);
err = cudaGetLastError();
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU time tiled %10.10f ms\n",milliseconds);
//printf("Copy output data from the CUDA device to the host memory\n");
cudaEventRecord(start);
err = cudaMemcpy(h_c, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
//printf("GPU memcpy time tiled %10.10f ms\n",milliseconds);
//-------------------------------------------------------------
/*
cudaEventRecord(start);
MatrixMult_naive<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
cudaEventRecord(stop);
err = cudaGetLastError();
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU time naive %10.10f ms\n",milliseconds);
printf("Copy output data from the CUDA device to the host memory\n");
cudaEventRecord(start);
err = cudaMemcpy(h_c, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU memcpy time naive %10.10f ms\n",milliseconds);
*/
//-------------------------------------------------------------
verify(h_c, h_c_check, matDim);
// Free device global memory
err = cudaFree(d_a);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaFree(d_b);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector B (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaFree(d_c);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector C (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Free host memory
free(h_a);
free(h_b);
free(h_c);
printf("Done\n");
}
void allocate_matrix(float *h_a, float *h_b, int matDim){
int i,j;
// Initialize the host input vectors
for (i = 0; i < matDim; i++)
{
for(j=0;j< matDim;j++){
h_a[i*matDim+j] = rand()%10;
h_b[i*matDim+j] = rand()%10;
}
}
}
void print_matrix(float *ha, int m,int n){
int i, j;
for(i=0;i<m;i++){
for(j=0;j<n;j++){
printf(" %.1f ",ha[i*m+j]);
}
printf("\n");
}
}
void transpose_matrix(float *h_a, int matDim){
int i, j;
int temp;
for(i=0;i<matDim;i++)
{
for(j=0;j<i;j++)
{
temp=h_a[i*matDim+j];
h_a[i*matDim+j]=h_a[j*matDim+i];
h_a[j*matDim+i]=temp;
}
}
}
void verify(float *h_c, float *h_c_check, int matDim){
int i,j;
//check the code
for (i = 0; i < matDim; i++)
{
for(j=0;j<matDim;j++){
if (fabs(h_c[i*matDim+j] - h_c_check[i*matDim+j]) > 1e-5)
{
printf("cpu : %f , gpu : %f\t",h_c[i*matDim+j],h_c_check[i*matDim+j]);
fprintf(stderr, "Result verification failed at element %d,%d !\n\n", i,j);
exit(EXIT_FAILURE);
}
}
}
printf("Test PASSED\n");
}
MatrixMult_tiling_coalesced
和void MatrixMult_tiling
分别是对 B 的元素有和没有合并内存访问的函数。
现在,问题是所花费的时间MatrixMult_tiling_coalesced
几乎是所花费时间的两倍MatrixMult_tiling
。我知道在MatrixMult_tiling
元素中以合并的方式(即按行主要顺序)为每个图块加载到图块中,但是图块沿列排列,而图块中的图块MatrixMult_tiling_coalesced
沿行排列,因此该功能MatrixMult_tiling_coalesced
应该比另一个更快。但在实践中,我可以看到相反的情况。如果有人能指出原因,我将不胜感激。提前谢谢..
编辑1:在罗伯特回答之后(见下文),我明白问题出在元素乘法期间的加载操作中。
tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];]
至
tb[tx][ty]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
和
res=res+ta[ty][j]*tb[tx][j];
至
res=res+ta[ty][j]*tb[j][tx];
这将MatrixMult_tiling_coalesced
函数的性能从 1500 毫秒提高到 1000 毫秒。但是,该函数MatrixMult_tiling
只需要 879 毫秒。因此,合并后的例程仍然较慢。我不明白现在问题出在哪里。
编辑 2:我意识到在编辑 1 中,我刚刚将银行冲突问题从元素乘法移到了元素加载部分。代码中的以下更改没有银行冲突,
tb[tx][ty]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
至
tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
和
res=res+ta[ty][j]*tb[j][tx];
至
res=res+ta[ty][j]*tb[ty][j];
但它仍然比功能稍慢MatrixMult_tiling
。该MatrixMult_tiling_coalesced
函数需要 982 ms 与 870 ms 的MatrixMult_tiling
函数。它应该至少类似于MatrixMult_tiling
如果不是更快的话。
最后编辑:
编辑 2 不会产生正确的结果。所以编辑 1 的代码将是最佳的。转置被乘数矩阵之一可能不是一个好主意。:-(
感谢大家的帮助。
解决方案
B
当然不是我要转置的矩阵C=AB
。但这既不是这里也不是那里。
我不确定你为什么认为:
在平铺实现中,对全局内存的访问也不是合并顺序
我没有在您的代码中看到任何MatrixMult_tiling
导致未合并访问的代码行。
只是为了确保我们不会被术语绊倒,“coalesced”或“uncoalesced”是我们应用于访问全局内存(不是共享内存)的模式的术语。您的全局内存访问模式位于平铺内核中的这些行中:
ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
tb[ty][tx]=d_b[(ty+i*TILE_WIDTH)*dim+col];
...
d_c[row*dim+col]=res;
并且这些模式都没有被合并到全局内存中。在每个生成的索引d_a
中,如果执行替换,您会发现变量存在于所有索引中,d_b
并且不乘以任何值,常量或其他值。因此,这些模式都会(很好地)合并。d_c
threadIdx.x
如果有人能指出原因,我将不胜感激。
当涉及到共享内存时,您做了一些坏事。
在您的平铺内核中,您的乘法运算如下所示:
res=res+ta[ty][j]*tb[j][tx];
对于这种情况:
ta[ty][j]
我们有一种情况,warp 中的所有线程(具有线性增加tx
的值但相同的ty
值)正在读取共享内存中的相同位置。这是一种“最佳”访问模式——它不会出现任何银行冲突,并且会在尽可能短的时间内得到服务。
对于这种情况:
tb[j][tx]
我们有一种情况,warp 中的相邻线程正在读取共享内存中的相邻位置。这也是一种“最佳”的、非银行冲突的模式,并将在尽可能短的时间内得到服务。
但是在您的MatrixMult_tiling_coalesced
内核中,相应的乘法运算是:
res=res+ta[ty][j]*tb[tx][j];
同样,在这种情况下:
ta[ty][j]
我们有一个共享内存“广播”模式(warp 中的所有线程从同一位置读取),这是最佳且快速的。但在这种情况下:
tb[tx][j]
您实际上已经创建了对共享内存的列访问。这是共享内存最糟糕的访问模式,它将导致加载过程的 32 路序列化(或者可能是 16 路序列化,在您的 16x16 线程块的情况下),并且性能肯定更差。为什么?请记住,对于给定的负载,j
在整个经线中是恒定的,并且tx
在整个经线中线性增加。因此,假设j
在特定循环迭代中为 1。经线 0 中的线程将读取:
tb[0][1], tb[1][1], tb[2][1], tb[3][1], ...
并且这些位置都属于共享内存的特定“列”,即它们都属于同一个共享内存库。这是共享内存的最坏情况模式。
为了完整起见,我声称MatrixMult_tiling_coalesced
内核中的所有全局内存访问模式也被合并:
ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
...
d_c[row*dim+col]=res;
因此,您的两个内核实现之间的全局内存访问模式/活动/效率应该没有重大差异。
作为旁注,我认为这都是一个学习练习。如果您对 GPU 上的高性能矩阵乘法感兴趣,我建议您考虑使用CUBLAS。
推荐阅读
- sql - 查询以输出每种类型的最大电影票
- python - 从张量流中的向量构造成对矩阵
- mysql - 尝试运行一个查询,向我们显示出站潜在客户与入站潜在客户的百分比
- angular - Angular6 *ngFor表与对象中的对象
- json - 跨多个字段收集给定键的值 - Apache Spark (Scala)
- mysql - 在一个表上回显两个 SQL 查询结果
- ios - 为什么在打开我的应用程序时调用它时我的数据没有取消存档
- authentication - 静默 SAML 身份验证?
- microsoft-graph-api - 需要从 AD userCertificate 属性中读取用户的证书
- html - 为什么斜体不适用于 CKEditor 5 中的日文字体