标签:
CUDA中使用shared_memory可以加速运算,在矩阵乘法中是一个体现。
矩阵C = A * B,正常运算时我们运用 C[i,j] = A[i,:] * B[:,j] 可以计算出结果。但是在CPU上完成这个运算我们需要大量的时间,设A[m,n],B[n,k],那么C矩阵为m*k,总体,我们需要做m*n*k次乘法运算,m*(b-1)*k次加法运算,并且是串行执行,总体的复杂度为O(m*n*k) 。
矩阵类:
1 class Matrix 2 { 3 public: 4 int cols; // x 5 int rows; // y 6 float *data; //数据,一位数组 7 }
CPU上的程序,一个三层循环
for(int i =0;i< C.rows;i++) { for(int j =0;j< C.cols;j++) { float *a = A.data; float *b = B.data; for(int k=0;k<A.cols;k++) C.data[i*C.cols+j]+=a[i*A.cols+k] * b[k*B.cols+j]; } } }
我们想到用GPU加速,在CUDA上实现,我们这么写kernel:
__global__ void matrixMulKernel(const Matrix A, const Matrix B, Matrix C) { // Each thread computes one element of C // by accumulating results into Cvalue float Cvalue = 0; int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; for (int e = 0; e < A.cols; ++e) Cvalue += A.data[row * A.cols + e]* B.data[e * B.cols + col]; C.data[row * C.cols + col] = Cvalue; }
此时,计算过程是并行的,但是访问A,B矩阵时,不能同时访问,因此主要的时间花在内存读取,每个线程读取A的一行,B的一列,计算C的对应值;所以这样需要从global memory中读n次A,m次B。时间复杂度是O(m+n)次内存访问,以及k次乘法运算。
实际上还有一种办法,可以用shared memory,这里我们把A,B矩阵按照blocksize划分为子矩阵subA[blocksize][blocksize]、subB[blocksize][blocksize]。并将子矩阵设置为__shared__。 thread block内所有threads共用(可读可写)shared memory。如此一来,A只从global memory读了n/block_size次,B只读了m/block_size次;时间复杂度是O(m/block_size+n/block_size)次内存访问,以及k次乘法运算。进一步减少的时间复杂度。代码如下:
__global__ void matrixMulKernel(const float *A, const float *B, float *C,int Aw ,int Bw) { const int bs = CUDA_LG::block_size; int tx = threadIdx.x; int ty = threadIdx.y; int bx = blockIdx.x; int by = blockIdx.y; int aBlockFisrt = by * bs * Aw ; int aBlockStep = bs ; int aBlockLast = by * bs * Aw + Aw - 1 ; int bBlockFisrt = bx * bs ; int bBlockStep = bs * Bw ; float subC=0; for(int a = aBlockFisrt,int b = bBlockFisrt; a <= aBlockLast ;a+=aBlockStep,b+=bBlockStep ) { //定义两个shared memory的子矩阵 __shared__ float subA[bs][bs]; __shared__ float subB[bs][bs]; subA[ty][tx] = A[a + ty * Aw + tx]; subB[ty][tx] = B[b + ty * Bw + tx]; __syncthreads(); for(int i = 0;i<bs;i++) { subC += subA[ty][i] * subB[i][tx]; } __syncthreads(); } C[ by*bs*Bw + bx*bs + ty * Bw +tx] = subC; }
参考sample_6.5\0_Simple\matrixMul程序。里面注释详细
参考Rachel zhang的博客CUDA学习系列之二:http://blog.csdn.net/abcjennifer/article/details/42528569
标签:
原文地址:http://www.cnblogs.com/jugg1024/p/4354672.html