CUDA 中更快的矩阵乘法
Faster Matrix Multiplication in CUDA
目前我在cuda c上做了一个神经网络程序。因为我需要操作矩阵乘法,所以我没有为 MM 使用 CUBLAS。我为MM使用以下代码。我想知道是否有人有一些建议可以让它更快,这可能非常有帮助,因为我需要在学习过程中使用 MM 数百万次。谢谢。
这是生成文件:
# cuda root
_CUDA_ROOT_ = /usr/local/cuda
NVCC = nvcc
# include and lib paths
INCLUDES=-I${_CUDA_ROOT_}/include
LIB_PATH=-L${_CUDA_ROOT_}/lib64
# libraries to link against
LIB= -lcudart -lcublas
CU_SRC= main.cu
EXE=$(CU_SRC:.cu=)
#------------------------------
# Choose your gpu arch
SM = sm_35
all: $(EXE)
$(EXE): $(CU_SRC)
$(NVCC) -arch $(SM) $(CU_SRC) -o $(EXE) $(LIB_PATH) $(LIB)
clean:
rm -f *.o *.cu_o $(EXE)
这是MM代码:
__global__
void matrixMulti(float* A_d, float* B_d, float* C_d, int m, int k, int n)
{
__shared__ float ds_A[TILE_WIDTH][TILE_WIDTH];
__shared__ float ds_B[TILE_WIDTH][TILE_WIDTH];
int col = blockIdx.x*blockDim.x + threadIdx.x;
int row = blockIdx.y*blockDim.y + threadIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float sum = 0;
for(int t=0; t<(n-1)/TILE_WIDTH+1; t++)
{
if(row<m && t*TILE_WIDTH+tx<n)
ds_A[ty][tx] = A_d[row*n + t*TILE_WIDTH+tx];
else
ds_A[ty][tx] = 0.0;
if(t*TILE_WIDTH+ty<n && col<k)
ds_B[ty][tx] = B_d[(t*TILE_WIDTH+ty)*k + col];
else
ds_B[ty][tx] = 0.0;
__syncthreads();
for(int i=0; i<TILE_WIDTH; i++)
sum += ds_A[ty][i] * ds_B[i][tx];
__syncthreads();
}
if(row<m && col<k)
C_d[col+row*k] = sum;
}
这是代码主要部分的示例:
const int TILE_WIDTH = 32;
int main()
{
int m, k, n;
m = 10000, k = 10000, n = 10000;
float *A, *B, *C;
A = new float[m*n];
B = new float[n*k];
C = new float[m*k];
float *A_d, *B_d, *C_d;
for (int i=0; i<m*n; i++)
{
A[i] = 2;
}
for (int i=0; i<n*k; i++)
{
B[i] = 3;
}
cudaMalloc(&A_d, sizeof(float)*m*n);
cudaMalloc(&B_d, sizeof(float)*n*k);
cudaMalloc(&C_d, sizeof(float)*m*k);
cudaMemcpy(A_d, A, sizeof(float)*m*n, cudaMemcpyHostToDevice);
cudaMemcpy(B_d, B, sizeof(float)*k*n, cudaMemcpyHostToDevice);
dim3 dimGrid((k-1)/TILE_WIDTH+1, (m-1)/TILE_WIDTH+1, 1);
dim3 dimBlock(TILE_WIDTH, TILE_WIDTH, 1);
matrixMulti<<<dimGrid,dimBlock>>>(A_d, B_d, C_d, m, k, n);
cudaMemcpy(C, C_d, sizeof(float)*m*k, cudaMemcpyDeviceToHost);
return 0;
}
首先,确定这就是您想要做的。如果不描述您想要执行的操作,就很难对此发表评论,但请注意矩阵乘法是 n 次方操作。如果您的操作不那么复杂,那么使用 cuBLAS 很可能会做得更好。
这是为什么? cuBLAS 可能比您将要编写的任何东西都快,并且由于它将遵循新的 GPU 架构,因此更易于维护。 GEMM 之类的最佳实现会因架构而异,因此您现在为硬件编写的任何代码都必须针对新硬件重新优化。
现在,进入正题。您应该考虑多种技术来优化此代码:
- 每个线程计算多个输出值。这减少了共享内存的压力,因为切片数据可用于多次计算。
- 修复共享内存中的库冲突。文档应该很好地涵盖这一点。
- 矢量化共享内存加载和存储。我注意到您正在为 sm_35 编译。该架构的共享内存组每个都有 64bits/clock 的带宽。加载单个浮点数仅为 32 位,因此如果不进行矢量化,您将无法在浮点数上获得完整带宽。你应该看看 float2/float4 类型。
- 考虑双缓冲。在对另一个共享内存块进行操作时将数据加载到一个共享内存块中。这允许更有效地隐藏全局内存操作的高延迟,减少同步开销,并且往往会表现得更好。不过,它使用两倍的共享内存,因为您一次需要两个图块。
在GPU上实现矩阵乘法的论文有很多,建议你去看看。你将从这些论文中获得比你在 SO 上提出广泛问题更多的细节。
最后...您确定不想使用 cuBLAS 吗?我不会指望获得 75% 的 cuBLAS 性能,即使那也是一个挑战。
目前我在cuda c上做了一个神经网络程序。因为我需要操作矩阵乘法,所以我没有为 MM 使用 CUBLAS。我为MM使用以下代码。我想知道是否有人有一些建议可以让它更快,这可能非常有帮助,因为我需要在学习过程中使用 MM 数百万次。谢谢。 这是生成文件:
# cuda root
_CUDA_ROOT_ = /usr/local/cuda
NVCC = nvcc
# include and lib paths
INCLUDES=-I${_CUDA_ROOT_}/include
LIB_PATH=-L${_CUDA_ROOT_}/lib64
# libraries to link against
LIB= -lcudart -lcublas
CU_SRC= main.cu
EXE=$(CU_SRC:.cu=)
#------------------------------
# Choose your gpu arch
SM = sm_35
all: $(EXE)
$(EXE): $(CU_SRC)
$(NVCC) -arch $(SM) $(CU_SRC) -o $(EXE) $(LIB_PATH) $(LIB)
clean:
rm -f *.o *.cu_o $(EXE)
这是MM代码:
__global__
void matrixMulti(float* A_d, float* B_d, float* C_d, int m, int k, int n)
{
__shared__ float ds_A[TILE_WIDTH][TILE_WIDTH];
__shared__ float ds_B[TILE_WIDTH][TILE_WIDTH];
int col = blockIdx.x*blockDim.x + threadIdx.x;
int row = blockIdx.y*blockDim.y + threadIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float sum = 0;
for(int t=0; t<(n-1)/TILE_WIDTH+1; t++)
{
if(row<m && t*TILE_WIDTH+tx<n)
ds_A[ty][tx] = A_d[row*n + t*TILE_WIDTH+tx];
else
ds_A[ty][tx] = 0.0;
if(t*TILE_WIDTH+ty<n && col<k)
ds_B[ty][tx] = B_d[(t*TILE_WIDTH+ty)*k + col];
else
ds_B[ty][tx] = 0.0;
__syncthreads();
for(int i=0; i<TILE_WIDTH; i++)
sum += ds_A[ty][i] * ds_B[i][tx];
__syncthreads();
}
if(row<m && col<k)
C_d[col+row*k] = sum;
}
这是代码主要部分的示例:
const int TILE_WIDTH = 32;
int main()
{
int m, k, n;
m = 10000, k = 10000, n = 10000;
float *A, *B, *C;
A = new float[m*n];
B = new float[n*k];
C = new float[m*k];
float *A_d, *B_d, *C_d;
for (int i=0; i<m*n; i++)
{
A[i] = 2;
}
for (int i=0; i<n*k; i++)
{
B[i] = 3;
}
cudaMalloc(&A_d, sizeof(float)*m*n);
cudaMalloc(&B_d, sizeof(float)*n*k);
cudaMalloc(&C_d, sizeof(float)*m*k);
cudaMemcpy(A_d, A, sizeof(float)*m*n, cudaMemcpyHostToDevice);
cudaMemcpy(B_d, B, sizeof(float)*k*n, cudaMemcpyHostToDevice);
dim3 dimGrid((k-1)/TILE_WIDTH+1, (m-1)/TILE_WIDTH+1, 1);
dim3 dimBlock(TILE_WIDTH, TILE_WIDTH, 1);
matrixMulti<<<dimGrid,dimBlock>>>(A_d, B_d, C_d, m, k, n);
cudaMemcpy(C, C_d, sizeof(float)*m*k, cudaMemcpyDeviceToHost);
return 0;
}
首先,确定这就是您想要做的。如果不描述您想要执行的操作,就很难对此发表评论,但请注意矩阵乘法是 n 次方操作。如果您的操作不那么复杂,那么使用 cuBLAS 很可能会做得更好。
这是为什么? cuBLAS 可能比您将要编写的任何东西都快,并且由于它将遵循新的 GPU 架构,因此更易于维护。 GEMM 之类的最佳实现会因架构而异,因此您现在为硬件编写的任何代码都必须针对新硬件重新优化。
现在,进入正题。您应该考虑多种技术来优化此代码:
- 每个线程计算多个输出值。这减少了共享内存的压力,因为切片数据可用于多次计算。
- 修复共享内存中的库冲突。文档应该很好地涵盖这一点。
- 矢量化共享内存加载和存储。我注意到您正在为 sm_35 编译。该架构的共享内存组每个都有 64bits/clock 的带宽。加载单个浮点数仅为 32 位,因此如果不进行矢量化,您将无法在浮点数上获得完整带宽。你应该看看 float2/float4 类型。
- 考虑双缓冲。在对另一个共享内存块进行操作时将数据加载到一个共享内存块中。这允许更有效地隐藏全局内存操作的高延迟,减少同步开销,并且往往会表现得更好。不过,它使用两倍的共享内存,因为您一次需要两个图块。
在GPU上实现矩阵乘法的论文有很多,建议你去看看。你将从这些论文中获得比你在 SO 上提出广泛问题更多的细节。
最后...您确定不想使用 cuBLAS 吗?我不会指望获得 75% 的 cuBLAS 性能,即使那也是一个挑战。